<jupyter_start><jupyter_text>Get the HPS model to run- get the inputs correct- run through all steps- get an accuracy estimate<jupyter_code>%matplotlib inline
import os
import sys
sys.path.append('/home/surchs/git/Proteus/')
from proteus.predic import high_confidence
#from proteus.predic import prediction

import numpy as np
import pandas as pd
import nibabel as nib
import sklearn as skl
import seaborn as sbn
import scipy.io as sio
from matplotlib import pyplot as plt
from sklearn.model_selection import StratifiedKFold
import warnings
#warnings.filterwarnings('ignore')
# Paths
root_p = '/home/surchs/sim_big/PROJECT/abide_hps/'
# Pheno
sample_p = os.path.join(root_p, 'pheno', 'psm_abide1.csv')
# Data
resid_ct_p = os.path.join(root_p, 'residual', 'ct_30mm_residual_psm.npy')
resid_fc_p = os.path.join(root_p, 'residual', 'sd_30mm_residual_psm.npy')
mask_p = os.path.join(root_p, 'mask', 'MIST_mask.nii.gz')
subtype_fc_p = os.path.join(root_p, 'subtypes', 'subtypes_fc.npz')
subtype_ct_p = os.path.join(root_p, 'subtypes', 'subtypes_ct.npz')
# Load data
sample = pd.read_csv(sample_p)
sample['DX_CODE'] = sample['DX_GROUP'].replace({'Autism':1, 'Control':0})
s_fc = np.load(subtype_fc_p)
s_ct = np.load(subtype_ct_p)
# Reshape FC to add to the table
fc_weights = np.reshape(s_fc['weights'], (370, 5*20))
fc_cols = ['fc_n{}_s{}'.format(nid+1, sid+1) for sid in range(5) for nid in range(20)]
# Same for CT
ct_weights = s_ct['weights']
ct_cols = ['ct_s{}'.format(sid+1) for sid in range(5)]
# Combine both
weights = np.concatenate((ct_weights, fc_weights),1)
cols = ct_cols + fc_cols
# Into a pandas DF
w_data = pd.DataFrame(data=weights, columns=cols)
# Combine both for the full dataset
dataset = sample.join(w_data)
# Select the features
scaler = skl.preprocessing.StandardScaler()
# Add BV to the subtype weights
col_features = ['BV'] + cols

# Build features
x_ = dataset.loc[:, col_features]
# Take the numeric diagnosis code
y_ = dataset.loc[:, ['DX_CODE']].values.squeeze()
# Normalize
x_ = scaler.fit_transform(x_)
# Prep store
store = {key:list() for key in ['accuracy_overall',
                                'accuracy_asd',
                                'precision_asd',
                                'precision_tdc',
                                'recall_asd',
                                'recall_tdc',
                                'f1_asd',
                                'f1_tdc'
                               ]}

# Run the model and see where that gets us
skf = StratifiedKFold(n_splits=3)
for train_index, val_index in skf.split(x_,y_):
    X_training, X_val = x_[train_index], x_[val_index]
    y_training, y_val = y_[train_index], y_[val_index]
    
    hpc = high_confidence.TwoStagesPrediction(
        n_iter=500,
        shuffle_test_split=0.5,
        min_gamma=.90,
        gamma=0.95,
        thresh_ratio=0.1,
        verbose=False)
    
    hpc.fit(X_training, X_training, y_training)
    _, dic_results = hpc.predict(X_val, X_val)
    acc = skl.metrics.accuracy_score(y_val, (dic_results['s1df'][:,0]>0).astype(float))
    store['accuracy_overall'].append(acc)
    
    # Get the guys we think are ASD
    pos_mask = (dic_results['s2df'][:,1]>0)
    acc_s2 = skl.metrics.accuracy_score(y_val[pos_mask], (dic_results['s1df'][:,0]>0).astype(float)[pos_mask])
    store['accuracy_asd'].append(acc_s2)
    
    print('Classifying TDC vs ASD...')
    print((dic_results['s1df'][:,0]>0).astype(float))
    
    y_pred = (dic_results['s1df'][:,0]>0).astype(float)
    # Compute precision of the classifier for ASD
    asd_p = skl.metrics.precision_score(y_val, y_pred, pos_label=1)
    store['precision_asd'].append(asd_p)
    # Compute precision of the classifier for TDC
    tdc_p = skl.metrics.precision_score(y_val, y_pred, pos_label=0)
    store['precision_tdc'].append(tdc_p)
    # Recall Ratio of ASD label
    asd_r = skl.metrics.recall_score(y_val, y_pred, pos_label=1)
    store['recall_asd'].append(asd_r)
    # Recall Ratio of TDC label
    tdc_r = skl.metrics.recall_score(y_val, y_pred, pos_label=0)
    store['recall_tdc'].append(tdc_r)
    # F1 Ratio of ASD label
    asd_f = skl.metrics.f1_score(y_val, y_pred, pos_label=1)
    store['f1_asd'].append(asd_f)
    # Recall Ratio of TDC label
    tdc_f = skl.metrics.f1_score(y_val, y_pred, pos_label=0)
    store['f1_tdc'].append(tdc_f)
#print(scores_ad_cn)
print('Mean stage 1 validation accuracy: ',np.mean(store['accuracy_overall']))
#print(scores_s2)
print('Mean stage 2 accuracy: ', np.mean(store['accuracy_asd']))
#print(ad_precision)
print('Mean precision for ASD: ',np.mean(store['precision_asd']))
#print(cn_precision)
print('Mean precision for TDC: ',np.mean(store['precision_tdc']))
#print(ad_recall)
print('Mean recall for ASD: ',np.mean(store['recall_asd']))
#print(cn_recall)
print('Mean recall for TDC: ',np.mean(store['recall_tdc']))
#print(ad_f1_score)
print('Mean f1 score for ASD: ',np.mean(store['f1_asd']))
#print(cn_f1_score)
print('Mean f1 score for TDC: ',np.mean(store['f1_tdc']))
# Run the model and see where that gets us
skf = StratifiedKFold(n_splits=3)
for train_index, val_index in skf.split(x_,y_):
    X_training, X_val = x_[train_index], x_[val_index]
    y_training, y_val = y_[train_index], y_[val_index]
    
    hpc = high_confidence.TwoStagesPrediction(
        n_iter=500,
        shuffle_test_split=0.5,
        min_gamma=.90,
        gamma=0.95,
        thresh_ratio=0.1,
        verbose=False)
    
    hpc.fit(X_training, X_training, y_training)
    _, dic_results = hpc.predict(X_val, X_val)
    acc = skl.metrics.accuracy_score(y_val, (dic_results['s1df'][:,0]>0).astype(float))
    store['accuracy_overall'].append(acc)
    
    # Get the guys we think are ASD
    pos_mask = (dic_results['s2df'][:,1]>0)
    acc_s2 = skl.metrics.accuracy_score(y_val[pos_mask], (dic_results['s1df'][:,0]>0).astype(float)[pos_mask])
    store['accuracy_asd'].append(acc_s2)
    
    print('Classifying TDC vs ASD...')
    print((dic_results['s1df'][:,0]>0).astype(float))
    array_results, dic_results = hpc.predict(X_val, X_val)
    y_pred = (dic_results['s1df'][:,0]>0).astype(float)
    lr_decision = dic_results['s2df'][:,1]
    predic_stats(y_val, y_pred, lr_decision)<jupyter_output>Stage 1
Stage 2
Classifying TDC vs ASD...
[ 0.  0.  0.  0.  1.  1.  0.  1.  0.  1.  0.  0.  0.  1.  1.  0.  1.  1.
  1.  0.  1.  1.  0.  0.  1.  0.  0.  0.  1.  1.  0.  1.  0.  0.  1.  1.
  1.  0.  0.  0.  0.  1.  0.  1.  1.  1.  1.  0.  1.  0.  1.  1.  0.  0.
  1.  0.  1.  1.  0.  0.  0.  1.  1.  0.  1.  0.  1.  0.  0.  0.  1.  0.
  0.  1.  0.  1.  0.  1.  0.  0.  1.  1.  0.  1.  1.  0.  0.  1.  0.  1.
  1.  0.  0.  0.  1.  1.  0.  0.  0.  1.  0.  0.  0.  1.  0.  1.  1.  1.
  1.  1.  1.  0.  0.  0.  0.  1.  0.  0.  0.  0.  1.  0.  1.  1.]
Total number of TARGET subjects:  61
Total number of NON-TARGET subjects:  63
Stage 1 number of hits (true and false positives):  58.0
Stage 1 TRUE positives:  30.0
Stage 1 FALSE positives:  28.0
Stage 1 TRUE negatives:  35.0
Total number of flagged HPC-AD subjects:  8.0
Number of flagged HPC-AD subjects that are TRUE positives:  7
Number of flagged HPC-AD subjects that are FALSE positives:  1.0
Number of true negatives:  62.0
#######################[...]<jupyter_text>Do it without crossvalidation<jupyter_code>#reload(high_confidence)
hpc = high_confidence.TwoStagesPrediction(
    n_iter=500,
    shuffle_test_split=0.5,
    min_gamma=.9,
    thresh_ratio=0.1, gamma=0.9)

