<jupyter_start><jupyter_text>ECM-MPT Data Analysis  The following notebook will go through prediction analysis for the Extracellular-Matrix Multiple Particle Tracking (ECM-MPT) study of pup age in P14, P21, P28, and P35 datasets.  Table of Contents  [1.   Load Data](1.-load-data)  [2.   Analysis](2.-analysis)  [3.   Modelling](modelling)  [4.   Evaluate Results](evaluate-results) ---  1.   Load Data Loading feature dataset from ECM_CSVS folder: p14, p21, p28, data are designated with their repective age while p35 data is designated with NT_brain_2.There are 15 total videos from each age group.<jupyter_code># libraries used
import boto3
import diff_classifier.aws as aws
import pandas as pd
import seaborn as sn
import numpy as np
import matplotlib.pyplot as pl


from os import listdir, getcwd, chdir
from os.path import isfile, join
import os

from matplotlib import colors as plt_colors

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn import metrics
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from sklearn import preprocessing

import xgboost as xgb
# from xgboost import cv

import shap
!pwd
dataset_path = workbookDir + '/diff_predictor/data/raw_data_age/'
filelist = [f for f in listdir(dataset_path) if isfile(join(dataset_path, f)) and ('features_' in f)]
filelist
fstats_tot = None
video_num = 0
for filename in filelist:
    try:
        fstats = pd.read_csv(dataset_path + filename, encoding = "ISO-8859-1", index_col='Unnamed: 0')
        print('{} size: {}'.format(filename, fstats.shape))
        if 'P14' in filename:
            fstats['age'] = pd.Series(fstats.shape[0]*[14], index=fstats.index)
        elif 'P21' in filename:
            fstats['age'] = pd.Series(fstats.shape[0]*[21], index=fstats.index)
        elif 'P28' in filename:
            fstats['age'] = pd.Series(fstats.shape[0]*[28], index=fstats.index)
        elif 'NT_brain_2' in filename:
            fstats['age'] = pd.Series(fstats.shape[0]*[35], index=fstats.index)
        elif 'P70' in filename:
            fstats['age'] = pd.Series(fstats.shape[0]*[70], index=fstats.index)
        else:
            print('Error, no target')
            continue
        fstats['Video Number'] = pd.Series(fstats.shape[0]*[video_num], index=fstats.index)
        if fstats_tot is None:
            fstats_tot = fstats
        else:
            fstats_tot = fstats_tot.append(fstats, ignore_index=True)
        video_num += 1
    except Exception:
            print('Skipped!: {}'.format(filename))
fstats_tot.head()<jupyter_output><empty_output><jupyter_text>2.   Analysis The following columns are present within the downloaded datasets:<jupyter_code>fstats_tot.columns<jupyter_output><empty_output><jupyter_text>Many of these features are not useful for prediction or have data which may negatively impact classification. The following features and the target feature are defined in the following cell. We also remove any datapoints that are empty or infinite:<jupyter_code>fstats_tot
features = [
    'alpha', # Fitted anomalous diffusion alpha exponenet
    'D_fit', # Fitted anomalous diffusion coefficient
    'kurtosis', # Kurtosis of track
    'asymmetry1', # Asymmetry of trajecory (0 for circular symmetric, 1 for linear)
    'asymmetry2', # Ratio of the smaller to larger principal radius of gyration
    'asymmetry3', # An asymmetric feature that accnts for non-cylindrically symmetric pt distributions
    'AR', # Aspect ratio of long and short side of trajectory's minimum bounding rectangle
    'elongation', # Est. of amount of extension of trajectory from centroid
    'boundedness', # How much a particle with Deff is restricted by a circular confinement of radius r
    'fractal_dim', # Measure of how complicated a self similar figure is
    'trappedness', # Probability that a particle with Deff is trapped in a region
    'efficiency', # Ratio of squared net displacement to the sum of squared step lengths
    'straightness', # Ratio of net displacement to the sum of squared step lengths
    'MSD_ratio', # MSD ratio of the track
#     'frames', # Number of frames the track spans
    'Deff1', # Effective diffusion coefficient at 0.33 s
    'Deff2', # Effective diffusion coefficient at 3.3 s
    'angle_mean', # Mean turning angle which is counterclockwise angle from one frame point to another
    'angle_mag_mean', # Magnitude of the turning angle mean
    'angle_var', # Variance of the turning angle
    'dist_tot', # Total distance of the trajectory
    'dist_net', # Net distance from first point to last point
    'progression', # Ratio of the net distance traveled and the total distance
    'Mean alpha', 
    'Mean D_fit', 
    'Mean kurtosis', 
    'Mean asymmetry1', 
    'Mean asymmetry2',
    'Mean asymmetry3', 
    'Mean AR',
    'Mean elongation', 
    'Mean boundedness',
    'Mean fractal_dim', 
    'Mean trappedness', 
    'Mean efficiency',
    'Mean straightness', 
    'Mean MSD_ratio', 
    'Mean Deff1', 
    'Mean Deff2',
    ]

target = 'age'           # prediction target (y)

ecm = fstats_tot
ecm = ecm[~ecm.isin([np.nan, np.inf, -np.inf]).any(1)]       # Removing nan and inf data points
# Showing a piece of our data:
ecm.head()<jupyter_output><empty_output><jupyter_text>Before prediction, it is required to balance data. As shown, The current dataset is highly imbalance with most datapoints belonging to P21 and P35 categories. The dataset is reduced using random sampling of each target category.<jupyter_code>def balance_data(df, target, **kwargs):
    if 'random_state' not in kwargs:
        random_state = 1
    else:
        random_state = kwargs['random_state']
    df_target = []
    bal_df = []
    for name in df[target].unique():
        df_target.append((name, df[df[target] == name]))
    print(f"Ratio before data balance ({':'.join([str(i[0]) for i in df_target])}) = {':'.join([str(len(i[1])) for i in df_target])}")
    for i in range(len(df_target)):
        ratio = min([len(i[1]) for i in df_target])/len(df_target[i][1])
        bal_df.append(df_target[i][1].sample(frac=ratio, random_state=random_state))
    print(f"Ratio after balance ({':'.join([str(i[0]) for i in df_target])}) = {':'.join([str(len(i)) for i in bal_df])}")
    return pd.concat(bal_df)
bal_ecm = balance_data(ecm, target, random_state=1)<jupyter_output>Ratio before data balance (35:14:21:28:70) = 17169:6416:20665:6194:14008
Ratio after balance (35:14:21:28:70) = 6194:6194:6194:6194:6194<jupyter_text>3.   Modelling The model used for this study is an extreme gradient boosting (XGBoost) decision tree which is a boosted decision tree. This model was used due to its past results within competitions and research. Due to the use of statistical surroundings in our feature analysis, binning is required in order to avoid data leakage between training/testing. The followingcode will implement binning and a checkerboard implementation to select certain bins for the training dataset.<jupyter_code># def bin_data(data, ):
#     pass

resolution = 128
assert not 2048%resolution and resolution >= 128, "resolution needs to be a factor of 2048 and > 128"
bins = list(range(0, 2048+1, resolution))
bin_labels = [int(i/resolution) for i in bins][:-1]
bal_ecm['binx'] = pd.cut(bal_ecm.X, bins, labels=bin_labels, include_lowest=True)
bal_ecm['biny'] = pd.cut(bal_ecm.Y, bins, labels=bin_labels, include_lowest=True)
bal_ecm['bins'] = (len(bins)-1)*bal_ecm['binx'].astype(np.int32) + bal_ecm['biny'].astype(np.int32)
bal_ecm = bal_ecm[np.isfinite(bal_ecm['bins'])]
bal_ecm['bins'] = bal_ecm['bins'].astype(int)

# cols = bal_ecm.columns.tolist()
# cols = cols[-3:] + cols[:-3]
# bal_ecm = bal_ecm[cols]
# Regular split

seed = 1234
np.random.seed(seed)
train_split = 0.9
test_split = 0.5

le = preprocessing.LabelEncoder()
bal_ecm['encoded_target'] = le.fit_transform(bal_ecm[target])

training_bins = np.random.choice(bal_ecm.bins.unique(), int(len(bal_ecm.bins.unique())*train_split), replace=False)

X_train = bal_ecm[bal_ecm.bins.isin(training_bins)]
X_test_val = bal_ecm[~bal_ecm.bins.isin(training_bins)]
X_val, X_test = train_test_split(X_test_val, test_size=test_split, random_state=seed)

y_train = X_train['encoded_target']
y_test = X_test['encoded_target']
y_val = X_val['encoded_target']

dtrain = xgb.DMatrix(X_train[features], label=y_train)
dtest = xgb.DMatrix(X_test[features], label=y_test)
dval = xgb.DMatrix(X_val[features], label=y_val)
#Check lengths of datasets:
def get_lengths(df, X_train, X_test, X_val=None):
    print(f'Tot before split: {len(df)}')
    print(f'Training: {len(X_train)} ({len(X_train)/len(bal_ecm):.3f}%)')
    print(f'Testing: {len(X_test)} ({len(X_test)/len(bal_ecm):.3f}%)')
    try:
        print(f'Evaluation: {len(X_val)} ({len(X_val)/len(bal_ecm):.3f}%)')
    except:
        pass
get_lengths(bal_ecm, X_train, X_test, X_val)<jupyter_output>Tot before split: 30970
Training: 27719 (0.895%)
Testing: 1626 (0.053%)
Evaluation: 1625 (0.052%)<jupyter_text>Model parameters are based on the best possible XGBoost parameters to minimize logloss error.<jupyter_code>param = {'max_depth': 10,
         'eta': 0.005,
         'min_child_weight': 8,
         'verbosity': 0,
         'objective': 'multi:softprob',
         'num_class': 5,
         'silent': 'True',
         'gamma': 5,
         'subsample': 0.6,
         'colsample_bytree': 0.6,
         'eval_metric': "mlogloss",
         # GPU integration will cut time in ~half:
         'gpu_id' : 0,
         'tree_method': 'gpu_hist',
         'predictor': 'gpu_predictor'
         }
