<jupyter_start><jupyter_text>Testing KERAS classification with MLP<jupyter_code>from fastespy.mlkeras.evaluation import SplitData, plot_metric, plot_sig_vs_thr, get_sig_bkg_rate_eff, get_tp_fp_fn, learning_curve, SignificanceMetric
from fastespy.mlkeras import models
from fastespy.ml import MLHyperParTuning
from fastespy.readpydata import convert_data_to_ML_format
from fastespy import feldman_cousins as fc
from tensorflow.keras.layers.experimental.preprocessing import Normalization
from tensorflow import keras
import tensorflow as tf
import matplotlib.pyplot as plt
import sys
import glob
import os
import tempfile
import logging
import numpy as np
import tqdm
import time
from pathlib import PosixPath
from fastespy.analysis import init_logging
sys.path.append("/Users/manuelmeyer/Python/fastespy/fastespy/scripts/")
sys.path.append("/Users/manuelmeyer/Python/TESdata/fastespy/fastespy/scripts/")
from ml_intrinsic_bkg import load_data
%matplotlib inline
init_logging("INFO", color=True)<jupyter_output><empty_output><jupyter_text>Read data<jupyter_code>path = "/Users/manuelmeyer/Documents/TESData/IntrinsicsData_NewSetup/" # May 2021 data
prefix = "IntrinsicsData_NewSetup"

path = "/Users/manuelmeyer/Documents/TESData/August2020/" # May 2021 data
prefix = "August2020"

#path = "/Users/manuelmeyer/Python/TESdata/data-01152021//" # August 2020 data
#prefix = "data_august2020"

in_dir = PosixPath(path)
files = glob.glob(os.path.join(path, '*.npy'))

print(len(files))
#files

if not len(files):
    path = "/Users/manuelmeyer/Python/TESdata/data-05122021/"
    in_dir = PosixPath(path)
    files = glob.glob(os.path.join(path, '*.npy'))
    print(len(files))
files = list(in_dir.glob("*.npy"))

logging.info("Using files:")
for f in files:
    logging.info(f)

# define the feature names
feature_names = []
remove = ['data', 'time', 'pulse integral raw', 'voltage error',
          'error', 'start time in hrs', 'end time in hrs',
          'trigger time'
        ]

remove.append('constant')  # this seems to help

x = np.load(files[0], allow_pickle=True).tolist()
for k in x[1].keys():
    if not k in remove and not 'error' in k:
        feature_names.append(k)

logging.info("Using features names {}".format(feature_names))

 # define light cleaning cuts
light_cleaning_cuts = {
     "chi2 reduced": "chi2 < 6.",
    "decay time": "decay < 10.e-6",
    "trigger time": "(trigger >= 29.5e-6) & (trigger <= 30.8e-6)",  # from gaussian fit, 5 sigma interval
}

# read the data
result, data, t_tot_hrs = load_data(files, feature_names, light_cleaning_cuts=light_cleaning_cuts)

# convert data to ML format
X, y = convert_data_to_ML_format(result,
                                 feature_names,
                                 bkg_type=0,
                                 signal_type=1)<jupyter_output>[0;36m3247362142.py:[0;35m   3[0;0m --- [1;36mINFO[1;0m: Using files:
[0;36m3247362142.py:[0;35m   5[0;0m --- [1;36mINFO[1;0m: /Users/manuelmeyer/Documents/TESData/August2020/0.3RN-1.5GHzGBWP-intrinsics-50MHz-20mV-1day-16-fit000.npy
[0;36m3247362142.py:[0;35m   5[0;0m --- [1;36mINFO[1;0m: /Users/manuelmeyer/Documents/TESData/August2020/0.3RN-1.5GHzGBWP-intrinsics-50MHz-20mV-72hrs-3-fit000.npy
[0;36m3247362142.py:[0;35m   5[0;0m --- [1;36mINFO[1;0m: /Users/manuelmeyer/Documents/TESData/August2020/0.3RN-1.5GHz-light-000-fit.npy
[0;36m3247362142.py:[0;35m   5[0;0m --- [1;36mINFO[1;0m: /Users/manuelmeyer/Documents/TESData/August2020/0.3RN-1.5GBWP-intrinsics-50MHz-20mV-1day-9-fit000.npy
[0;36m3247362142.py:[0;35m   5[0;0m --- [1;36mINFO[1;0m: /Users/manuelmeyer/Documents/TESData/August2020/0.3RN-1.5GHzGBWP-intrinsics-50MHz-20mV-1day-13-fit000.npy
[0;36m3247362142.py:[0;35m   5[0;0m --- [1;36mINFO[1;0m: /Users/manuelmeyer/Documents/TESData/August2020/0.3RN-[...]<jupyter_text>Do a log transform<jupyter_code>X_log, y_log = MLHyperParTuning.transform_data_log(X.copy(), y.copy(), feature_names)<jupyter_output>Discarded 234 triggers in log transform<jupyter_text>Get a split of the data (will be repeated later on):<jupyter_code>sd = SplitData(X_log, y_log, n_splits=5, stratify=True, random_state=42)
i_split = 0
#X_train, X_test, y_train, y_test = sd.get_split(i_split)
X_train, X_val, y_train, y_val = sd.get_split_with_test_set(m=i_split, n=0)<jupyter_output><empty_output><jupyter_text>Normalize the data<jupyter_code># normalize the input data
normalizer = Normalization(axis=-1)

# adapt to training data
normalizer.adapt(X_train)

# normalize training data
X_train = normalizer(X_train)
X_val = normalizer(X_val)
X_test = normalizer(sd.X_test)<jupyter_output><empty_output><jupyter_text>Define the model<jupyter_code>initial_bias, initial_loss = models.initial_output_bias(y_train)
print(initial_bias, initial_loss)
n_nodes = 100
n_layers = 5
l2_regularizer = None
dropout = 0.1
sm = SignificanceMetric(N_tot=y_train.size + y_val.size + sd.y_test.size,
                        t_obs=t_tot_hrs * 3600.,
                        #dtype=tf.float
                       )
metrics = list(models.metrics)
metrics.append(sm)
model = models.make_mlp_model(n_features=X_train.shape[1],
                              n_layers=n_layers,
                              n_nodes=n_nodes,
                              l2_regularizer=l2_regularizer,
                              dropout=dropout,
                              output_bias=initial_bias,
                              learning_rate=3e-4,
                              metrics=metrics
                             )<jupyter_output><empty_output><jupyter_text>Initial predictions<jupyter_code>model.predict(X_train[10:]), y_train.sum() / y_train.size
results = model.evaluate(X_train, y_train, verbose=0, return_dict=True)
print(results)
results['loss'], initial_loss<jupyter_output><empty_output><jupyter_text>Save the inital weights<jupyter_code>initial_weights = os.path.join(tempfile.mkdtemp(), f'initial_weights_{n_layers}_{n_nodes}')
model.save_weights(initial_weights)<jupyter_output><empty_output><jupyter_text>Train the full model<jupyter_code>model.summary()
model.load_weights(initial_weights)
keras.utils.plot_model(model, "model.png")
EPOCHS = 200
BATCH_SIZE = 2048  # large enough so that you have enough signal samples in each batch

baseline_history = models.train_model(model, X_train, y_train, epochs=EPOCHS, batch_size=BATCH_SIZE, X_val=X_val, y_val=y_val)
fig=plt.figure(figsize=(8,10))
ax = plt.subplot(311)
plot_metric(baseline_history, ax=ax, lw=2)
ax.tick_params(direction='in', labelbottom=False)
ax.set_xlabel("")
ax.grid(which='both')
plt.legend()

ax = plt.subplot(312)
plot_metric(baseline_history, ax=ax, metric='accuracy', lw=2)
ax.tick_params(direction='in', labelbottom=False)
ax.set_xlabel("")
ax.grid(which='both')
plt.legend()

ax = plt.subplot(313)
plot_metric(baseline_history, ax=ax, metric='significance', lw=2)
ax.grid(which='both')
plt.legend()

plt.subplots_adjust(hspace=0.)

plt.savefig(f"{prefix:s}_metrics_vs_epochs.png")
results = model.evaluate(X_train, y_train, batch_size=BATCH_SIZE, verbose=0, return_dict=True)
print(results)<jupyter_output>{'loss': 0.0027402015402913094, 'tp': 680.0, 'fp': 16.0, 'tn': 25164.0, 'fn': 3.0, 'accuracy': 0.9992653727531433, 'precision': 0.977011501789093, 'recall': 0.9956076145172119, 'auc': 0.99991774559021, 'prc': 0.9957250952720642, 'significance': 4.281110763549805}<jupyter_text>Evaluate metrics on training and validation set<jupyter_code>thr = 0.5
y_pred = model.predict(X_train)

sig, bkg_rate, eff = get_sig_bkg_rate_eff(y_train,
                                          y_pred,
                                          N_tot=y_train.size + y_val.size + sd.y_test.size,
                                          t_obs=t_tot_hrs * 3600.,
                                          thr=thr
                                         )
print(sig, bkg_rate, eff)
thr = 0.5
y_pred = model.predict(X_val)
#class_pred_test, tp, fp, fn = get_tp_fp_fn(y_val, y_pred, thr=thr)


sig, bkg_rate, eff = get_sig_bkg_rate_eff(y_val,
                                          y_pred,
                                          N_tot=y_train.size + y_val.size + sd.y_test.size,
                                          t_obs=t_tot_hrs * 3600.,
                                          thr=thr
                                         )
print(sig, bkg_rate, eff)<jupyter_output>3.702799436482806 2.0112795085573605e-05 1.0<jupyter_text>Retrain with class weights<jupyter_code># weighting following
# https://www.tensorflow.org/tutorials/structured_data/imbalanced_data#class_weights
# Scaling by total/2 helps keep the loss to a similar magnitude.
# The sum of the weights of all examples stays the same.

# original: increasing weight for class 1:
# false negatives go down, false positives go up
# significance goes down
weight_for_0 = 1 / (y.size - y.sum()) * (y.size / 2.)
weight_for_1 = 1 / y.sum() * (y.size / 2.)

# inverting weights: fp down, fn up, significance up!
power = 1.
weight_for_0 = 1. / (1 / (y.size - y.sum()) * (y.size / 2.))**power
weight_for_1 = 1. / (1 / y.sum() * (y.size / 2.))**power

print('Weight for class 0: {:.3f}'.format(weight_for_0))
print('Weight for class 1: {:.3f}'.format(weight_for_1))


use_class_weights = False

if use_class_weights:
    model.load_weights(initial_weights)
    print({0: weight_for_0, 1: weight_for_1})

    class_weight_history = model.fit(
        X_train,
        y_train,
        batch_size=BATCH_SIZE,
        epochs=EPOCHS,
        validation_data=(X_val, y_val), 
        verbose=0,
        class_weight={0: weight_for_0, 1: weight_for_1},
    )
if use_class_weights:
    result = model.evaluate(X_val, return_dict=True)
    print(result)
    
    plot_metric(class_weight_history)
    plt.legend()
    
    plot_metric(class_weight_history, metric='significance')
    plt.legend()
    
    thr = 0.5
    y_pred = model.predict(X_train)

    sig, bkg_rate, eff = get_sig_bkg_rate_eff(y_train,
                                          y_pred,
                                          N_tot=y_train.size + y_val.size + sd.y_test.size,
                                          t_obs=t_tot_hrs * 3600.,
                                          thr=thr
                                         )
    print(sig, bkg_rate, eff)
    
    y_pred = model.predict(X_val)

    sig, bkg_rate, eff = get_sig_bkg_rate_eff(y_val,
                                          y_pred,
                                          N_tot=y_train.size + y_val.size + sd.y_test.size,
                                          t_obs=t_tot_hrs * 3600.,
                                          thr=thr
                                         )
    print(sig, bkg_rate, eff)<jupyter_output><empty_output><jupyter_text>results as function of power in weight definition for May 2021 data| power | weight 0 | weight 1 | sig val | bkg rate (Hz) val | eff val | |---|---|---|---|---|---|| 0 | 1 | 1 | 1.95 | 8.20e-5 | 0.99 | | 1 | 1.53 | 0.47 | 2.26 | 5.7e-5 | 0.97 | | 2 | 2.33 | 0.22 | 2.34 | 4.99e-5 | 0.95 | | 3 | 3.558 | 0.106 | 2.91 | 2.49e-5 | 0.86 || 4 | 5.43 | 0.05 | 3.73 | 1.07e-5 | 0.80 | | 5 | 8.29 | 0.02 | 0. | 0. | 0. | As can be seen, class weights have similar effet as setting the threshold to a higher value than 0.5. FP goes down, TP also down. For August 2020 data:| power | weight 0 | weight 1 | sig val | bkg rate (Hz) val | eff val | |---|---|---|---|---|---|| 0 | 1 | 1 | 3.97 | 1.68e-05 | 1.0 | | 1 | 1.95 | 0.05 | 4.27 | 6.7e-6 | 0.75 | | 1.2 | 2.225 | 0.03 | 5.18 | 3.35e-6 | 0.75 | | 1.5 | 2.72 | 0.01 | 0 | 0| 0. | | 2 | 3.97 | 0.003 | 0. | 0. | 0. |   Significance as function of threshold For validation set<jupyter_code>fig=plt.figure(dpi=120, figsize=(6,2*3))
ax, threshold, significance, bkg_rate, eff = plot_sig_vs_thr(model,
                                                             X_val,
                                                             y_val,
                                                             t_obs_hours=t_tot_hrs,
                                                             N_tot=y_train.size + y_val.size + sd.y_test.size)
plt.subplots_adjust(hspace=0.)<jupyter_output>Max significance: 6.02 for threshold 0.9989 background rate 0.00e+00 and analysis efficiency 0.35<jupyter_text>For training set<jupyter_code>fig=plt.figure(dpi=120, figsize=(6,2*3))
ax, threshold, significance, bkg_rate, eff = plot_sig_vs_thr(model,
                                                             X_train,
                                                             y_train,
                                                             t_obs_hours=t_tot_hrs,
                                                             N_tot=y_train.size + y_val.size + sd.y_test.size)
plt.subplots_adjust(hspace=0.)
plt.savefig(f"{prefix:s}_perfomance_vs_thr.png")<jupyter_output>Max significance: 7.53 for threshold 0.9972 background rate 0.00e+00 and analysis efficiency 0.54<jupyter_text>Significance for test set with the threshold obtained from training set:<jupyter_code>imax = np.argmax(significance)
y_pred = model.predict(X_val)
thr = threshold[imax]
#thr = 0.5
class_pred_val, tp, fp, fn = get_tp_fp_fn(y_val, y_pred, thr=thr)
print(thr)
sig, bkg_rate, eff = get_sig_bkg_rate_eff(y_val,
                                          y_pred,
                                          N_tot=y_train.size + y_val.size + sd.y_test.size,
                                          t_obs=t_tot_hrs * 3600.,
                                          thr=thr
                                         )
print(sig, bkg_rate, eff)
y_pred = model.predict(normalizer(sd.X_test))
thr = threshold[imax]
thr = 0.5
class_pred_test, tp, fp, fn = get_tp_fp_fn(sd.y_test, y_pred, thr=thr)

sig, bkg_rate, eff = get_sig_bkg_rate_eff(sd.y_test,
                                          y_pred,
                                          N_tot=y_train.size + y_val.size + sd.y_test.size,
                                          t_obs=t_tot_hrs * 3600.,
                                          thr=thr
                                         )
print(sig, bkg_rate, eff)<jupyter_output>4.262981638138017 1.3407700629234085e-05 0.9906103286384976<jupyter_text>Plot false positives for test set<jupyter_code>d = {}
d['y_test'] = sd.y_test
scorer = 'keras_loss'
d['y_pred_test'] = {scorer: class_pred_test}
d['idx_test'] = sd.idx_test
d['prob_test'] = {scorer: np.hstack([1. - y_pred, y_pred])}
d['classifier'] = "keras_mlp"

if fp.sum() < 50:
    ax = MLHyperParTuning.plot_misidentified_time_lines(
        d, scorer,
        data['time'],
        data['data'],
        X=X,  # give original data, for right values in legend
        feature_names=feature_names,
        plot_false_positive=True, save_plot=False
    )
plt.savefig(f"{prefix:s}_{thr:.2f}_miss_ids.png")<jupyter_output><empty_output><jupyter_text>Learning curve<jupyter_code>np.random.seed(42)
idx = np.arange(y_log.size)
np.random.shuffle(idx)
print(X_log[idx].shape)

sample_sizes, results_train, results_test = learning_curve(model, X_log[idx], y_log[idx], random_state=42, normalizer=normalizer, epochs=100)
plt.figure(figsize=(8,6))
for i, k in enumerate(['accuracy', 'significance']):
    ax = plt.subplot(2,1,i+1)
    if not i:
        ax.tick_params(direction='in', labelbottom=False)
    plt.errorbar(sample_sizes, results_train[k].mean(axis=-1), yerr=np.sqrt(results_train[k].var(axis=-1)), marker='o', ls='-', color="C0", label="Train")
    plt.errorbar(sample_sizes, results_test[k].mean(axis=-1), yerr=np.sqrt(results_test[k].var(axis=-1)), marker='o', label="Test", ls='--', color="C1")
    plt.legend()
    ax.grid(which='both')
    ax.set_ylabel(k)
plt.subplots_adjust(hspace=0.)
ax.set_xlabel("Sample size")
plt.savefig(f"{prefix:s}_learning_curve.png")
model.evaluate(X_train, y_train, batch_size=BATCH_SIZE, verbose=0, return_dict=True)<jupyter_output><empty_output><jupyter_text>Repeat analysis over different folds<jupyter_code>n_nodes = 100
n_layers = 5
l2_regularizer = None
dropout = 0.1
sig, br, ef = np.zeros(5), np.zeros(5), np.zeros(5)

for i in range(5):
    print(i)
    X_train, X_val, y_train, y_val = sd.get_split_with_test_set(m=i, n=0)
    
    # normalize the input data
    normalizer = Normalization(axis=-1)

    # adapt to training data
    normalizer.adapt(X_train)

    # normalize training data
    X_train = normalizer(X_train)
    X_val = normalizer(X_val)
    
    # calculate initial bias
    initial_bias, initial_loss = models.initial_output_bias(y_train)
    
    # the model
    model = models.make_mlp_model(n_features=X_train.shape[1],
                              n_layers=n_layers,
                              n_nodes=n_nodes,
                              l2_regularizer=l2_regularizer,
                              dropout=dropout,
                              output_bias=initial_bias,
                              learning_rate=3e-4
                             )
    
    # fit the model
    models.train_model(model,
                       X_train,
                       y_train,
                       epochs=EPOCHS,
                       batch_size=BATCH_SIZE,
                       X_val=X_val,
                       y_val=y_val,
                       class_weight=None
                      )
    
    # compute best threshold 
    ax, threshold, significance, bkg_rate, eff = plot_sig_vs_thr(model,
                                                                 X_train,
                                                                 y_train,
                                                                 t_obs_hours=t_tot_hrs,
                                                                 N_tot=y_train.size + y_val.size + sd.y_test.size)
    
    imax = np.argmax(significance)
    y_pred = model.predict(X_val)
    class_pred_test, tp, fp, fn = get_tp_fp_fn(y_val, y_pred, thr=threshold[imax])

    sig[i], br[i], ef[i] = get_sig_bkg_rate_eff(y_val,
                                                y_pred,
                                                N_tot=y_train.size + y_val.size + sd.y_test.size,
                                                t_obs=t_tot_hrs * 3600.,
                                               # thr=threshold[imax]
                                                 thr=0.5
                                               )
    
    del model
mean = np.vstack([sig, br, ef]).mean(axis=1)
std = np.sqrt(np.vstack([sig, br, ef]).var(axis=1))
maximum = np.vstack([sig, br, ef]).max(axis=1)
minimum = np.sqrt(np.vstack([sig, br, ef]).min(axis=1))

print(mean)
print(std)
print(minimum)
print(maximum)<jupyter_output>[5.42258357e+00 1.34085301e-05 9.78933609e-01]
[2.55901644e+00 1.01675213e-05 2.23698488e-02]
[1.80836806 0.         0.96730167]
[1.01578297e+01 2.68170601e-05 1.00000000e+00]
{'fraction_non_alphanumeric': 0.12184494393963548, 'fraction_numerical': 0.055635262234975294, 'mean_word_length': 2.3333333333333335, 'pattern_counts': {'":': 3, '<': 65, '<?xml version=': 0, '>': 62, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 93, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '584319', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7223, 'n_tokens_neox': 6438, 'n_words': 1169, 'path': 'notebooks/keras_workflow.ipynb', 'repo_name': 'me-manu/fastespy', 'size': '18819'}