Commit bc5187d7 authored by JGlombitza's avatar JGlombitza
Browse files

Add DNN's

Add evaluation scripts
New ReadMe
parent 86dafa50
import numpy as np
import glob
import tensorflow as tf
import norm
Folder = './RawData' # Path to simulated shower packages
log_dir = './Data_preprocessed' # Path to preprocessed shower packages
filenames =glob.glob(Folder + "/showers-*")
a = 10000 # packagesize
for i in range(0,len(filenames)):
data = np.load(filenames[i])
if i==0:
signal = data['signal']
time = data['time']
logE = data['logE']
Xmax = data['Xmax']
showeraxis = data['showeraxis']
mass = data['mass']
else:
signal = np.append(signal, data['signal'], axis = 0)
time = np.append(time, data['time'], axis=0)
logE = np.append(logE, data['logE'], axis=0)
Xmax = np.append(Xmax, data['Xmax'], axis = 0)
showeraxis = np.append(showeraxis, data['showeraxis'], axis = 0)
mass = np.append(mass, data['mass'], axis = 0)
# Normalization of arrival times
time = norm.NormCenter(time)
# Preprocessing of signal
signal[signal<0] = 0
SignalSum = norm.NormPhysicalEnergylog10(np.sum(signal, axis=-1))
signal = norm.NormPhysicalTimeTracelog10(signal)
Energy = 10**logE/10**18
# Preprocessed time traces (input1)
Input1 = signal.reshape(len(showeraxis), 9, 9, signal.shape[2], 1)
# Additional preprocessed input maps (arrival times, total signals) (input2)
Input2 = np.stack((time,SignalSum), axis=-1)
# Shuffle data
idx = np.arange(Xmax.shape[0])
np.random.shuffle(idx)
Xmax = Xmax[idx]
showeraxis = showeraxis[idx]
Energy = Energy[idx]
Input1 = Input1[idx]
Input2 = Input2[idx]
mass = mass[idx]
# Store new packages of preprocessed showers
for i in range(len(filenames)):
np.savez_compressed(log_dir+ "/PlanarWaveShower_PrePro_%i" %i, Xmax=Xmax[a*i:a*(i+1)], Energy = Energy[a*i:a*(i+1)], showeraxis=showeraxis[a*i:a*(i+1)], Input1 = Input1[a*i:a*(i+1)], Input2=Input2[a*i:a*(i+1)], mass=mass[a*i:a*(i+1)])
import numpy as np
def NormCenter(time):
u = np.isnan(time)
time -= np.nanmean(time, axis=1)[:,np.newaxis]
time /= np.nanstd(time)
time[u] = 0
try:
time = time.reshape(time.shape[0],11,11)
except:
time = time.reshape(time.shape[0],9,9)
return time
def NormPhysicalEnergylog10(time):
u = np.isnan(time)
time[time<0] = 0
time[~u] = np.log10(time[~u]-np.nanmin(time)+1.0)
time /= np.nanstd(time)
time[u] = 0
try:
time = time.reshape(time.shape[0], 9, 9)
except:
time = time.reshape(time.shape[0], 11,11)
return time
def NormPhysicalTimeTracelog10(signal, log_dir="."):
signal[signal<0] = 0
x = np.arange(80)
signal = np.log10(signal + 1)
signal[np.isnan(signal)] = 0
return signal
import numpy as np
from airshower import shower, utils, plotting
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--mass", default=1, type=int, help='mass number')
parser.add_argument("--nfiles", default=10, type=int, help='number of files')
parser.add_argument("--nevents", default=10000, type=int, help='number of events per file')
args = parser.parse_args()
print 'Simulating'
print 'A', args.mass
print 'nfiles', args.nfiles
print 'nevents', args.nevents
v_stations = utils.station_coordinates(9, layout='offset')
for i in range(args.nfiles):
logE = 18.5 + 1.5 * np.random.rand(args.nevents)
mass = args.mass * np.ones(args.nevents)
data = shower.rand_events(logE, mass, v_stations)
print "save package:", i
np.savez_compressed(
'./RawData/showers-A%i-%i.npz' % (args.mass, i),
detector=data['detector'],
logE=data['logE'],
mass=data['mass'],
Xmax=data['Xmax'],
showercore=data['showercore'],
showeraxis=data['showeraxis'],
showermax=data['showermax'],
time=data['time'],
signal=data['signal'])
import numpy as np
import keras
import glob
import tensorflow as tf
from keras.callbacks import TensorBoard, EarlyStopping, ReduceLROnPlateau#, NanChecker
from keras.models import Sequential, Model
from keras.layers import *
from keras.optimizers import Adam
from keras.utils import plot_model
import tools
max_steps=15
lr = 0.001
log_dir="."
nbatch = 132
# ----------------------------------------------------------------------
# load and prepare data
# ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*")
N = 10000*len(filenames)
a = 10000 # packagesize
Input1 = np.zeros(N*9*9*80).reshape(N,9,9,80,1)
Input2 = np.zeros(N*9*9*2).reshape(N,9,9,2)
showeraxis = np.zeros(N*3).reshape(N,3)
Energy = np.zeros(N)
for i in range(0,len(filenames)):
data = np.load(filenames[i])
Input1[a*i:a*(i+1)] = data['Input1'].reshape(a,9,9,80,1)
showeraxis[a*i:a*(i+1)] = data['showeraxis'].reshape(a,3)
Input2[a*i:a*(i+1)] = data['Input2'].reshape(a,9,9,2)
Energy[a*i:a*(i+1)] = data['Energy'].reshape(a)
print('Loaded %i showers' %N )
### Reshape, Split to test and train sets
ntest = 40000
nval = 2000.
nsplit = N - ntest
Input1, x_test = Input1[:-ntest], Input1[-ntest:]
x_train = Input1
x_train_t0, x_test_t0 = Input2[:-ntest], Input2[-ntest:]
y_train, y_test = showeraxis[:-ntest], showeraxis[-ntest:]
## Callbacks & Optimizers
Adam = Adam(lr=lr, beta_1=0.9, beta_2=0.999, epsilon=1e-08, decay=0.0, clipnorm = 5)
reduceLearning = ReduceLROnPlateau(monitor='val_MeanElong', factor=0.6, patience=17, verbose=1, mode='min', epsilon=0.04, cooldown=0, min_lr=0.0)
# ----------------------------------------------------------------------
# Model
# ----------------------------------------------------------------------
TimeProcessingNodes = 10
kwargs = dict(activation='relu', kernel_initializer='he_normal')
inputTimeTraces = Input(shape=(x_train.shape[1], x_train.shape[2], x_train.shape[3], 1), name='TraceInput') # Input Cube aus time traces
inputT0 = Input(shape=(x_train_t0.shape[1], x_train_t0.shape[2], x_train_t0.shape[3]), name='T0Input') # T0 input
## Time trace characterization
TimeProcess = Conv3D(64, (1,1,7), padding='valid', activation='relu', strides=(1,1,4))(inputTimeTraces)
TimeProcess = Conv3D(32, (1,1,7), padding='valid', activation='relu', strides=(1,1,4))(TimeProcess)
TimeProcess = Conv3D(10, (1,1,4), padding='valid', activation='relu', strides=(1,1,1))(TimeProcess)
TimeProcess = Reshape((9,9,TimeProcessingNodes))(TimeProcess)
x = concatenate([TimeProcess, inputT0])
# Densely connected convolutions
nfilter=12
for i in range(4):
x = tools.DenselyConnectedSepConv(x, nfilter, **kwargs)
nfilter = 2*nfilter
x = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(x)
x = Flatten()(x)
output = Dense(3)(x)
model = Model(inputs=[inputTimeTraces, inputT0], outputs=output)
# ----------------------------------------------------------------------
model.compile(loss='mse', optimizer=Adam, metrics=[tools.MeanElong])
fit = model.fit([x_train, x_train_t0],[y_train], epochs=max_steps, batch_size=nbatch, verbose=1, callbacks=[reduceLearning], validation_split = nval/nsplit)
#------------#
# Evaluation #
#------------#
# Test evaluation
y_pred = model.predict([x_test, x_test_t0])
x_test = x_test.reshape(x_test.shape[0], x_test.shape[1], x_test.shape[2], x_test.shape[3])
angulardistance = tools.PlotHistosAngularReconstruction(Signal = np.sum(x_test, axis=-1), y_true=y_test, y_pred=y_pred, log_dir=log_dir, name="test", Energy = Energy[-ntest:])
## Train evaluation
#y_pred = model.predict([x_train, x_train_t0])
#x_train = x_train.reshape(x_train.shape[0], x_train.shape[1], x_train.shape[2], x_train.shape[3])
#angulardistance = tools.PlotHistosAngularReconstruction(Signal = np.sum(x_train, axis=-1), y_true=y_train, y_pred=y_pred, log_dir=log_dir, name="train", Energy = Energy[:-ntest])
import numpy as np
import keras
import glob
import tensorflow as tf
from keras.callbacks import TensorBoard, EarlyStopping, ReduceLROnPlateau#, NanChecker
from keras.models import Sequential, Model
from keras.layers import *
from keras.optimizers import Adam
import tools
max_steps=15
lr = 0.001
log_dir="."
nbatch = 132
# ----------------------------------------------------------------------
# load and prepare data
# ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*")
N = 10000*len(filenames)
a = 10000 # packagesize
Input2 = np.zeros(N*9*9*2).reshape(N,9,9,2)
showeraxis = np.zeros(N*3).reshape(N,3)
Energy = np.zeros(N)
for i in range(0,len(filenames)):
data = np.load(filenames[i])
showeraxis[a*i:a*(i+1)] = data['showeraxis'].reshape(a,3)
Input2[a*i:a*(i+1)] = data['Input2'][:,:,:,0].reshape(a,9,9,1)
Energy[a*i:a*(i+1)] = data['Energy'].reshape(a)
## Reshape, Split to test and train sets
ntest = 40000
nval = 2000.
nsplit = N - ntest
x_train, x_test = np.split(Input2, [nsplit])
y_train, y_test = showeraxis[:nsplit], showeraxis[nsplit:]
## Callbacks & Optimizers
Adam = Adam(lr=lr, beta_1=0.9, beta_2=0.999, epsilon=1e-08, decay=0.0, clipnorm = 2)
reduceLearning = ReduceLROnPlateau(monitor='val_loss', factor=0.6, patience=17, verbose=1, mode='min', epsilon=0.07, cooldown=0, min_lr=0.0)
# ----------------------------------------------------------------------
# Model
# ----------------------------------------------------------------------
input1 = Input(shape=(Input2.shape[1], Input2.shape[2], Input2.shape[3]))
kwargs = dict(activation='relu', kernel_initializer='he_normal')
x = Conv2D(16, (3, 3), padding = 'same', **kwargs)(input1)
# Densely connected convolutions
nfilter=16
for i in range(4):
x = tools.DenselyConnectedSepConv(x, nfilter, **kwargs)
nfilter = 2*nfilter
x = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(x)
x = Flatten()(x)
output = Dense(3)(x)
model = Model([input1], [output])
# ----------------------------------------------------------------------
model.compile(loss='mse', optimizer=Adam, metrics=[tools.MeanElong])
fit = model.fit(x_train, y_train, batch_size=nbatch, shuffle=True, epochs = max_steps, verbose=1, validation_split=nval/nsplit, callbacks=[reduceLearning])
#------------#
# Evaluation #
#------------#
# Test evaluation
y_pred = model.predict(x_test)
angulardistance = tools.PlotHistosAngularReconstruction(y_true=y_test, y_pred=y_pred, log_dir=log_dir, name="OnlyTime_test", Energy = Energy[-ntest:])
## Train evaluation
#y_pred = model.predict(x_train)
#angulardistance = tools.PlotHistosAngularReconstruction(y_true=y_train, y_pred=y_pred, log_dir=log_dir, name="OnlyTime_train", Energy = Energy[:-ntest])
import numpy as np
import keras
import glob
import tensorflow as tf
from keras.callbacks import TensorBoard, EarlyStopping, ReduceLROnPlateau#, NanChecker
from keras.models import Sequential, Model
from keras.layers import *
from keras.optimizers import Adam
import tools
max_steps=15
lr = 0.001
log_dir="."
nbatch = 132
def DenselyConnectedSepConv(z, nfilter, **kwargs):
c = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(z)
return concatenate([z, c], axis=-1)
# ----------------------------------------------------------------------
# load and prepare data
# ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*")
N = 10000*len(filenames)
a = 10000 # packagesize
Input1 = np.zeros(N*9*9*80).reshape(N,9,9,80,1)
Input2 = np.zeros(N*9*9*2).reshape(N,9,9,2)
Energy = np.zeros(N)
for i in range(0,len(filenames)):
data = np.load(filenames[i])
Input1[a*i:a*(i+1)] = data['Input1'].reshape(a,9,9,80,1)
Energy[a*i:a*(i+1)] = data['Energy']
Input2[a*i:a*(i+1)] = data['Input2'].reshape(a,9,9,2)
## Reshape, Split to test and train sets
ntest = 40000
nval = 2000.
nsplit = N - ntest
Input1, x_test = Input1[:-ntest], Input1[-ntest:]
x_train = Input1
x_train_t0, x_test_t0 = Input2[:-ntest], Input2[-ntest:]
y_train, y_test = Energy[:-ntest], Energy[-ntest:]
## Callbacks & Optimizers
Adam = Adam(lr=lr, beta_1=0.9, beta_2=0.999, epsilon=1e-08, decay=0.0, clipnorm = 5)
reduceLearning = ReduceLROnPlateau(monitor='val_Distance', factor=0.4, patience=17, verbose=1, mode='min', epsilon=1, cooldown=0, min_lr=0.0)
# ----------------------------------------------------------------------
# Model
# ----------------------------------------------------------------------
kwargs = dict(activation='relu', kernel_initializer='he_normal')
inputTimeTraces = Input(shape=(x_train.shape[1], x_train.shape[2], x_train.shape[3], 1), name='TraceInput') # Input Cube aus time traces
inputT0 = Input(shape=(x_train_t0.shape[1], x_train_t0.shape[2], x_train_t0.shape[3]), name='T0Input') # T0 input
## Time trace characterization
TimeProcess = Conv3D(64, (1,1,7), padding='valid', activation='relu', strides=(1,1,4))(inputTimeTraces)
TimeProcess = Conv3D(32, (1,1,7), padding='valid', activation='relu', strides=(1,1,4))(TimeProcess)
TimeProcess = Conv3D(10, (1,1,4), padding='valid', activation='relu', strides=(1,1,1))(TimeProcess)
TimeProcess = Reshape((9,9,10))(TimeProcess)
x = concatenate([TimeProcess, inputT0])
# Densely connected convolutions
nfilter =12
for i in range(4):
x = DenselyConnectedSepConv(x, nfilter, **kwargs)
nfilter = 2*nfilter
x = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(x)
x = Flatten()(x)
output = Dense(1)(x)
model = Model(inputs=[inputTimeTraces, inputT0], outputs=output)
# ----------------------------------------------------------------------
model.compile(loss='mse', optimizer=Adam, metrics=[tools.Distance])
fit = model.fit([x_train, x_train_t0],[y_train], epochs=max_steps, batch_size=nbatch, verbose=1, callbacks=[reduceLearning], validation_split = nval/nsplit)
#------------#
# Evaluation #
#------------#
y_pred = model.predict([x_test, x_test_t0])
# Test evaluation
reco = tools.PlotHistosEnergyReconstruction(y_true=y_test, y_pred=y_pred, log_dir=log_dir, name="test")
## Train evaluation
#y_pred = model.predict([x_train, x_train_t0])
#tools.PlotHistosEnergyReconstruction(y_true=y_train, y_pred=y_pred, log_dir=log_dir, name="train")
import numpy as np
import keras
import glob
import tensorflow as tf
from keras.callbacks import TensorBoard, EarlyStopping, ReduceLROnPlateau#, NanChecker
from keras.models import Sequential, Model
from keras.layers import *
from keras.optimizers import Adam
import tools
max_steps=15
lr = 0.001
log_dir="."
nbatch = 132
# ----------------------------------------------------------------------
# load and prepare data
# ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*")
N = 10000*len(filenames)
a = 10000 # packagesize
Input2 = np.zeros(N*9*9*2).reshape(N,9,9,2)
Energy = np.zeros(N)
for i in range(0,len(filenames)):
data = np.load(filenames[i])
Energy[a*i:a*(i+1)] = data['Energy']
Input2[a*i:a*(i+1)] = data['Input2'][:,:,:,1].reshape(a, 9, 9, 1)
print('Loaded %i showers' %N )
## Split to test and train sets
ntest = 40000
nval = 2000.
nsplit = N - ntest
x_train, x_test = np.split(Input2, [nsplit])
y_train, y_test = np.split(Energy, [nsplit])
## Callbacks & Optimizers
Adam = Adam(lr=lr, beta_1=0.9, beta_2=0.999, epsilon=1e-08, decay=0.0, clipnorm = 2)
reduceLearning = ReduceLROnPlateau(monitor='val_Distance', factor=0.4, patience=17, verbose=1, mode='min', epsilon=1, cooldown=0, min_lr=0.0)
# ----------------------------------------------------------------------
# Model
# ----------------------------------------------------------------------
input1 = Input(shape=Input2.shape[1:])
kwargs = dict(activation='relu', kernel_initializer='he_normal')
x = Conv2D(16, (3, 3), padding = 'same', **kwargs)(input1)
# Densely connected convolutions
nfilter=16
for i in range(4):
x = tools.DenselyConnectedSepConv(x, nfilter)
nfilter = 2*nfilter
x = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(x)
x = Flatten()(x)
output = Dense(1)(x)
model = Model([input1], [output])
# ----------------------------------------------------------------------
model.compile(loss='mse', optimizer=Adam, metrics=[tools.Distance])
fit = model.fit(x_train, y_train, batch_size=nbatch, shuffle=True, epochs = max_steps, verbose=1, validation_split=nval/nsplit, callbacks=[reduceLearning])
#------------#
# Evaluation #
#------------#
y_pred = model.predict(x_test)
# Test evaluation
reco = tools.PlotHistosEnergyReconstruction(y_true=y_test, y_pred=y_pred, log_dir=log_dir, name="OnlySignal_test")
## Train evaluation
#y_pred = model.predict(x_train)
#tools.PlotHistosEnergyReconstruction(y_true=y_train, y_pred=y_pred, log_dir=log_dir, name="OnlySignal_train")
import numpy as np
import keras
import glob
import tensorflow as tf
from keras.callbacks import TensorBoard, EarlyStopping, ReduceLROnPlateau#, NanChecker
from keras.models import Sequential, Model
from keras.layers import *
from keras.optimizers import Adam
import tools
max_steps=15
lr = 0.001
log_dir="."
nbatch = 132
# ----------------------------------------------------------------------
# load and prepare data
# ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*")
N = 10000*len(filenames)
a = 10000 # packagesize
Input1 = np.zeros(N*9*9*80).reshape(N,9,9,80,1)
Input2 = np.zeros(N*9*9*2).reshape(N,9,9,2)
Xmax = np.zeros(N)
Energy = np.zeros(N)
for i in range(0,len(filenames)):
data = np.load(filenames[i])
Input1[a*i:a*(i+1)] = data['Input1'].reshape(a,9,9,80,1)
Xmax[a*i:a*(i+1)] = data['Xmax']
Energy[a*i:a*(i+1)] = data['Energy']
Input2[a*i:a*(i+1)] = data['Input2'].reshape(a,9,9,2)
## Reshape, Split to test and train sets
ntest = 40000
nval = 2000.
nsplit = N - ntest
Input1, x_test = Input1[:-ntest], Input1[-ntest:]
x_train = Input1
x_train_t0, x_test_t0 = Input2[:-ntest], Input2[-ntest:]
y_train, y_test = Xmax[:-ntest], Xmax[-ntest:]
## Callbacks & Optimizers
Adam = Adam(lr=lr, beta_1=0.9, beta_2=0.999, epsilon=1e-08, decay=0.0, clipnorm = 5)
reduceLearning = ReduceLROnPlateau(monitor='val_Distance', factor=0.4, patience=17, verbose=1, mode='min', epsilon=4, cooldown=0, min_lr=0.0)
def DenselyConnectedSepConv(z, nfilter, **kwargs):
c = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(z)
return concatenate([z, c], axis=-1)
# ----------------------------------------------------------------------
# Model
# ----------------------------------------------------------------------
kwargs = dict(activation='relu', kernel_initializer='he_normal')
inputTimeTraces = Input(shape=(x_train.shape[1], x_train.shape[2], x_train.shape[3], 1), name='TraceInput') # Input Cube aus time traces
inputT0 = Input(shape=(x_train_t0.shape[1], x_train_t0.shape[2], x_train_t0.shape[3]), name='T0Input') # T0 input
## Time trace characterization
TimeProcess = Conv3D(64, (1,1,7), padding='valid', activation='relu', strides=(1,1,4))(inputTimeTraces)
TimeProcess = Conv3D(32, (1,1,7), padding='valid', activation='relu', strides=(1,1,4))(TimeProcess)
TimeProcess = Conv3D(10, (1,1,4), padding='valid', activation='relu', strides=(1,1,1))(TimeProcess)
TimeProcess = Reshape((9,9,10))(TimeProcess)
x = concatenate([TimeProcess, inputT0])
# Densely connected convolutions
nfilter=12
for i in range(4):
x = DenselyConnectedSepConv(x, nfilter, **kwargs)
nfilter = 2*nfilter
x = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(x)
x = Flatten()(x)
output = Dense(1)(x)
model = Model(inputs=[inputTimeTraces, inputT0], outputs=output)
# ----------------------------------------------------------------------
model.compile(loss='mse', optimizer=Adam, metrics=[tools.Distance])
fit = model.fit([x_train, x_train_t0],[y_train], epochs=max_steps, batch_size=nbatch, verbose=1, callbacks=[reduceLearning], validation_split = nval/nsplit)
#------------#
# Evaluation #
#------------#
# Test evaluation
y_pred = model.predict([x_test, x_test_t0])
reco = tools.PlotHistosXmaxReconstruction(y_true=Xmax[-ntest:], y_pred=y_pred, log_dir=log_dir, name="test", Energy=Energy[-ntest:])
## Train evaluation
#y_pred = model.predict([x_train, x_train_t0])
#tools.PlotHistosXmaxReconstruction(y_true=Xmax[:-ntest], y_pred=y_pred, log_dir=log_dir, name="train", Energy=Energy[:-ntest])