# from xgboost import XGBClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
# model = XGBClassifier()
features
from xgboost.training import CVPack
from xgboost import callback
from xgboost.core import CallbackEnv
from xgboost.core import EarlyStopException

def cv(params, X_train, y_train, features=None, num_boost_round=20, nfold=3, stratified=False, folds=None,
       metrics=(), obj=None, feval=None, maximize=False, early_stopping_rounds=None,
       fpreproc=None, as_pandas=True, verbose_eval=None, show_stdv=True,
       seed=0, callbacks=None, shuffle=True):
    # pylint: disable = invalid-name
    """Cross-validation with given parameters.

    Parameters
    ----------
    params : dict
        Booster params.
    dtrain : DMatrix
        Data to be trained.
    num_boost_round : int
        Number of boosting iterations.
    nfold : int
        Number of folds in CV.
    stratified : bool
        Perform stratified sampling.
    folds : a KFold or StratifiedKFold instance or list of fold indices
        Sklearn KFolds or StratifiedKFolds object.
        Alternatively may explicitly pass sample indices for each fold.
        For ``n`` folds, **folds** should be a length ``n`` list of tuples.
        Each tuple is ``(in,out)`` where ``in`` is a list of indices to be used
        as the training samples for the ``n`` th fold and ``out`` is a list of
        indices to be used as the testing samples for the ``n`` th fold.
    metrics : string or list of strings
        Evaluation metrics to be watched in CV.
    obj : function
        Custom objective function.
    feval : function
        Custom evaluation function.
    maximize : bool
        Whether to maximize feval.
    early_stopping_rounds: int
        Activates early stopping. Cross-Validation metric (average of validation
        metric computed over CV folds) needs to improve at least once in
        every **early_stopping_rounds** round(s) to continue training.
        The last entry in the evaluation history will represent the best iteration.
        If there's more than one metric in the **eval_metric** parameter given in
        **params**, the last metric will be used for early stopping.
    fpreproc : function
        Preprocessing function that takes (dtrain, dtest, param) and returns
        transformed versions of those.
    as_pandas : bool, default True
        Return pd.DataFrame when pandas is installed.
        If False or pandas is not installed, return np.ndarray
    verbose_eval : bool, int, or None, default None
        Whether to display the progress. If None, progress will be displayed
        when np.ndarray is returned. If True, progress will be displayed at
        boosting stage. If an integer is given, progress will be displayed
        at every given `verbose_eval` boosting stage.
    show_stdv : bool, default True
        Whether to display the standard deviation in progress.
        Results are not affected, and always contains std.
    seed : int
        Seed used to generate the folds (passed to numpy.random.seed).
    callbacks : list of callback functions
        List of callback functions that are applied at end of each iteration.
        It is possible to use predefined callbacks by using
        :ref:`Callback API <callback_api>`.
        Example:

        .. code-block:: python

            [xgb.callback.reset_learning_rate(custom_rates)]
    shuffle : bool
        Shuffle data before creating folds.

    Returns
    -------
    evaluation history : list(string)
    """
    if stratified is True and not SKLEARN_INSTALLED:
        raise XGBoostError('sklearn needs to be installed in order to use stratified cv')
    if isinstance(metrics, str):
        metrics = [metrics]
    if not features:
        features = X_train.columns
    if isinstance(params, list):
        _metrics = [x[1] for x in params if x[0] == 'eval_metric']
        params = dict(params)
        if 'eval_metric' in params:
            params['eval_metric'] = _metrics
    else:
        params = dict((k, v) for k, v in params.items())
    if (not metrics) and 'eval_metric' in params:
        if isinstance(params['eval_metric'], list):
            metrics = params['eval_metric']
        else:
            metrics = [params['eval_metric']]
    params.pop("eval_metric", None)
    results = {}
    # create folds in data
    cvfolds, wt_list = mknfold(X_train, y_train, nfold, params, metrics, features)
    
    # setup callbacks
    callbacks = [] if callbacks is None else callbacks
    if early_stopping_rounds is not None:
        callbacks.append(callback.early_stop(early_stopping_rounds,
                                             maximize=maximize,
                                             verbose=False))
    if isinstance(verbose_eval, bool) and verbose_eval:
        callbacks.append(callback.print_evaluation(show_stdv=show_stdv))
    elif isinstance(verbose_eval, int):
        callbacks.append(callback.print_evaluation(verbose_eval, show_stdv=show_stdv))
    callbacks_before_iter = [
        cb for cb in callbacks if
        cb.__dict__.get('before_iteration', False)]
    callbacks_after_iter = [
        cb for cb in callbacks if
        not cb.__dict__.get('before_iteration', False)]
    for i in range(num_boost_round):
        for cb in callbacks_before_iter:
            cb(CallbackEnv(model=None,
                           cvfolds=cvfolds,
                           iteration=i,
                           begin_iteration=0,
                           end_iteration=num_boost_round,
                           rank=0,
                           evaluation_result_list=None))
        for fold in cvfolds:
            fold.update(i, obj)
        res = aggcv([f.eval(i, feval) for f in cvfolds], wt_list)
        for key, mean, std in res:
            if key + '-mean' not in results:
                results[key + '-mean'] = []
            if key + '-std' not in results:
                results[key + '-std'] = []
            results[key + '-mean'].append(mean)
            results[key + '-std'].append(std)
        try:
            for cb in callbacks_after_iter:
                cb(CallbackEnv(model=None,
                               cvfolds=cvfolds,
                               iteration=i,
                               begin_iteration=0,
                               end_iteration=num_boost_round,
                               rank=0,
                               evaluation_result_list=res))
        except EarlyStopException as e:
            for k in results:
                results[k] = results[k][:(e.best_iteration + 1)]
            break
    
    if as_pandas:
        try:
            import pandas as pd
            results = pd.DataFrame.from_dict(results)
        except ImportError:
            pass
    return results
def bin_fold(X_train, nfold):
    bin_list = [X_train[X_train['bins'] == i_bin].index.to_numpy() for i_bin in X_train.bins.unique()]
    bin_list = sorted(bin_list, key=len)
    i = 0
    while(len(bin_list) > nfold):
        if (i >= len(bin_list)-1):
            i = 0
        bin_list[i] = np.concatenate([bin_list[i], bin_list.pop()])
        i += 1
    wt_list = [len(i)/sum(len(s) for s in bin_list) for i in bin_list]
    return bin_list, wt_list

def mknfold(X_train, y_train, nfold, param, evals=(), features=None):
    if not features:
        features = X_train.columns
    dall = xgb.DMatrix(X_train[features], label=y_train)
    out_idset, wt_list = bin_fold(X_train, nfold)
    in_idset = [np.concatenate([out_idset[i] for i in range(nfold) if k != i]) for k in range(nfold)]
    evals = list(evals)
    ret = []
    for k in range(nfold):
        # perform the slicing using the indexes determined by the above methods
        x_train_snip = X_train.loc[in_idset[k]][features]
        y_train_snip = X_train.loc[in_idset[k]]['encoded_target']        
        x_test_snip = X_train.loc[out_idset[k]][features]
        y_test_snip = X_train.loc[out_idset[k]]['encoded_target']
        dtrain = xgb.DMatrix(x_train_snip, label=y_train_snip)
        dtest = xgb.DMatrix(x_test_snip, label=y_test_snip)
        tparam = param
        plst = list(tparam.items()) + [('eval_metric', itm) for itm in evals]
        ret.append(CVPack(dtrain, dtest, plst))
    return ret, wt_list
from xgboost.core import STRING_TYPES

def aggcv(rlist, wt_list):
    # pylint: disable=invalid-name
    """
    Aggregate cross-validation results.

    If verbose_eval is true, progress is displayed in every call. If
    verbose_eval is an integer, progress will only be displayed every
    `verbose_eval` trees, tracked via trial.
    """
    cvmap = {}
    idx = rlist[0].split()[0]
    for line in rlist:
        arr = line.split()
        assert idx == arr[0]
        for metric_idx, it in enumerate(arr[1:]):
            if not isinstance(it, STRING_TYPES):
                it = it.decode()
            k, v = it.split(':')
            if (metric_idx, k) not in cvmap:
                cvmap[(metric_idx, k)] = []
            cvmap[(metric_idx, k)].append(float(v))
    msg = idx
    results = []
    for (metric_idx, k), v in sorted(cvmap.items(), key=lambda x: x[0][0]):
        v = np.array(v)
        if not isinstance(msg, STRING_TYPES):
            msg = msg.decode()
        mean = np.average(v, weights=wt_list)
        std = np.average((v-mean)**2, weights=wt_list)
        results.extend([(k, mean, std)])
    return results
cv_results = cv(param, X_train, y_train, features, num_boost_round=200, nfold=5, early_stopping_rounds=3, metrics={'mlogloss'})
cv_results
seed = 1234
np.random.seed(seed)
param
import operator
import numpy as np

