In [1]:
import sys
sys.path.append('/home/jbourbeau/cr-composition')
sys.path
Out[1]:
In [2]:
import argparse
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import seaborn.apionly as sns
from sklearn.metrics import accuracy_score
from sklearn.model_selection import validation_curve, GridSearchCV, cross_val_score, ParameterGrid
from sklearn.neighbors import KNeighborsClassifier
import composition as comp
%matplotlib inline
In [3]:
sns.set_palette('muted')
sns.set_color_codes()
In [4]:
df, cut_dict = comp.load_sim(return_cut_dict=True)
selection_mask = np.array([True] * len(df))
standard_cut_keys = ['reco_exists', 'reco_zenith', 'NStations', 'IT_signal',
'StationDensity', 'reco_containment', 'energy_range']
for key in standard_cut_keys:
selection_mask *= cut_dict[key]
df = df[selection_mask]
feature_list, feature_labels = comp.get_training_features()
print(feature_list)
X_train, X_test, y_train, y_test, le = comp.get_train_test_sets(
df, feature_list)
print('number training events = ' + str(y_train.shape[0]))
Get error in charge vs. energy distribution
In [ ]:
In [5]:
pipeline = comp.get_pipeline('RF')
param_range = np.arange(1, 20)
param_grid = {'classifier__max_depth': param_range}
gs = GridSearchCV(estimator=pipeline,
param_grid=param_grid,
scoring='accuracy',
cv=5,
verbose=1,
n_jobs=10)
gs = gs.fit(X_train, y_train)
print('best GS CV score = {}'.format(gs.best_score_))
print('best GS CV depths = {}'.format(gs.best_params_))
print('Grid scores on development set:')
means = gs.cv_results_['mean_test_score']
stds = gs.cv_results_['std_test_score']
for mean, std, params in zip(means, stds, gs.cv_results_['params']):
print("%0.3f (+/-%0.03f) for %r"
% (mean, std * 2, params))
pipeline.set_params(**gs.best_params_)
pipeline.fit(X_train, y_train)
scaler = pipeline.named_steps['scaler']
clf = pipeline.named_steps['classifier']
In [6]:
clf_name = clf.__class__.__name__
print('=' * 30)
print(clf_name)
test_predictions = pipeline.predict(X_test)
test_acc = accuracy_score(y_test, test_predictions)
print('Test accuracy: {:.4%}'.format(test_acc))
train_predictions = pipeline.predict(X_train)
train_acc = accuracy_score(y_train, train_predictions)
print('Train accuracy: {:.4%}'.format(train_acc))
scores = cross_val_score(
estimator=pipeline, X=X_test, y=y_test, cv=10, n_jobs=10)
print('CV score: {:.2%} (+/- {:.2%})'.format(scores.mean(), scores.std()))
print('=' * 30)
In [7]:
correctly_identified_mask = (test_predictions == y_test)
energy_bin_width = 0.1
energy_bins = np.arange(6.2, 8.1, energy_bin_width)
# energy_bins = np.arange(6.2, 9.51, energy_bin_width)
energy_midpoints = (energy_bins[1:] + energy_bins[:-1]) / 2
log_energy = X_test[:, 0]
MC_proton_mask = (le.inverse_transform(y_test) == 'P')
MC_iron_mask = (le.inverse_transform(y_test) == 'Fe')
# Get number of MC proton and iron as a function of MC energy
num_MC_protons_energy = np.histogram(log_energy[MC_proton_mask],
bins=energy_bins)[0]
num_MC_protons_energy_err = np.sqrt(num_MC_protons_energy)
num_MC_irons_energy = np.histogram(log_energy[MC_iron_mask],
bins=energy_bins)[0]
num_MC_irons_energy_err = np.sqrt(num_MC_irons_energy)
num_MC_total_energy = np.histogram(log_energy, bins=energy_bins)[0]
num_MC_total_energy_err = np.sqrt(num_MC_total_energy)
# Get number of reco proton and iron as a function of MC energy
num_reco_proton_energy = np.histogram(
log_energy[MC_proton_mask & correctly_identified_mask],
bins=energy_bins)[0]
num_reco_proton_energy_err = np.sqrt(num_reco_proton_energy)
num_reco_iron_energy = np.histogram(
log_energy[MC_iron_mask & correctly_identified_mask],
bins=energy_bins)[0]
num_reco_iron_energy_err = np.sqrt(num_reco_iron_energy)
num_reco_total_energy = np.histogram(
log_energy[correctly_identified_mask],
bins=energy_bins)[0]
num_reco_total_energy_err = np.sqrt(num_reco_total_energy)
# Calculate reco proton and iron fractions as a function of MC energy
reco_proton_frac, reco_proton_frac_err = comp.ratio_error(
num_reco_proton_energy, num_reco_proton_energy_err,
num_MC_protons_energy, num_MC_protons_energy_err)
reco_iron_frac, reco_iron_frac_err = comp.ratio_error(
num_reco_iron_energy, num_reco_iron_energy_err,
num_MC_irons_energy, num_MC_irons_energy_err)
reco_total_frac, reco_total_frac_err = comp.ratio_error(
num_reco_total_energy, num_reco_total_energy_err,
num_MC_total_energy, num_MC_total_energy_err)
In [8]:
# Plot fraction of events vs energy
fig, ax = plt.subplots()
ax.errorbar(energy_midpoints, reco_proton_frac,
yerr=reco_proton_frac_err,
# xerr=energy_bin_width / 2,
marker='.', markersize=10,
label='Proton')
ax.errorbar(energy_midpoints, reco_iron_frac,
yerr=reco_iron_frac_err,
# xerr=energy_bin_width / 2,
marker='.', markersize=10,
label='Iron')
ax.errorbar(energy_midpoints, reco_total_frac,
yerr=reco_total_frac_err,
# xerr=energy_bin_width / 2,
marker='.', markersize=10,
label='Total')
plt.xlabel('$\log_{10}(E_{\mathrm{reco}}/\mathrm{GeV})$')
ax.set_ylabel('Fraction correctly identified')
ax.set_ylim([0.0, 1.0])
ax.set_xlim([6.2, 8.0])
# ax.set_xlim([6.2, 9.5])
plt.grid()
plt.legend(loc=3)
# place a text box in upper left in axes coords
textstr = 'Training features: \n'
for i, label in enumerate(feature_labels):
if (i == len(feature_labels)-1):
textstr += '{}) '.format(i+1) + label
else:
textstr += '{}) '.format(i+1) + label + '\n'
print('textstr = \n' + textstr)
props = dict(facecolor='white')
ax.text(0.65, 0.3, textstr, transform=ax.transAxes, fontsize=8,
verticalalignment='top', bbox=props)
outfile = '/home/jbourbeau/public_html/figures/composition' + \
'/fraction-reco-correct_vs_reco-energy_RF.png'
plt.savefig(outfile)
plt.show()
print(reco_total_frac)
In [9]:
num_features = len(feature_list)
importances = pipeline.named_steps['classifier'].feature_importances_
indices = np.argsort(importances)[::-1]
fig, ax = plt.subplots()
# feature_labels = np.array(['$\\log_{10}({\mathrm{E/GeV})$', 'InIce charge',
# '$\cos(\\theta)$', '$\mathrm{Laputop}\ \chi^2/\mathrm{n.d.f.}$', 'NChannels'])
for f in range(num_features):
print('{}) {}'.format(f + 1, importances[indices[f]]))
plt.ylabel('Feature Importances')
plt.bar(range(num_features),
importances[indices],
align='center')
plt.xticks(range(num_features),
feature_labels[indices], rotation=90)
plt.xlim([-1, len(feature_list)])
plt.ylim([0, .40])
plt.show()
In [ ]: