<jupyter_start><jupyter_text>Photometry Classification Pipeline For PLASTICC Dataset<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import math
import seaborn as sns

import sys

from sklearn.model_selection import cross_val_score
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import roc_curve, roc_auc_score, accuracy_score, make_scorer, f1_score, recall_score, precision_score
from sklearn.model_selection import cross_val_score, StratifiedKFold, cross_val_predict, train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV
from sklearn.calibration import CalibratedClassifierCV
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import MinMaxScaler

import copy
from sklearn.model_selection import StratifiedKFold

from sklearn.cluster import KMeans, DBSCAN
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA

from sklearn.calibration import CalibratedClassifierCV
from sklearn.model_selection import train_test_split

from IPython.display import clear_output
from sklearn.metrics import silhouette_score, silhouette_samples
from matplotlib import cm

from utils import binarization, interpolation, rolling_diff, rolling_ratio, log10_, KFoldsClassifier

from scipy import stats

from IPython.display import Markdown, display

def printmd(string):
    display(Markdown(string))
    
import copy
import datetime
    
pd.options.display.max_columns = 500<jupyter_output><empty_output><jupyter_text>Data Download<jupyter_code>!python download_data.py<jupyter_output>Checksum test for ../data/plasticc/note1_dataRelease.pdf PASSED.
Checksum test for ../data/plasticc/note2_modelNames.pdf PASSED.
Checksum test for ../data/plasticc/plasticc_modelpar.tar.gz PASSED.
Checksum test for ../data/plasticc/plasticc_test_lightcurves_01.csv.gz PASSED.
Checksum test for ../data/plasticc/plasticc_test_lightcurves_06.csv.gz PASSED.
Checksum test for ../data/plasticc/plasticc_test_lightcurves_02.csv.gz PASSED.
Checksum test for ../data/plasticc/plasticc_test_lightcurves_04.csv.gz PASSED.
Checksum test for ../data/plasticc/plasticc_test_lightcurves_08.csv.gz PASSED.
Checksum test for ../data/plasticc/plasticc_test_lightcurves_07.csv.gz PASSED.
Checksum test for ../data/plasticc/plasticc_test_lightcurves_03.csv.gz PASSED.
Checksum test for ../data/plasticc/plasticc_test_lightcurves_05.csv.gz PASSED.
Checksum test for ../data/plasticc/plasticc_train_metadata.csv.gz PASSED.
Checksum test for ../data/plasticc/plasticc_test_metadata.csv.gz PASSED.
Checksum test for ../dat[...]<jupyter_text>Data Import<jupyter_code>metadata = pd.read_csv('../data/plasticc/plasticc_test_metadata.csv.gz')

def set_sn_type(x):
#     if x in {90, 67, 52}:
    if x == 90:
        return 1
    else:
        return 0
    
metadata['is_ia'] = metadata.true_target.apply(set_sn_type)
def agg_to_arr(arr):
    return arr

def get_rms(row):
    weights = 1 / row.flux_err
    mean = np.average(row.flux, weights=weights)
    rms_for_pval = np.sum(np.square((row.flux - mean) * weights))
    return rms_for_pval

def get_pval(row):
    pval = 1.0 - stats.chi2.cdf(row.rms, row.size_r - 1)
    return pval

def normalize(x):
    return np.array(MinMaxScaler().fit_transform(np.array(x).reshape(-1, 1))).reshape(1, -1)[0]

def shifting(row):
    x = row.mjd - row.mjd[row.y_normed.argmax()]
    return x

def get_bins(row):
    return binarization(row.x_shifted, row.y_normed, n=16, x_ranges=[-50, 100])

def interpolate(row):
    return interpolation(np.arange(len(row.binarized)), row.binarized)

def generate_features(row):
    return rolling_ratio(row.interpolated)
X = 0
y = 0
flg = 0
for i in range(12):
    if i == 11:
        print('Starting File with Train')
        filename = '../data/plasticc/plasticc_train_lightcurves.csv.gz'
        metadata = pd.read_csv('../data/plasticc/plasticc_train_metadata.csv.gz')
        metadata['is_ia'] = metadata.true_target.apply(set_sn_type)
    else:
        print('Starting File #{0}'.format(i + 1))
        filename = '../data/plasticc/plasticc_test_lightcurves_{:02d}.csv.gz'.format(i + 1)
    data = pd.read_csv(filename)
    print('Read File OK ' + str(datetime.datetime.now().time()))
    
    data = data[(data.passband == 2) & (data.detected_bool == 1)]
    data = data.groupby(by='object_id').agg({'mjd' : list, 
                                         'flux' : list, 
                                         'flux_err' : list, 
                                         'detected_bool' : 'count'})
    print('Group Points By Object OK ' + str(datetime.datetime.now().time()))
    data.rename({'detected_bool' : 'size_r'}, axis=1, inplace=True)
    data.mjd = data.mjd.apply(lambda x: np.array(x))
    data.flux = data.flux.apply(lambda x: np.array(x))
    data.flux_err = data.flux_err.apply(lambda x: np.array(x))
    data = data.join(metadata.set_index('object_id')[['true_target', 'is_ia', 'true_z']])
    print('Join Data with Metadata OK ' + str(datetime.datetime.now().time()))
    data = data[data.true_target.isin({90, 67, 52, 42, 62, 95, 15, 88})]
    data.drop(['true_target'], axis=1, inplace=True)
    data['rms'] = data.apply(get_rms, axis=1)
    data['pval'] = data.apply(get_pval, axis=1)
    print('Calculate P-value OK ' + str(datetime.datetime.now().time()))
    data.pval.fillna(0, inplace=True)
    data = data[(data.pval < 0.001) & (data.size_r >= 4)].reset_index()
    data['y_normed'] = data.flux.apply(normalize)
    data['x_shifted'] = data.apply(shifting, axis=1)
    data['binarized'] = data.apply(get_bins, axis=1)
    data['interpolated'] = data.apply(interpolate, axis=1)
    print('Preprocess OK ' + str(datetime.datetime.now().time()))
    data['rolling_ratios'] = data.apply(generate_features, axis=1)
    print('Generate Features OK ' + str(datetime.datetime.now().time()))
    features_df = pd.DataFrame(np.concatenate((np.stack(data.interpolated.values),
                                               np.stack(data.rolling_ratios.values)), axis=1))
    features_df.columns = ['x_' + str(i) for i in range(16)] + ['d_' + str(i) for i in range(16)]
    xdata = pd.concat([data.reset_index(drop=True), features_df.reset_index(drop=True)], axis=1)
    print('Concatenate Features to Initial Dataset OK ' + str(datetime.datetime.now().time()))
    _X = xdata.loc[:, 'x_0':'d_15'].values
    _y = xdata['is_ia'].values
    if flg == 1:
        X = np.concatenate((X, _X))
        y = np.concatenate((y, _y))
    elif flg == 0:
        X = copy.deepcopy(_X)
        y = copy.deepcopy(_y)
        flg = 1
    del _X
    del _y
    print('Append to X, y OK ' + str(datetime.datetime.now().time()))
print(np.histogram(data['true_z'], bins=10, range=(0,2)))
print(np.count_nonzero(data['true_z'] <= 0.6) / data.shape[0])
print(np.histogram(data['size_r'], bins=10, range=(0,100)))
# import pickle

# with open('y_plasticc.pickle', 'rb') as f:
#     y = pickle.load(f)
    
# with open('X_plasticc.pickle', 'rb') as f:
#     X = pickle.load(f)
y.sum() / len(y)
y.sum()
X.shape
y.shape
a = np.concatenate((X, y.reshape(-1, 1)), axis=1)
not_ia = a[a[:, -1] == 0]
ia = a[a[:, -1] == 1]
shape = not_ia[:, -1].shape
n_samples = 59047
inds = np.random.randint(0, shape[0], size=n_samples)
notia_sample = not_ia[inds]
notia_sample.shape
b = np.concatenate((ia, notia_sample), axis=0)
b.shape
X = b[:, :-1]
y = b[:, -1]
X.shape
y.shape<jupyter_output><empty_output><jupyter_text>Train Test Split<jupyter_code>X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)<jupyter_output><empty_output><jupyter_text>Models  Random Forest  Train: PLASTiCC  Test: PLASTiCC<jupyter_code>rf = RandomForestClassifier(n_estimators=5000,
                            min_samples_split=10, 
                            min_samples_leaf=2, 
                            max_features='auto',
                            bootstrap=True, 
                            max_depth=30, 
                            n_jobs=-1, 
                            random_state=25, verbose=1)

rf.fit(X_train, y_train)
rf_pred = rf.predict_proba(X_test)<jupyter_output>[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done  42 tasks      | elapsed:   16.0s
[Parallel(n_jobs=-1)]: Done 192 tasks      | elapsed:  1.2min
[Parallel(n_jobs=-1)]: Done 442 tasks      | elapsed:  2.7min
[Parallel(n_jobs=-1)]: Done 792 tasks      | elapsed:  4.9min
[Parallel(n_jobs=-1)]: Done 1242 tasks      | elapsed:  7.3min
[Parallel(n_jobs=-1)]: Done 1792 tasks      | elapsed: 10.9min
[Parallel(n_jobs=-1)]: Done 2442 tasks      | elapsed: 14.8min
[Parallel(n_jobs=-1)]: Done 3192 tasks      | elapsed: 19.6min
[Parallel(n_jobs=-1)]: Done 4042 tasks      | elapsed: 24.9min
[Parallel(n_jobs=-1)]: Done 4992 tasks      | elapsed: 31.5min
[Parallel(n_jobs=-1)]: Done 5000 out of 5000 | elapsed: 31.6min finished
[Parallel(n_jobs=4)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done  42 tasks      | elapsed:    0.4s
[Parallel(n_jobs=4)]: Done 192 tasks      | elapsed:    1.8s
[Parallel(n_jobs=4[...]<jupyter_text>Random Forest  Train: PLASTiCC  Test: OSC<jupyter_code>import pickle
with open('X_osc.pickle', 'rb') as f:
    X_osc = pickle.load(f)
    
with open('y_osc.pickle', 'rb') as f:
    y_osc = pickle.load(f)
plasticc_osc = rf.predict_proba(X_osc)<jupyter_output>[Parallel(n_jobs=4)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done  42 tasks      | elapsed:    0.1s
[Parallel(n_jobs=4)]: Done 192 tasks      | elapsed:    0.4s
[Parallel(n_jobs=4)]: Done 442 tasks      | elapsed:    0.8s
[Parallel(n_jobs=4)]: Done 792 tasks      | elapsed:    1.5s
[Parallel(n_jobs=4)]: Done 1242 tasks      | elapsed:    2.3s
[Parallel(n_jobs=4)]: Done 1792 tasks      | elapsed:    3.3s
[Parallel(n_jobs=4)]: Done 2442 tasks      | elapsed:    4.4s
[Parallel(n_jobs=4)]: Done 3192 tasks      | elapsed:    5.8s
[Parallel(n_jobs=4)]: Done 4042 tasks      | elapsed:    7.2s
[Parallel(n_jobs=4)]: Done 4992 tasks      | elapsed:    8.5s
[Parallel(n_jobs=4)]: Done 5000 out of 5000 | elapsed:    8.5s finished<jupyter_text>Random Forest  Train: OSC Test: PLASTiCC<jupyter_code>with open('rf_osc_model.sav', 'rb') as f:
    osc_rf_model = pickle.load(f)

osc_plasticc = osc_rf_model.predict_proba(X)<jupyter_output><empty_output><jupyter_text>Random Forest  Train: OSC Test: OSC<jupyter_code>with open('OSC_rf_y_pred.pickle', 'rb') as f:
    osc_osc = pickle.load(f)
    
with open('OSC_y_true.pickle', 'rb') as f:
    osc_y_true = pickle.load(f)<jupyter_output><empty_output><jupyter_text>Model Evaluation<jupyter_code>scores_table = []

fig, ax = plt.subplots(1, 1, figsize=(5, 4), dpi=300)
ax.set_xlabel('False positive rate (FPR)')
ax.set_ylabel('True positive rate (TPR)')
ax.set_xlim(0., 1)
ax.set_ylim(0., 1)
ax.set_title(r'ROC Curves', fontsize=10)
ax.grid(b=0.1)

colors = ['C0', 'C1', 'C2', 'C3']
labels = ['Train: $Plasticc$ | Test: $Plasticc$',
         'Train: $Plasticc$ | Test: $OSC$',
         'Train: $OSC$ | Test: $Plasticc$',
         'Train: $OSC$ | Test: $OSC$']

for i, train_test in enumerate([(y_test, rf_pred),
                              (y_osc, plasticc_osc),
                              (y, osc_plasticc),
                              (osc_y_true, osc_osc)]):
    try:
        fpr, tpr, thresholds = roc_curve(train_test[0], train_test[1][:, 1])
        ax.plot(1-tpr, 1-fpr, linewidth=1.5, label=labels[i], color=colors[i])
        auc = roc_auc_score(train_test[0], train_test[1][:, 1])
        acc = accuracy_score(train_test[0], train_test[1][:, 1]>=0.5)
        f1 = f1_score(train_test[0], train_test[1][:, 1]>=0.5)
        pr = precision_score(train_test[0], train_test[1][:, 1]>=0.5)
        rec = recall_score(train_test[0], train_test[1][:, 1]>=0.5)
    except:
        fpr, tpr, thresholds = roc_curve(train_test[0], train_test[1])
        ax.plot(1-tpr, 1-fpr, linewidth=1.5, label=labels[i], color=colors[i])
        auc = roc_auc_score(train_test[0], train_test[1])
        acc = accuracy_score(train_test[0], train_test[1]>=0.5)
        f1 = f1_score(train_test[0], train_test[1]>=0.5)
        pr = precision_score(train_test[0], train_test[1]>=0.5)
        rec = recall_score(train_test[0], train_test[1]>=0.5)
    scores_table.append([auc, acc, f1, pr, rec])
    
ax.legend()
plt.tight_layout()
plt.savefig('images/roc_plasticc.pdf')
plt.show()
std_errors = {}
for i, train_test in enumerate([(y_test, rf_pred),
                              (y_osc, plasticc_osc),
                              (y, osc_plasticc),
                              (osc_y_true, osc_osc)]):
    bs_auc = []
    bs_acc = []
    bs_f1 = []
    bs_pr = []
    bs_rec = []
    
    for j in range(1000):
        idx = np.random.choice([i for i in range(len(train_test[0]))], len(train_test[0]))
        try:
            bs_auc.append(roc_auc_score(train_test[0][idx], train_test[1][idx][:, 1]))
            bs_acc.append(accuracy_score(train_test[0][idx], train_test[1][idx][:, 1]>=0.5))
            bs_f1.append(f1_score(train_test[0][idx], train_test[1][idx][:, 1]>=0.5))
            bs_pr.append(precision_score(train_test[0][idx], train_test[1][idx][:, 1]>=0.5))
            bs_rec.append(recall_score(train_test[0][idx], train_test[1][idx][:, 1]>=0.5))
        except:
            bs_auc.append(roc_auc_score(train_test[0][idx], train_test[1][idx]))
            bs_acc.append(accuracy_score(train_test[0][idx], train_test[1][idx]>=0.5))
            bs_f1.append(f1_score(train_test[0][idx], train_test[1][idx]>=0.5))
            bs_pr.append(precision_score(train_test[0][idx], train_test[1][idx]>=0.5))
            bs_rec.append(recall_score(train_test[0][idx], train_test[1][idx]>=0.5))
    std_errors[labels[i]] = [np.std(bs_auc), np.std(bs_acc), np.std(bs_f1), np.std(bs_pr), np.std(bs_rec)]
    print(std_errors)
    print('model done')
errs = pd.DataFrame(std_errors)
errs.columns = ['Plasticc-Plasticc', 'Plasticc-OSC', 'OSC-Plasticc', 'OSC-OSC']
errs.index = ['AUC', 'Acc', 'F1', 'Pr', 'Rec']
errs
pd.options.display.float_format = "{:.3f}".format

scores_df = pd.DataFrame(scores_table).T
scores_df.index = ['AUC', 'Acc', 'F1-Score', 'Precision', 'Recall']
scores_df.columns = ['Plasticc-Plasticc', 'Plasticc-OSC', 'OSC-Plasticc', 'OSC-OSC']
scores_df
np.sum(y) / len(y)
np.sum(y_osc) / len(y_osc)
from sklearn.metrics import confusion_matrix

cf_matrix = confusion_matrix(y_test, rf_pred[:, 1]>=0.5)

fig, ax = plt.subplots(figsize=(5, 4), dpi=200)

group_names = ['True Neg','False Pos','False Neg','True Pos']
group_counts = ["{0:0.0f}".format(value) for value in
                cf_matrix.flatten()]
group_percentages = ["{0:.2%}".format(value) for value in
                     cf_matrix.flatten()/np.sum(cf_matrix)]
labels = [f"{v1}\n{v2}\n{v3}" for v1, v2, v3 in
          zip(group_names,group_counts,group_percentages)]
labels = np.asarray(labels).reshape(2,2)
sns.heatmap(cf_matrix, annot=labels, fmt="", cmap='Blues', ax=ax)
plt.title('Plasticc-Plasticc',fontsize=16)
cf_matrix = confusion_matrix(y_osc, plasticc_osc[:, 1]>=0.5)

fig, ax = plt.subplots(figsize=(5, 4), dpi=200)

group_names = ['True Neg','False Pos','False Neg','True Pos']
group_counts = ["{0:0.0f}".format(value) for value in
                cf_matrix.flatten()]
group_percentages = ["{0:.2%}".format(value) for value in
                     cf_matrix.flatten()/np.sum(cf_matrix)]
labels = [f"{v1}\n{v2}\n{v3}" for v1, v2, v3 in
          zip(group_names,group_counts,group_percentages)]
labels = np.asarray(labels).reshape(2,2)
sns.heatmap(cf_matrix, annot=labels, fmt="", cmap='Blues', ax=ax)
plt.title('Plasticc-OSC',fontsize=16)
cf_matrix = confusion_matrix(y, osc_plasticc[:, 1]>=0.5)

fig, ax = plt.subplots(figsize=(5, 4), dpi=200)

group_names = ['True Neg','False Pos','False Neg','True Pos']
group_counts = ["{0:0.0f}".format(value) for value in
                cf_matrix.flatten()]
group_percentages = ["{0:.2%}".format(value) for value in
                     cf_matrix.flatten()/np.sum(cf_matrix)]
labels = [f"{v1}\n{v2}\n{v3}" for v1, v2, v3 in
          zip(group_names,group_counts,group_percentages)]
labels = np.asarray(labels).reshape(2,2)
sns.heatmap(cf_matrix, annot=labels, fmt="", cmap='Blues', ax=ax)
plt.title('OSC-Plasticc',fontsize=16)
cf_matrix = confusion_matrix(osc_y_true, osc_osc>=0.5)

fig, ax = plt.subplots(figsize=(5, 4), dpi=200)

group_names = ['True Neg','False Pos','False Neg','True Pos']
group_counts = ["{0:0.0f}".format(value) for value in
                cf_matrix.flatten()]
group_percentages = ["{0:.2%}".format(value) for value in
                     cf_matrix.flatten()/np.sum(cf_matrix)]
labels = [f"{v1}\n{v2}\n{v3}" for v1, v2, v3 in
          zip(group_names,group_counts,group_percentages)]
labels = np.asarray(labels).reshape(2,2)
sns.heatmap(cf_matrix, annot=labels, fmt="", cmap='Blues', ax=ax)
plt.title('OSC-OSC',fontsize=16)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15549065420560748, 'fraction_numerical': 0.032535046728971964, 'mean_word_length': 3.7426592797783935, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '458642', 'license': 'MIT', 'n_tokens_mistral': 6870, 'n_tokens_neox': 6449, 'n_words': 1205, 'path': 'notebooks/plasticc_pipeline.ipynb', 'repo_name': 'HSE-LAMBDA/supernovae_classification', 'size': '17120'}