def xgb_paramsearch(X_train, y_train, features, init_params, nfold=5, num_boost_round=2000, early_stopping_rounds=3, **kwargs):
    params = {**init_params}
    if 'use_gpu' in kwargs and kwargs['use_gpu']:
        # GPU integration will cut cv time in ~half:
        params.update({'gpu_id' : 0,
                       'tree_method': 'gpu_hist',
                       'predictor': 'gpu_predictor'})
    if 'metrics' not in kwargs:
        metrics = {params['eval_metric']}
        print('test')
    else:
        metrics = kwargs['metrics']
        metrics.add(params['eval_metric'])
        print(metrics)
    if params['eval_metric'] in ['map', 'auc', 'aucpr']:
        eval_f = operator.gt
    else: 
        eval_f = operator.lt
    if 'early_break' not in kwargs:
        early_break = 5
    else: 
        early_break = kwargs['early_break']
    if 'thresh' not in kwargs:
        thresh = 0.01
    else: 
        thresh = kwargs['thresh']
    if 'seed' not in kwargs:
        seed = 1111
    else: 
        seed = kwargs['seed']
    best_param = params
    best_model = cv(params, 
                    X_train, 
                    y_train, 
                    features, 
                    nfold=nfold, 
                    num_boost_round=num_boost_round, 
                    early_stopping_rounds=early_stopping_rounds, 
                    metrics=metrics)
    best_eval = best_model[f"test-{params['eval_metric']}-mean"].min()
    best_boost_rounds = best_model[f"test-{params['eval_metric']}-mean"].idxmin()
    
    def _gs_helper(var1n, var2n, best_model, best_param, best_eval, best_boost_rounds):
        local_param = {**best_param}
        for var1, var2 in gs_params:
            print(f"Using CV with {var1n}={{{var1}}}, {var2n}={{{var2}}}")
            local_param[var1n] = var1
            local_param[var2n] = var2
            cv_model = cv(local_param, 
                          X_train, 
                          y_train, 
                          features, 
                          nfold=nfold, 
                          num_boost_round= num_boost_round, 
                          early_stopping_rounds=early_stopping_rounds, 
                          metrics=metrics)
            cv_eval = cv_model[f"test-{local_param['eval_metric']}-mean"].min()
            boost_rounds = cv_model[f"test-{local_param['eval_metric']}-mean"].idxmin()
            if(eval_f(cv_eval, best_eval)):
                best_model = cv_model
                best_param[var1n] = var1
                best_param[var2n] = var2
                best_eval = cv_eval
                best_boost_rounds = boost_rounds
                print(f"New best param found: "
                      f"{local_param['eval_metric']} = {{{best_eval}}}, "
                      f"boost_rounds = {{{best_boost_rounds}}}")
        return best_model, best_param, best_eval, best_boost_rounds
    
    while(early_break >= 0):
        np.random.seed(seed)
        best_eval_init = best_eval
        gs_params = {
            (subsample, colsample)
            for subsample in np.random.choice([i/10. for i in range(5,11)], 3)
            for colsample in np.random.choice([i/10. for i in range(5,11)], 3)
        }
        best_model, best_param, best_eval, best_boost_rounds = _gs_helper('subsample', 
                                                                          'colsample_bytree', 
                                                                          best_model, 
                                                                          best_param, 
                                                                          best_eval, 
                                                                          best_boost_rounds)
        gs_params = {
            (max_depth, min_child_weight)
            for max_depth in [10] + list(np.random.randint(1, 10, 3))
            for min_child_weight in [0, 10] + list(np.random.randint(0, 10, 3))
        }
        best_model, 
        best_param, 
        best_eval, 
        best_boost_rounds = _gs_helper('max_depth', 
                                                                          'min_child_weight', 
                                                                          best_model, 
                                                                          best_param, 
                                                                          best_eval, 
                                                                          best_boost_rounds)
        gs_params = {
            (eta, gamma)
            for eta in np.random.choice([.005, .01, .05, .1, .2, .3], 3)
            for gamma in [0] + list(np.random.choice([0.01, 0.001, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0], 3))
        }
        best_model, best_param, best_eval, best_boost_rounds = _gs_helper('eta', 
                                                                          'gamma', 
                                                                          best_model, 
                                                                          best_param, 
                                                                          best_eval, 
                                                                          best_boost_rounds)
        if (abs(best_eval_init - best_eval) < thresh):
            early_break-=1
        seed+=1
    return best_model, best_param, best_eval, best_boost_rounds

best_model, best_param, best_eval, best_boost_rounds = xgb_paramsearch(X_train, y_train, features, init_params=param, nfold=5, num_boost_round=2000, early_stopping_rounds=3, metrics={'mlogloss', 'merror'}, use_gpu='True', seed=3456)
tmp = (best_model, best_param, best_eval, best_boost_rounds)
best_param
cv_model = cv(best_param, 
              X_train, 
              y_train, 
              features=features, 
              nfold=5, 
              num_boost_round=best_boost_rounds,
              early_stopping_rounds=3, 
              metrics={'mlogloss'}) 
cv_model[f"test-mlogloss-mean"].min()
best_boost_rounds
# evals = [(dtrain, 'train'), (dval, 'eval')]
# num_round = best_boost_rounds
# bst = xgb.train(param, dtrain, num_round, evals, )


# ######
# label = dtest.get_label()
# ypred1 = booster.predict(dtest)
# # by default, we predict using all the trees
# pred = [np.where(x == np.max(x))[0][0] for x in ypred1]
# print("Accuracy:",metrics.accuracy_score(label, pred))
def train(param, dtrain, dtest, dval=None, evals=[(dtrain, 'train')], num_round=40):
    if dval is not None and (dval, 'eval') not in evals:
        evals += [(dval, 'eval')]
    model = xgb.train(param, dtrain, num_round, evals, )
    true_label = dtest.get_label()
    ypred = model.predict(dtest)
    preds = [np.where(x == np.max(x))[0][0] for x in ypred]
    acc = metrics.accuracy_score(true_label, preds)
    print("Accuracy:",acc)
    return model, acc
booster, acc = train(best_param, dtrain, dtest, dval, evals=[(dtrain, 'train'), (dval, 'eval')], num_round=32)

# To load in a model, use this
bst = xgb.XGBClassifier()
booster = xgb.Booster({'nthread':4})
from datetime import date
import json
# bst.save_model(f'model_xgboost_P14_P21_P28_P32_50-50-split_{str(date.today())}')
# with open(f'config_xgboost_P14_P21_P28_P32_50-50-split_{str(date.today())}', 'w', encoding='utf-8') as f:
#     json.dump(bst.save_config(), f, ensure_ascii=False, indent=4)
    
def save(model, folder, filename):
    filename = 'model_' + filename.strip("/")
    model_file = '/'.join([folder, filename])
    config_file = model_file.replace('model_', 'config_')
    model.save_model(model_file)
    with open(config_file, 'w', encoding='utf-8') as f:
        json.dump(model.save_config(), f, ensure_ascii=False, indent=4)
    f.close()

save(booster, 'diff_predictor/data/saved_models/', f'xgboost_age_P14_P21_P28_P35_P70_9050split_{str(date.today())}')
from datetime import date
import json
bst_model = f'xgboost_age_P14_P21_P28_P35_P70_9050split_2020-11-04'

def load(folder, filename):
    booster = xgb.Booster({'nthread':4})
    filename = 'model_' + filename
    model_file = '/'.join([folder, filename])
    config_file = model_file.replace('model_', 'config_')
    booster.load_model(model_file)
    with open(config_file, 'r', encoding='utf-8') as f:
        config = f.read()
        config = json.loads(config)
    booster.load_config(config)
    metadata = json.loads(config)
    return booster, metadata

booster, metadata = load('diff_predictor/data/saved_models', bst_model)
def search_nested_dict(value, keyword):
    '''
    Function which searches a nested dictionary
    recurrently for a keyword.
    Parameters
    ----------
    value : dict
        dictionary to search through
    keyword : str
        keyword to find withing dictionary
    Returns
    -------
    result
        resulting value of search
    '''
    result = None
    if keyword in value.keys() and not isinstance(value[keyword], dict):
        return(value[keyword])
    else:
        for key in value.keys():
            if isinstance(value[key], dict):
                result = search_nested_dict(value[key], keyword)
            if result is not None:
                break
    return result
search_nested_dict(metadata, 'objective')
def is_numeric(value):
    try:
        return float(value)
    except Exception:
        return False
def get_params(params, metadata):
    '''
    Gets parameters from a loaded configuration JSON file
    '''
    for key in params.keys():
        if params[key] is None:
            value = search_nested_dict(metadata, key)
            if is_numeric(value):
                value = float(value)
                if value % 1 == 0:
                    params[key] = int(value)
                else: 
                    params[key] = float(value)
            else:
                params[key] = value
    return params

param = {'max_depth': None,
         'eta': None,
         'min_child_weight': None,
         'verbosity': None,
         'objective': None,
         'num_class': None,
         'silent': 'True',
         'gamma': None,
         'subsample': None,
         'colsample_bytree': None,
         'eval_metric': 'mlogloss'}
param = get_params(param, metadata)
param
# Chack loaded model for accuracy
ypred1 = booster.predict(dtest)
# by default, we predict using all the trees
pred = [np.where(x == np.max(x))[0][0] for x in ypred1]

print("Accuracy:",metrics.accuracy_score(y_test, pred))
def perf_meas(y_actual, y_pred, cls):
    TP = 0
    FP = 0
    TN = 0
    FN = 0
    for i in range(len(y_actual)): 
        if (y_actual[i]==y_pred[i]) and (y_pred[i]==cls):
           TP += 1
        if (y_pred[i]==cls) and (y_actual[i]!=y_pred[i]):
           FP += 1
        if (y_actual[i]==y_pred[i]) and (y_pred[i]!=cls):
           TN += 1
        if (y_pred[i]!=cls) and (y_actual[i]!=y_pred[i]):
           FN += 1
    return(TP, FP, TN, FN)
perf_meas(list(y_test), pred, 3)
bst = booster
# Need to run this on models to run Shap plots
model_bytearray = bst.save_raw()[4:]
def myfun(self=None):
    return model_bytearray
bst.save_raw = myfun
results = X_test[features]
results['predicted'] = pred
results['actual'] = y_test<jupyter_output>A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy<jupyter_text>4.   Evaluate Results<jupyter_code># print('0 == {}'.format(le.inverse_transform([0]))) # Checks to make sure encoding is right for each class
# print('1 == {}'.format(le.inverse_transform([1])))
# print('2 == {}'.format(le.inverse_transform([2])))
# print('3 == {}'.format(le.inverse_transform([3])))
# print('4 == {}'.format(le.inverse_transform([4])))
class_names = ['P14', 'P21', 'P28', 'P35', 'P70']
class_results = classification_report(y_test, pred, digits=4, target_names = ['P14', 'P21', 'P28', 'P35', 'P70'])
print(str(class_results))
# import matplotlib.pyplot as plt

# model = booster

# ax = plt.subplot(211)
# ax.plot(model.history.history['loss'], label='train')
# ax.plot(model.history.history['val_loss'], label='test')
# ax.legend()
# ax.set_ylabel('Loss')
# ax.set_xlabel('Number of Boosting Rounds')

# ax2 = plt.subplot(212)
# ax2.plot(model.history.history['accuracy'], label='train')
# ax2.plot(model.history.history['val_accuracy'], label='test')
# ax2.legend()
# ax2.set_ylabel('Accuracy')
# ax2.set_xlabel('Number of Boosting Rounds')
# py.show()
confusion_matrix(y_test, pred)
pl.figure(figsize=(12,10))
cm_array = confusion_matrix(y_test, pred)
df_cm = pd.DataFrame(cm_array, index = class_names, columns = class_names)

sn.set(font_scale=1.4) # for label size
ax = sn.heatmap(df_cm, annot=True, annot_kws={"size": 16}, cmap="YlGnBu")
ax.set(xlabel='Predicted', ylabel='Actual')

pl.show()
explainer = shap.TreeExplainer(bst)
shap_values = explainer.shap_values(X_test[features])
%matplotlib inline
import matplotlib.pyplot as pl
cP14 = '#808080'
cP21 = '#daa520'
cP28 = '#6495ed'
cP35 = '#9400d3'
cP70 = '#006400'
colors = [cP14, cP21, cP28, cP35, cP70]
class_inds = np.argsort([-np.abs(shap_values[i]).mean() for i in range(len(shap_values))])
cmap = plt_colors.ListedColormap(np.array(colors)[class_inds])
# sn.reset_orig() # Reset matplot lib to no longer use seaborn
figsize = (7.5, 5)
fig = pl.figure(figsize=figsize)
ax = fig.gca()
shap.summary_plot(shap_values, X_test[features], class_names=np.array(class_names), max_display=15, title='Total SHAP Values', color=cmap, show=False)
handles, labels = ax.get_legend_handles_labels()
ax.set_xlabel(r'Average Impact on Model Output Magnitude (Mean(|SHAP|))')
handles, labels = ax.get_legend_handles_labels()
labels, handles = zip(*sorted(zip(labels, handles), key=lambda t: t[0]))
ax.legend(handles, labels)
labels
# This is optional, but right units of diffusion are in pixels^2/frame.
# To convert, There are 0.07 um/pixel and 30.3 frames/s in videos.
# *This may be different depending on zoom of camera, make sure to ask 
# whoever took video.
X_test['Mean Deff1'] = X_test['Mean Deff1']*(0.07**2)/30.3
X_test['Mean D_fit'] = X_test['Mean D_fit']*(0.07**2)/30.3
X_test['progression'] = X_test['progression']*(0.07)
pl.ioff()
%matplotlib inline

figsize = (7.5, 5)
bottom = -1.5
top = 1.5
for i in range(len(shap_values)):
    fig = pl.figure(figsize=figsize)
    ax = fig.gca()
    ax.set_ylim(bottom, top)
    shap.dependence_plot("Mean fractal_dim", shap_values[i], X_test[features], interaction_index = None, color=colors[i], alpha=0.5, ax=ax, show=False)
    ax.set_xlabel(r'Mean $\alpha$')
    ax.set_ylabel('SHAP Value')
# Use for removing nan values from shap graph, this is just for cleaning up graph, prediction is still based off those vals.
# Only works for dependence_plots
df = X_test.reset_index(drop=True)
indeces = df['Mean Deff1'].index[df['Mean Deff1'].apply(pd.isna)]
df.loc[~df.index.isin(indeces)]['Mean Deff1']
not_index = [i for i in df.index if i not in indeces]
removed_values = np.array([i[not_index] for i in np.array(shap_values)]) # Use this in place of shap_values
df = df.loc[~df.index.isin(indeces)] # Use this in place of X_test
pl.ioff()
%matplotlib inline