hpc.fit(x_, x_, y_)
hpc.gamma
array_results, dic_results = hpc.predict(x_, x_)
y_pred = (dic_results['s1df'][:,0]>0).astype(float)
lr_decision = dic_results['s2df'][:,1]
predic_stats(y_, y_pred, lr_decision)
sample.shape
sample['DX_GROUP'].value_counts()
array_results.shape
def predic_stats(y_, y_pred, lr_decision):
    # number of AD subjects
    n_ad = sum(y_)
    print('Total number of TARGET subjects: ', n_ad)

    # number of CN subjects
    n_cn = len(y_) - sum(y_)
    print('Total number of NON-TARGET subjects: ', n_cn)
    
    # number of subjects predicted as AD at stage 1
    n_pos = sum(y_pred)
    print('Stage 1 number of hits (true and false positives): ', n_pos)
    
    # true positives at stage 1
    n_pos_ad = sum(y_pred[y_.astype(bool)])
    print('Stage 1 TRUE positives: ', n_pos_ad)
    
    # false positives at stage 1
    n_pos_cn = n_pos - n_pos_ad
    print('Stage 1 FALSE positives: ', n_pos_cn)
    
    # number of CN subjects not identified as positive (true negatives)
    n_neg1_cn = n_cn - n_pos_cn
    print('Stage 1 TRUE negatives: ', n_neg1_cn)

    # number of all flagged HPC-AD subjects
    n_flag = sum(y_pred[lr_decision>0])
    print('Total number of flagged HPC-AD subjects: ', n_flag)

    # number of flagged HPC-AD subjects who are actually AD (true positives)
    n_flag_ad = sum(y_[lr_decision>0])
    print('Number of flagged HPC-AD subjects that are TRUE positives: ', n_flag_ad)

    # number of flagged HPC-AD subjects that are actually CN (false positives)
    n_flag_cn = n_flag - n_flag_ad
    print('Number of flagged HPC-AD subjects that are FALSE positives: ', n_flag_cn)

    # number of CN subjects that were not flagged (true negatives)
    n_neg_cn = n_cn - n_flag_cn
    print('Number of true negatives: ', n_neg_cn)
    
    print('#############################')
    print('Stage 1 stats for TARGET vs NON-TARGET')
    print('Precision for AD: ', n_pos_ad/(n_pos_ad + n_pos_cn))
    print('Recall (or sensitivity)  for AD: ', n_pos_ad/n_ad)
    sens = n_pos_ad/n_ad
    print('Specificity: ', n_neg1_cn/n_cn)
    spec = n_neg1_cn/n_cn
    fp = (1-spec)*664
    tp = sens*336
    adj_prec = tp/(tp+fp)
    print('Adjusted precision for 33.6% baseline rate: ', adj_prec)
    print('Accuracy: ', (n_pos_ad + n_neg1_cn)/(n_ad + n_cn))

    print('#############################')
    print('Stage 2 stats for TARGET vs NON-TARGET')
    print('Precision for HPC-AD: ', n_flag_ad/n_flag)
    print('Recall (or sensitivity) for HPC-AD: ', n_flag_ad/n_ad)
    sens_2 = n_flag_ad/n_ad
    print('Specificity: ', n_neg_cn/n_cn)
    spec_2 = n_neg_cn/n_cn
    fp_2 = (1-spec_2)*664
    tp_2 = sens_2*336
    adj_prec_2 = tp_2/(tp_2 + fp_2)
    print('Adjusted precision for 33.6% baseline rate: ', adj_prec_2)
    print('Accuracy: ', (n_flag_ad + n_neg_cn)/(n_ad + n_cn))
plt.hist(hpc.training_hit_probability,20)
plt.ylim(0,300)
plt.rcParams.update({'font.size': 18})
# Make some figures

predic_stats()
dic_results.keys()
def stats_mask(y_true, y_pred, mask_selected=None):
    if mask_selected is None:
        mask_selected = np.ones(y_pred.shape).astype(bool)
    print('------------------------')
    print('Ratio:', y_true[mask_selected].sum()/y_true.sum()) 
    print('#    : ', y_true[mask_selected].sum()) 
    print('# true values: ',mask_selected.sum())
    print('ACC  : ', np.mean((y_true == y_pred)[mask_selected]))
print('Level 1')
stats_mask(y_, (dic_results['s1df'][:,0]>0).astype(float))
pos_mask = (dic_results['s2df'][:,1]>0)
acc_s2 = metrics.accuracy_score(y_val[tmp_mask], (dic_results['s1df'][:,0]>0).astype(float)[pos_mask])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14760746147607462, 'fraction_numerical': 0.030638911417500227, 'mean_word_length': 3.3761829652996846, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '491115', 'license': 'MIT', 'n_tokens_mistral': 4545, 'n_tokens_neox': 4347, 'n_words': 984, 'path': 'Analysis/hps_v1.ipynb', 'repo_name': 'surchs/ABIDE_HPS', 'size': '11097'}