figsize = (7.5, 5)
bottom = -3.0
top = 3.0
fig = pl.figure(figsize=figsize)
ax = fig.gca()
for i in range(len(shap_values)):
    ax.set_ylim(bottom, top)
    ax.set_xlim(-0.0005, 0.015)
    # Alternative shap values used
    shap.dependence_plot("progression", shap_values[i], X_test[features], interaction_index = None, color=colors[i], alpha=0.5, ax=ax, show=False)
    ax.set_xlabel(r'Progression')
    ax.set_ylabel('SHAP Value')
pl.ioff()
%matplotlib inline

figsize = (7.5, 5)
bottom = -3.0
top = 3.0
for i in range(len(shap_values)):
    fig = pl.figure(figsize=figsize)
    ax = fig.gca()
    ax.set_xlim(-0.0005, 0.015)
    ax.set_ylim(bottom, top)
    shap.dependence_plot("progression", shap_values[i], X_test[features], interaction_index = None, color=colors[i], alpha=0.5, ax=ax, show=False)
    ax.set_xlabel(r'Progression')
    ax.set_ylabel('SHAP Value')
figsize = (7.5, 5)
bottom = -1.5
top = 1.5
for i in range(len(shap_values)):
    fig = pl.figure(figsize=figsize)
    ax = fig.gca()
    ax.set_ylim(bottom, top)
    shap.dependence_plot("Mean fractal_dim", shap_values[i], X_test[features], interaction_index = None, color=colors[i], alpha=0.5, ax=ax, show=False)
    ax.set_xlabel('Mean Fractal Dimension')
    ax.set_ylabel('SHAP Value')
figsize = (7.5, 5)
bottom = -1.5
top = 1.5
for i in range(len(shap_values)):
    fig = pl.figure(figsize=figsize)
    ax = fig.gca()
    ax.set_ylim(bottom, top)
    shap.dependence_plot("Mean D_fit", shap_values[i], X_test[features], interaction_index = None, color=colors[i], alpha=0.5, ax=ax, show=False)
    ax.set_xlabel(r'Mean $D_{fit}$ ($\mu^2/s$)')
    ax.set_ylabel('SHAP Value')
figsize = (7.5, 5)
bottom = -1.5
top = 1.5
for i in range(len(shap_values)):
    fig = pl.figure(figsize=figsize)
    ax = fig.gca()
    ax.set_ylim(bottom, top)
    shap.dependence_plot("Mean MSD_ratio", shap_values[i], X_test[features], interaction_index = None, color=colors[i], alpha=0.5, ax=ax, show=False)
    ax.set_xlabel('Mean MSD Ratio')
    ax.set_ylabel('SHAP Value')
figsize = (7.5, 5)
bottom = -1.5
top = 1.5
for i in range(len(shap_values)):
    fig = pl.figure(figsize=figsize)
    ax = fig.gca()
    ax.set_ylim(bottom, top)
    shap.dependence_plot("Mean kurtosis", shap_values[i], X_test[features], interaction_index = None, color=colors[i], alpha=0.5, ax=ax, show=False)
    ax.set_xlabel('Mean Kurtosis')
    ax.set_ylabel('SHAP Value')
    
figsize = (7.5, 5)
fig = pl.figure(figsize=figsize)
ax = fig.gca()
shap.summary_plot(shap_values[0], X_test[features], max_display=5, class_names = class_names, title = 'SHAP Value for P14', show=False)
ax.set_title('Top 5 Features for P14 Pup')
figsize = (7.5, 5)
fig = pl.figure(figsize=figsize)
ax = fig.gca()
shap.summary_plot(shap_values[1], X_test[features], max_display=5, class_names = class_names, title = 'SHAP Value for P21', show=False)
ax.set_title('Top 5 Features for P21 Pup')
figsize = (7.5, 5)
fig = pl.figure(figsize=figsize)
ax = fig.gca()
shap.summary_plot(shap_values[2], X_test[features], max_display=5, class_names = class_names, title='SHAP Value for P28', show=False)
ax.set_title('Top 5 Features for P28 Pup')
figsize = (7.5, 5)
fig = pl.figure(figsize=figsize)
ax = fig.gca()
shap.summary_plot(shap_values[3], X_test[features], max_display=5, class_names=class_names, title='SHAP Value for P35', show=False)
ax.set_title('Top 5 Features for P35 Pup')
figsize = (7.5, 5)
fig = pl.figure(figsize=figsize)
ax = fig.gca()
shap.summary_plot(shap_values[4], X_test[features], max_display=5, class_names=class_names, title='SHAP Value for P70', show=False)
ax.set_title('Top 5 Features for P70 Pup')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11809716494176847, 'fraction_numerical': 0.01906466422846471, 'mean_word_length': 2.8179837591602297, 'pattern_counts': {'":': 1, '<': 28, '<?xml version=': 0, '>': 32, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 50, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '533439', 'license': 'MIT', 'n_tokens_mistral': 13058, 'n_tokens_neox': 12366, 'n_words': 3253, 'path': 'notebooks/XGBoost_age_analysis_P14_P21_P28_P35_P70.ipynb', 'repo_name': 'dash2927/diff_predictor', 'size': '38553'}