diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/data_generation.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/data_generation.py index fbb47f1a54..418e6276ca 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/data_generation.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/data_generation.py @@ -241,7 +241,7 @@ def generate_data( data['inputs'].append(data_run[:input_width]) data['labels'].append(data_run[input_width:]) data['contact_matrix'].append(np.array(damped_contact_matrix)) - data['damping_day'].append(damping_day) + data['damping_day'].append([damping_day]) bar.next() bar.finish() @@ -271,13 +271,13 @@ def getBaselineMatrix(): """ loads the baselinematrix""" baseline_contact_matrix0 = os.path.join( - "./data/contacts/baseline_home.txt") + "./data/Germany/contacts/baseline_home.txt") baseline_contact_matrix1 = os.path.join( - "./data/contacts/baseline_school_pf_eig.txt") + "./data/Germany/contacts/baseline_school_pf_eig.txt") baseline_contact_matrix2 = os.path.join( - "./data/contacts/baseline_work.txt") + "./data/Germany/contacts/baseline_work.txt") baseline_contact_matrix3 = os.path.join( - "./data/contacts/baseline_other.txt") + "./data/Germany/contacts/baseline_other.txt") baseline = np.loadtxt(baseline_contact_matrix0) \ + np.loadtxt(baseline_contact_matrix1) + \ @@ -329,7 +329,7 @@ def get_population(path): os.path.dirname(os.path.realpath(path)))), 'data') path_population = os.path.abspath( - r"data//pydata//Germany//county_population.json") + r"data//Germany//pydata//county_current_population.json") input_width = 5 label_width = 30 diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/model.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/model.py index 81023326f0..221291f36f 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/model.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/model.py @@ -42,7 +42,7 @@ def plot_compartment_prediction_model( :param plot_col: string name of compartment to be plotted. :param max_subplots: Number of the simulation runs to be plotted and compared against. (Default value = 8) :param plot_compartment: (Default value = 'InfectedSymptoms') - + :returns: No return """ num_groups = 6 num_compartments = 8 @@ -137,166 +137,278 @@ def plot_compartment_prediction_model( plt.savefig('plots/evaluation_secir_groups_' + plot_compartment + '.png') +#################### +# Helper functions # +#################### +def calc_split_index(n, split_train=0.7, + split_valid=0.2, split_test=0.1): + """ + Calculating the indixes for a split_train:split_valid:split_test decomposition of a set with size n + + It must hold split_train + split_valid + split_test = 1 + + :param n: integer value + :param split_train: value between 0 and 1 + :param split_valid: value between 0 and 1 + :param split_test: value between 0 and 1 + :returns: a list of the form [i_train, i_valid, i_test] + """ + if split_train + split_valid + split_test > 1 + 1e-10: + raise ValueError( + "Summed data set shares are greater than 1. Please adjust the values.") + n_train = int(n * split_train) + n_valid = int(n * split_valid) + n_test = n - n_train - n_valid + + return [n_train, n_valid, n_test] + + +def flat_input(input): + """ Flatten input dimension + + :param input: input array of size (n,k,l) + :returns: reshaped array of size (n, k*l) + + """ + dim = tf.reduce_prod(tf.shape(input)[1:]) + return tf.reshape(input, [-1, dim]) + + +def prepare_data_classic(data): + """ + Transforming data to be processable by "classic" network, simply flattening and concatenating for each data instance. + + :param data: dictionary produces by data_generation + :returns: dictionary with entries { + "train_inputs", "train_labels", "valid_inputs", + "valid_labels", "test_inputs", "test_labels"} + """ + # Getting number of samples + n = data["inputs"].shape[0] + + # Calculate split inidces + split_indices = calc_split_index(n) + + # Flattening all inputs + inputs = flat_input(data["inputs"]) + labels = data["labels"] + cmatrix = flat_input(data["contact_matrix"]) + dampdays = data["damping_day"] + + # Splitting the data + compinputs_train, compinputs_valid, compinputs_test = tf.split( + inputs, split_indices, 0) + labels_train, labels_valid, labels_test = tf.split( + labels, split_indices, 0) + cmatrix_train, cmatrix_valid, cmatrix_test = tf.split( + cmatrix, split_indices, 0) + dampdays_train, dampdays_valid, dampdays_test = tf.split( + dampdays, split_indices, 0) + + # Combining the ingredients to one input object + inputs_train = tf.concat( + [tf.cast(compinputs_train, tf.float32), + tf.cast(cmatrix_train, tf.float32), + tf.cast(dampdays_train, tf.float32)], + axis=1, name='concat') + inputs_valid = tf.concat( + [tf.cast(compinputs_valid, tf.float32), + tf.cast(cmatrix_valid, tf.float32), + tf.cast(dampdays_valid, tf.float32)], + axis=1, name='concat') + inputs_test = tf.concat( + [tf.cast(compinputs_test, tf.float32), + tf.cast(cmatrix_test, tf.float32), + tf.cast(dampdays_test, tf.float32)], + axis=1, name='concat') + + return { + "train_inputs": inputs_train, + "train_labels": labels_train, + "valid_inputs": inputs_valid, + "valid_labels": labels_valid, + "test_inputs": inputs_test, + "test_labels": labels_test + } + + +def prod_time_series(obj, n, length_input): + """ + Repeating static informations to fit into a time series framework + + :param obj: an array of objects, which should be repeated + :param n: total number of samples + :param length_input: number of days observed per input + :returns: a tensor of shape [n, length_input, -1], where for each sample the static object is repeated length_input times + """ + new_obj = [] + for i in obj: + new_obj.extend([i for _ in range(length_input)]) + + new_obj = tf.reshape( + tf.stack(new_obj), + [n, length_input, -1]) + return new_obj + + +def prepare_data_timeseries(data): + """ + Transforming data to be processable by "time_series" network, simply repeating static values, flattening and concatenating for each data instance. + + :param data: dictionary produces by data_generation + :returns: dictionary with entries { + "train_inputs", "train_labels", "valid_inputs", + "valid_labels", "test_inputs", "test_labels" + } + """ + # Getting the number of samples + n = data["inputs"].shape[0] + + # number of days per input sample + input_width = 5 + + # Reshaping the matrix input + cmatrix = flat_input(tf.stack(data["contact_matrix"])) + + # Repeat data (contact matrix and dampinhg day) to produce time series + cmatrix_repeated = prod_time_series(cmatrix, n, input_width) + dampdays_repeated = prod_time_series(data["damping_day"], n, input_width) + + # Calculate split inidces + split_indices = calc_split_index(n) + + # Splitting the data + compinputs_train, compinputs_valid, compinputs_test = tf.split( + data["inputs"], split_indices, 0) + labels_train, labels_valid, labels_test = tf.split( + data["labels"], split_indices, 0) + cmatrix_train, cmatrix_valid, cmatrix_test = tf.split( + cmatrix_repeated, split_indices, 0) + dampdays_train, dampdays_valid, dampdays_test = tf.split( + dampdays_repeated, split_indices, 0) + + # Combining the ingredients to one input object + inputs_train = tf.concat( + [tf.cast(compinputs_train, tf.float16), + tf.cast(cmatrix_train, tf.float16), + tf.cast(dampdays_train, tf.float16)], + axis=2, name='concat') + inputs_valid = tf.concat( + [tf.cast(compinputs_valid, tf.float16), + tf.cast(cmatrix_valid, tf.float16), + tf.cast(dampdays_valid, tf.float16)], + axis=2, name='concat') + inputs_test = tf.concat( + [tf.cast(compinputs_test, tf.float16), + tf.cast(cmatrix_test, tf.float16), + tf.cast(dampdays_test, tf.float16)], + axis=2, name='concat') + + return { + "train_inputs": inputs_train, + "train_labels": labels_train, + "valid_inputs": inputs_valid, + "valid_labels": labels_valid, + "test_inputs": inputs_test, + "test_labels": labels_test + } + + +######################################### +# Initialization and Training of Models # +######################################### + +def initialize_model(parameters): + """ Initialize model from given list of parameters + + :param parameters: tuple of parameters describing the model architecture, it should be of the form + (number_of_output_days, number_age_groups, number_compartments, + hidden_layers, neurons_in_hidden_layer, activation_function, modelname) + :returns: tensor flow keras model with the given architecture, see network_architectures.py + """ + label_width, number_age_groups, number_compartments, hidden_layers, neurons_in_hidden_layer, activation_function, modelname = parameters + + if modelname == "Dense": + return network_architectures.mlp_multi_input_multi_output( + label_width=label_width, + num_age_groups=number_age_groups, + num_outputs=number_compartments, + num_hidden_layers=hidden_layers, + num_neurons_per_layer=neurons_in_hidden_layer, + activation=activation_function + ) + elif modelname == "LSTM": + return network_architectures.lstm_multi_input_multi_output( + label_width=label_width, + num_age_groups=number_age_groups, + num_outputs=number_compartments, + num_hidden_layers=hidden_layers, + num_neurons_per_layer=neurons_in_hidden_layer, + activation=activation_function + ) + elif modelname == "CNN": + return network_architectures.cnn_multi_input_multi_output( + label_width=label_width, + num_age_groups=number_age_groups, + num_outputs=number_compartments, + num_hidden_layers=hidden_layers, + num_neurons_per_layer=neurons_in_hidden_layer, + activation=activation_function + ) + else: + raise ValueError( + "name_architecture must be one of 'Dense', 'LSTM' or 'CNN'" + ) + + def network_fit( - path, model, modeltype, max_epochs=30, early_stop=500, plot=True): + model, modeltype, training_parameter, path, filename='data_secir_groups_30days_Germany_10k_damp.pickle', plot=True): """ Training and evaluation of a given model with mean squared error loss and Adam optimizer using the mean absolute error as a metric. - :param path: path of the dataset. :param model: Keras sequential model. :param modeltype: type of model. Can be 'classic' or 'timeseries'. Data preparation is made based on the modeltype. - :param max_epochs: int maximum number of epochs in training. (Default value = 30) - :param early_stop: Integer that forces an early stop of training if the given number of epochs does not give a significant reduction of validation loss. (Default value = 500) + :param training_parameter: tuple of parameters used for the training process, it should be of the form + (early_stop, max_epochs, loss, optimizer, metrics), where loss is a loss-function implemented in keras, optimizer is the name of the used optimizer, + metrics is a list of used training metrics, e.g. [tf.keras.metrics.MeanAbsoluteError(), tf.keras.metrics.MeanAbsolutePercentageError()] + :param path: path of the dataset. + :param filename: name of the file containing the data :param plot: (Default value = True) - + :returns: training history as returned by the keras fit() method. """ + # Unpacking training parameters + early_stop, max_epochs, loss, optimizer, metrics = training_parameter - if not os.path.isfile(os.path.join(path, 'data_secir_groups.pickle')): + # Getting data and loading it + if not os.path.isfile(os.path.join(path, filename)): ValueError("no dataset found in path: " + path) - - file = open(os.path.join(path, 'data_secir_groups.pickle'), 'rb') + file = open(os.path.join( + path, filename), 'rb') data = pickle.load(file) - data_splitted = split_data(data['inputs'], data['labels']) + # preprocessing the data if modeltype == 'classic': - - train_inputs_compartments = flat_input(data_splitted["train_inputs"]) - train_labels = (data_splitted["train_labels"]) - valid_inputs_compartments = flat_input(data_splitted["valid_inputs"]) - valid_labels = (data_splitted["valid_labels"]) - test_inputs_compartments = flat_input(data_splitted["test_inputs"]) - test_labels = (data_splitted["test_labels"]) - - contact_matrices = split_contact_matrices( - tf.stack(data["contact_matrix"])) - contact_matrices_train = flat_input(contact_matrices['train']) - contact_matrices_valid = flat_input(contact_matrices['valid']) - contact_matrices_test = flat_input(contact_matrices['test']) - - damping_days = data['damping_day'] - damping_days_splitted = split_damping_days(damping_days) - damping_days_train = damping_days_splitted['train'] - damping_days_valid = damping_days_splitted['valid'] - damping_days_test = damping_days_splitted['test'] - - train_inputs = tf.concat( - [tf.cast(train_inputs_compartments, tf.float32), - tf.cast(contact_matrices_train, tf.float32), - tf.cast(damping_days_train, tf.float32)], - axis=1, name='concat') - valid_inputs = tf.concat( - [tf.cast(valid_inputs_compartments, tf.float32), - tf.cast(contact_matrices_valid, tf.float32), - tf.cast(damping_days_valid, tf.float32)], - axis=1, name='concat') - test_inputs = tf.concat( - [tf.cast(test_inputs_compartments, tf.float32), - tf.cast(contact_matrices_test, tf.float32), - tf.cast(damping_days_test, tf.float32)], - axis=1, name='concat') + data_prep = prepare_data_classic(data) elif modeltype == 'timeseries': + data_prep = prepare_data_timeseries(data) - train_inputs_compartments = (data_splitted["train_inputs"]) - train_labels = (data_splitted["train_labels"]) - valid_inputs_compartments = (data_splitted["valid_inputs"]) - valid_labels = (data_splitted["valid_labels"]) - test_inputs_compartments = (data_splitted["test_inputs"]) - test_labels = (data_splitted["test_labels"]) - - contact_matrices = split_contact_matrices( - tf.stack(data["contact_matrix"])) - contact_matrices_train = flat_input(contact_matrices['train']) - contact_matrices_valid = flat_input(contact_matrices['valid']) - contact_matrices_test = flat_input(contact_matrices['test']) - - n = np.array(data['damping_day']).shape[0] - train_days = data['damping_day'][:int(n*0.7)] - valid_days = data['damping_day'][int(n*0.7):int(n*0.9)] - test_days = data['damping_day'][int(n*0.9):] - - # concatenate the compartment data with contact matrices and damping days - # to receive complete input data - new_contact_train = [] - for i in contact_matrices_train: - new_contact_train.extend([i for j in range(5)]) - - new_contact_train = tf.reshape( - tf.stack(new_contact_train), - [train_inputs_compartments.shape[0], - 5, np.asarray(new_contact_train).shape[1]]) - - new_damping_days_train = [] - for i in train_days: - new_damping_days_train.extend([i for j in range(5)]) - new_damping_days_train = tf.reshape( - tf.stack(new_damping_days_train), - [train_inputs_compartments.shape[0], - 5, 1]) - - train_inputs = tf.concat( - (tf.cast(train_inputs_compartments, tf.float16), - tf.cast(new_contact_train, tf.float16), - tf.cast(new_damping_days_train, tf.float16)), - axis=2) - - new_contact_test = [] - for i in contact_matrices_test: - new_contact_test.extend([i for j in range(5)]) - - new_contact_test = tf.reshape(tf.stack(new_contact_test), [ - contact_matrices_test.shape[0], 5, contact_matrices_test.shape[1]]) - - new_damping_days_test = [] - for i in test_days: - new_damping_days_test.extend([i for j in range(5)]) - new_damping_days_test = tf.reshape( - tf.stack(new_damping_days_test), - [test_inputs_compartments.shape[0], - 5, 1]) - - test_inputs = tf.concat( - (tf.cast(test_inputs_compartments, tf.float16), - tf.cast(new_contact_test, tf.float16), - tf.cast(new_damping_days_test, tf.float16)), - axis=2) - - new_contact_val = [] - for i in contact_matrices_valid: - new_contact_val.extend([i for j in range(5)]) - - new_contact_val = tf.reshape( - tf.stack(new_contact_val), - [contact_matrices_valid.shape[0], - 5, contact_matrices_valid.shape[1]]) - - new_damping_days_valid = [] - for i in valid_days: - new_damping_days_valid.extend([i for j in range(5)]) - new_damping_days_valid = tf.reshape( - tf.stack(new_damping_days_valid), - [valid_inputs_compartments.shape[0], - 5, 1]) - - valid_inputs = tf.concat( - (tf.cast(valid_inputs_compartments, tf.float16), - tf.cast(new_contact_val, tf.float16), - tf.cast(new_damping_days_valid, tf.float16)), - axis=2) - + # Setting up the training parameters batch_size = 32 - early_stopping = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=early_stop, mode='min') - model.compile( - loss=tf.keras.losses.MeanAbsolutePercentageError(), - optimizer=tf.keras.optimizers.Adam(learning_rate=0.001), - metrics=[tf.keras.metrics.MeanSquaredError()]) - - history = model.fit(train_inputs, train_labels, epochs=max_epochs, - validation_data=(valid_inputs, valid_labels), + loss=loss, + optimizer=optimizer, + metrics=metrics) + + history = model.fit(data_prep["train_inputs"], + data_prep["train_labels"], + epochs=max_epochs, + validation_data=( + data_prep["valid_inputs"], + data_prep["valid_labels"]), batch_size=batch_size, callbacks=[early_stopping]) @@ -310,11 +422,15 @@ def network_fit( return history +##################### +# Plots etc. +##################### + def plot_losses(history): """ Plots the losses of the model training. :param history: model training history. - + :returns: No return """ plt.plot(history.history['loss']) plt.plot(history.history['val_loss']) @@ -334,7 +450,7 @@ def get_test_statistic(test_inputs, test_labels, model): :param test_inputs: inputs from test data. :param test_labels: labels (output) from test data. :param model: trained model. - + :returns: dataframe containing the MAPE for the different compartments """ pred = model(test_inputs) @@ -358,120 +474,8 @@ def get_test_statistic(test_inputs, test_labels, model): return mean_percentage -def split_data(inputs, labels, split_train=0.7, - split_valid=0.2, split_test=0.1): - """ Split data set in training, validation and testing data sets. - - :param inputs: input dataset - :param labels: label dataset - :param split_train: Share of training data sets. (Default value = 0.7) - :param split_valid: Share of validation data sets. (Default value = 0.2) - :param split_test: Share of testing data sets. (Default value = 0.1) - - """ - - if split_train + split_valid + split_test > 1 + 1e-10: - raise ValueError( - "Summed data set shares are greater than 1. Please adjust the values.") - elif inputs.shape[0] != labels.shape[0] or inputs.shape[2] != labels.shape[2]: - raise ValueError( - "Number of batches or features different for input and labels") - - n = inputs.shape[0] - n_train = int(n * split_train) - n_valid = int(n * split_valid) - n_test = n - n_train - n_valid - - inputs_train, inputs_valid, inputs_test = tf.split( - inputs, [n_train, n_valid, n_test], 0) - labels_train, labels_valid, labels_test = tf.split( - labels, [n_train, n_valid, n_test], 0) - - data = { - 'train_inputs': inputs_train, - 'train_labels': labels_train, - 'valid_inputs': inputs_valid, - 'valid_labels': labels_valid, - 'test_inputs': inputs_test, - 'test_labels': labels_test - } - - return data - - -def flat_input(input): - """ Flatten input dimension - - :param input: input array - - """ - dim = tf.reduce_prod(tf.shape(input)[1:]) - return tf.reshape(input, [-1, dim]) - - -def split_contact_matrices(contact_matrices, split_train=0.7, - split_valid=0.2, split_test=0.1): - """ Split dampings in train, valid and test - - :param contact_matrices: contact matrices - :param labels: label dataset - :param split_train: ratio of train datasets (Default value = 0.7) - :param split_valid: ratio of validation datasets (Default value = 0.2) - :param split_test: ratio of test datasets (Default value = 0.1) - - """ - - if split_train + split_valid + split_test != 1: - ValueError("summed Split ratios not equal 1! Please adjust the values") - - n = contact_matrices.shape[0] - n_train = int(n * split_train) - n_valid = int(n * split_valid) - n_test = n - n_train - n_valid - - contact_matrices_train, contact_matrices_valid, contact_matrices_test = tf.split( - contact_matrices, [n_train, n_valid, n_test], 0) - data = { - "train": contact_matrices_train, - "valid": contact_matrices_valid, - "test": contact_matrices_test - } - - return data - - -def split_damping_days(damping_days, split_train=0.7, - split_valid=0.2, split_test=0.1): - """ Split damping days in train, valid and test - - :param damping_days: damping days - :param split_train: ratio of train datasets (Default value = 0.7) - :param split_valid: ratio of validation datasets (Default value = 0.2) - :param split_test: ratio of test datasets (Default value = 0.1) - - """ - - if split_train + split_valid + split_test != 1: - ValueError("summed Split ratios not equal 1! Please adjust the values") - damping_days = np.asarray(damping_days) - n = damping_days.shape[0] - n_train = int(n * split_train) - n_valid = int(n * split_valid) - n_test = n - n_train - n_valid - - damping_days_train, damping_days_valid, damping_days_test = tf.split( - damping_days, [n_train, n_valid, n_test], 0) - data = { - "train": tf.reshape(damping_days_train, [n_train, 1]), - "valid": tf.reshape(damping_days_valid, [n_valid, 1]), - "test": tf.reshape(damping_days_test, [n_test, 1]) - } - - return data - - def get_input_dim_lstm(path): - """ Extract the dimensiond of the input data + """ Extract the dimension of the input data :param path: path to the data @@ -489,29 +493,31 @@ def get_input_dim_lstm(path): path = os.path.dirname(os.path.realpath(__file__)) path_data = os.path.join(os.path.dirname(os.path.realpath( os.path.dirname(os.path.realpath(path)))), 'data') - max_epochs = 100 - label_width = 30 - - input_dim = get_input_dim_lstm(path_data) - model = "CNN" - if model == "Dense_Single": - model = network_architectures.mlp_multi_input_single_output() - modeltype = 'classic' - - elif model == "Dense": - model = network_architectures.mlp_multi_input_multi_output(label_width) - modeltype = 'classic' - - elif model == "LSTM": - model = network_architectures.lstm_multi_input_multi_output( - label_width) - modeltype = 'timeseries' - - elif model == "CNN": - model = network_architectures.cnn_multi_input_multi_output(label_width) - modeltype = 'timeseries' + # Defining model parameters + label_width = 30 + number_age_groups = 6 + number_compartments = 8 + hidden_layers = 3 + neurons_in_hidden_layer = 512 + activation_function = 'relu' + modelname = "Dense" + modeltype = "classic" # or "classic" + + model_parameters = (label_width, number_age_groups, number_compartments, + hidden_layers, neurons_in_hidden_layer, activation_function, modelname) + + # Defining training parameters + early_stop = 100 + max_epochs = 200 + loss = tf.keras.losses.MeanAbsolutePercentageError() + optimizer = "AdamW" + metrics = [tf.keras.metrics.MeanAbsoluteError()] + training_parameters = (early_stop, max_epochs, loss, optimizer, metrics) + + # input_dim = get_input_dim_lstm(path_data) -> Warum? + model = initialize_model(model_parameters) model_output = network_fit( - path_data, model=model, modeltype=modeltype, - max_epochs=max_epochs) + model=model, modeltype=modeltype, + training_parameter=training_parameters, path=path_data) diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/network_architectures.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/network_architectures.py index 5c623be085..5b2eab266c 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/network_architectures.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_groups/network_architectures.py @@ -20,87 +20,203 @@ import tensorflow as tf -def mlp_multi_input_single_output(num_age_groups=6): - """ Simple MLP Network which takes the compartments for multiple time steps as input and - returns the 8 compartments for all six age groups for one single time step. +def mlp_multi_input_single_output(num_age_groups=6, num_outputs=8, num_hidden_layers=3, num_neurons_per_layer=32, activation='relu'): + """ Simple MLP Network which takes the compartments for multiple time steps as input and returns the compartments for one single time step for all age groups. - Reshaping adds an extra dimension to the output, so the shape of the output is 1x48. - This makes the shape comparable to that of the multi-output models. - - :param num_age_groups: Number of age groups in population. (Default value = 6) + Reshaping adds an extra dimension to the output, so the shape of the output is 1x8. This makes the shape comparable to that of the multi-output models. + :param num_age_groups: Number of the age groups, default 6 + :param num_outputs: Default: 8 Number of compartments. Default value is reached when aggregating the confirmed compartments. + :param num_hidden_layers: Number of hidden dense layers in the MLP architecture + :param num_neurons_per_layer: Number of neurons per hidden layer + :param activation: name of the used activation function + :returns: Tensorflow keras model with given architecture """ - model = tf.keras.Sequential([ - tf.keras.layers.Flatten(), - tf.keras.layers.Dense(units=32, activation='relu'), - tf.keras.layers.Dense(units=32, activation='relu'), - tf.keras.layers.Dense(units=8*num_age_groups), - tf.keras.layers.Reshape([1, -1]), ]) + + # Catching unallowed inputs + if num_age_groups < 1: + raise ValueError( + "Number of age groups has to be positive, here %d" % ( + num_age_groups) + ) + if num_outputs < 1: + raise ValueError( + "Output dimension must be at least 1, here %d" % (num_outputs)) + if num_hidden_layers < 0: + raise ValueError( + "Number of layers must be at least 0, here %d" % (num_hidden_layers)) + if num_neurons_per_layer < 1: + raise ValueError("Number of neurons per layer must be at least 1, here %d" % ( + num_neurons_per_layer)) + + # Setting up basic model + model = tf.keras.Sequential() + model.add(tf.keras.layers.Flatten()) + + # Adding new hidden layers + for _ in range(num_hidden_layers): + model.add(tf.keras.layers.Dense( + units=num_neurons_per_layer, activation=activation)) + + # Adding output layer and reshaping output + model.add(tf.keras.layers.Flatten()) + model.add(tf.keras.layers.Dense(units=num_outputs*num_age_groups)) + return model -def mlp_multi_input_multi_output(label_width, num_age_groups=6): - """ Simple MLP Network which takes the compartments for multiple time steps as input and - returns the 8 compartments for all age groups for multiple time steps in the future. +def mlp_multi_input_multi_output(label_width, num_age_groups=6, num_outputs=8, num_hidden_layers=3, num_neurons_per_layer=32, activation='relu'): + """ Simple MLP Network which takes the compartments for multiple time steps as input and returns the compartments for multiple time step. - Reshaping adds an extra dimension to the output, so the shape of the output is (label_width)x48. - This makes the shape comparable to that of the multi-output models. + Reshaping adds an extra dimension to the output, so the shape of the output is label_widthx(num_age_groups*num_outputs). This makes the shape comparable to that of the multi-output models. :param label_width: Number of time steps in the output. - :param num_age_groups: Number of age groups in population. (Default value = 6) - + :param num_age_groups: Number of age groups in the population, default 6. + :param num_outputs: Default: 8 Number of compartments. Default value is reached when aggregating the confirmed compartments. + :param num_hidden_layers: Number of hidden dense layers in the MLP architecture + :param num_neurons_per_layer: Number of neurons per hidden layer + :param activation: name of the used activation function + :returns: MLP network with given architecture """ - model = tf.keras.Sequential([ - tf.keras.layers.Flatten(), - tf.keras.layers.Dense(units=32, activation='relu'), - tf.keras.layers.Dense(units=32, activation='relu'), - tf.keras.layers.Dense(units=label_width*num_age_groups*8), - tf.keras.layers.Reshape([label_width, num_age_groups*8]) - ]) + + # Catching unallowed inputs + if label_width < 1: + raise ValueError("label width has to be a positive integer") + if num_age_groups < 1: + raise ValueError( + "Number of age groups must be positive, here %d" % (num_age_groups)) + if num_outputs < 1: + raise ValueError( + "Output dimension must be at least 1, here %d" % (num_outputs)) + if num_hidden_layers < 0: + raise ValueError( + "Number of layers must be at least 0, here %d" % (num_hidden_layers)) + if num_neurons_per_layer < 1: + raise ValueError("Number of neurons per layer must be at least 1, here %d" % ( + num_neurons_per_layer)) + + # Setting up basic model + model = tf.keras.Sequential() + model.add(tf.keras.layers.Flatten()) + + # Adding new hidden layers + for _ in range(num_hidden_layers): + model.add(tf.keras.layers.Dense( + units=num_neurons_per_layer, activation=activation)) + + # Adding output layer and reshaping output + model.add(tf.keras.layers.Dense(label_width*num_outputs*num_age_groups)) + model.add(tf.keras.layers.Reshape( + [label_width, num_outputs*num_age_groups])) + return model -def cnn_multi_input_multi_output(label_width, num_age_groups=6): - """ CNN Network which uses multiple time steps as input and returns the 8 compartments for - each age group for multiple time steps in the future. +def cnn_multi_input_multi_output(label_width, num_age_groups=6, num_outputs=8, conv_size=3, num_filters=256, num_hidden_layers=1, num_neurons_per_layer=256, activation="relu"): + """ CNN Network which uses multiple time steps as input and returns the compartments for multiple age groups and time steps in the future. - Input and output have shape [number of expert model simulations, time points in simulation, - number of individuals in infection states]. + The parameter conv_size describes the kernel_size of the 1d Conv layer. + We also use the parameter in combination with a lambda layer to transform the input to shape [batch, CONV_WIDTH, features]. :param label_width: Number of time steps in the output. - :param num_age_groups: Number of age groups in population. (Default value = 6) - + :param num_age_groups: Number of age groups in the population. + :param num_outputs: Default: 8 Number of compartments. Default value is reached when aggregating the confirmed compartments. + :param conv_size: Default: 3 Convolution kernel width which is 3 per default. + :param num_filters: Number of different filters used in the Conv1D-Layer + :param num_hidden_layers: Number of layers in the dense network following the convolution layer + :param num_neurons_per_layer: Number of neurons in each of the hidden layers (except the output layer) + :param activation: activation function used in the hidden MLP-layers. + :returns: Tensorflow keras CNN model with given architecture """ - + # Catching unallowed inputs + if label_width < 1: + raise ValueError("label width has to be a positive integer") + if num_age_groups < 1: + raise ValueError( + "Number of age groups has to be positive, here %d" % (num_age_groups)) + if conv_size < 2: + raise ValueError( + "Size of the convolution kernel has to be larger than 1, here %d" % (conv_size)) + if num_outputs < 1: + raise ValueError( + "Output dimension must be at least 1, here %d" % (num_outputs)) + if num_filters < 1: + raise ValueError( + "Number of filters must be at least 1, here %d" % (num_filters) + ) + if num_hidden_layers < 0: + raise ValueError( + "Number of hidden layers must be at least 0, here %d" % (num_hidden_layers)) + if num_neurons_per_layer < 1: + raise ValueError("Number of neurons per layer must be at least 1, here %d" % ( + num_neurons_per_layer)) + + # Defining convolutional layer model = tf.keras.Sequential([ - tf.keras.layers.Conv1D(filters=64, kernel_size=3, activation='relu'), - tf.keras.layers.Conv1D(filters=64, kernel_size=3, activation='relu'), - # tf.keras.layers.MaxPooling1D(pool_size=2), - tf.keras.layers.Flatten(), - tf.keras.layers.GaussianNoise(0.35), - tf.keras.layers.Dense(512, activation='relu'), - # tf.keras.layers.Dense(512, activation='relu'), - tf.keras.layers.Dense( - label_width * num_age_groups * 8, activation='linear'), - tf.keras.layers.Reshape([label_width, 8 * num_age_groups]) + tf.keras.layers.Lambda(lambda x: x[:, -conv_size:, :]), + tf.keras.layers.Conv1D(num_filters, activation='relu', + kernel_size=(conv_size)) ]) - return model + # Constructing the hidden layers + for _ in range(num_hidden_layers): + model.add(tf.keras.layers.Dense( + units=num_neurons_per_layer, activation=activation)) -def lstm_multi_input_multi_output(label_width, num_age_groups=6): - """ LSTM Network which uses multiple time steps as input and returns the 8 compartments for - multiple time steps in the future. + # Adding the output layer + model.add(tf.keras.layers.Dense(label_width*num_outputs*num_age_groups, + kernel_initializer=tf.initializers.zeros())) + # Adding reshaping layer + model.add(tf.keras.layers.Reshape( + [label_width, num_outputs*num_age_groups])) - Input and output have shape [number of expert model simulations, time points in simulation, - number of individuals in infection states]. + return model - :param label_width: Number of time steps in the output. - :param num_age_groups: (Default value = 6) +def lstm_multi_input_multi_output(label_width, num_age_groups=6, num_outputs=8, internal_dimension=32, num_hidden_layers=1, num_neurons_per_layer=32, activation="relu"): + """ LSTM Network which uses multiple time steps as input and returns the compartments for multiple age groups and time steps in the future. + + :param label_width: Number of time steps in the output. + :param num_age_groups: Number of age groups in the population, default 6. + :param num_outputs: Default: 8 Number of compartments. Default value is reached when aggregating the confirmed compartments. + :param internal_dimension: Output dimension of the LSTM-layer. + :param num_hidden_layers: Number of hidden layers in the dense network + :param num_neurons_per_layer: Number of neurons per hidden layer + :param activation: Name of the used activation function + :retruns: Tensorflow keras LSTM model with given architecture """ + # Catching unallowed inputs + if label_width < 1: + raise ValueError("label width has to be a positive integer") + if num_age_groups < 1: + raise ValueError( + "Number of age groups has to be positive, here %d" % (num_age_groups)) + if num_outputs < 1: + raise ValueError( + "Output dimension must be at least 1, here %d" % (num_outputs)) + if internal_dimension < 1: + raise ValueError( + "Internal dimension must be at least 1, here %d" % (internal_dimension)) + if num_hidden_layers < 0: + raise ValueError( + "Number of hidden layers must be at least 0, here %d" % (num_hidden_layers)) + if num_neurons_per_layer < 1: + raise ValueError("Number of neurons per layer must be at least 1, here %d" % ( + num_neurons_per_layer)) + + # Defining the model with one LSTM layer model = tf.keras.Sequential([ - tf.keras.layers.LSTM(32, return_sequences=False), - tf.keras.layers.Dense(label_width * 8 * num_age_groups, - kernel_initializer=tf.initializers.zeros()), - tf.keras.layers.Reshape([label_width, 8 * num_age_groups])]) + tf.keras.layers.LSTM(internal_dimension, return_sequences=False) + ]) + + # Constructing the hidden layers + for _ in range(num_hidden_layers): + model.add(tf.keras.layers.Dense( + units=num_neurons_per_layer, activation=activation)) + + # Adding output and reshape layer + model.add(tf.keras.layers.Dense(label_width*num_outputs*num_age_groups, + kernel_initializer=tf.initializers.zeros())) + model.add(tf.keras.layers.Reshape( + [label_width, num_outputs*num_age_groups])) return model diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/data_generation.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/data_generation.py index f10d7709fe..f7541f5e9f 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/data_generation.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/data_generation.py @@ -29,31 +29,42 @@ from progress.bar import Bar from sklearn.preprocessing import FunctionTransformer -from memilio.simulation import (AgeGroup, ContactMatrix, Damping, LogLevel, - UncertainContactMatrix, set_log_level) +from memilio.simulation import (AgeGroup, LogLevel, set_log_level) from memilio.simulation.osecir import (Index_InfectionState, - InfectionState, Model, Simulation, + InfectionState, Model, interpolate_simulation_result, simulate) def remove_confirmed_compartments(result_array): + """Aggregates the confirmed and non-confirmed infection compartments. + + :param result_array: Array containing simulation results with compartment populations + :returns: Modified array with aggregated compartments """ - :param result_array: + # Defining relevant indices + indices_inf_no_symp = [2, 3] + index_inf_no_symp = 2 + indices_inf_symp = [4, 5] + index_inf_symp = 4 + indices_removed = [3, 5] - """ - sum_inf_no_symp = np.sum(result_array[:, [2, 3]], axis=1) - sum_inf_symp = np.sum(result_array[:, [2, 3]], axis=1) - result_array[:, 2] = sum_inf_no_symp - result_array[:, 4] = sum_inf_symp - return np.delete(result_array, [3, 5], axis=1) + # Calculating the values for the merged compartments + sum_inf_no_symp = np.sum(result_array[:, indices_inf_no_symp], axis=1) + sum_inf_symp = np.sum(result_array[:, indices_inf_symp], axis=1) + + # Seting the new values + result_array[:, index_inf_no_symp] = sum_inf_no_symp + result_array[:, index_inf_symp] = sum_inf_symp + + return np.delete(result_array, indices_removed, axis=1) def run_secir_simple_simulation(days): """ Uses an ODE SECIR model allowing for asymptomatic infection. The model is not stratified by region or demographic properties such as age. Virus-specific parameters are fixed and initial number of persons in the particular infection states are chosen randomly from defined ranges. - :param days: Describes how many days we simulate within a single run. + :param days: Describes how many days we simulate within a single run. :returns: List containing the populations in each compartment for each day of the simulation. """ @@ -129,8 +140,6 @@ def run_secir_simple_simulation(days): result_array = remove_confirmed_compartments( result_array[1:, :].transpose()) - dataset = [] - dataset_entries = copy.deepcopy(result_array) return dataset_entries.tolist() @@ -199,7 +208,13 @@ def generate_data( os.mkdir(path) # save dict to json file - with open(os.path.join(path, 'data_secir_simple.pickle'), 'wb') as f: + if num_runs > 1000: + filename = "data_secir_simple_%ddays_%dk.pickle" % ( + label_width, num_runs//1000) + else: + filename = "data_secir_simple_%ddays_%d.pickle" % ( + label_width, num_runs) + with open(os.path.join(path, filename), 'wb') as f: pickle.dump(data, f) return data @@ -212,6 +227,6 @@ def generate_data( input_width = 5 label_width = 30 - num_runs = 1000 + num_runs = 10000 data = generate_data(num_runs, path_data, input_width, label_width) diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/grid_search.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/grid_search.py new file mode 100644 index 0000000000..869c60404d --- /dev/null +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/grid_search.py @@ -0,0 +1,174 @@ +############################################################################# +# Copyright (C) 2020-2025 MEmilio +# +# Authors: Manuel Heger, Henrik Zunker +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# + +import os +import tensorflow as tf +import pickle + +import pandas as pd +import time +from sklearn.model_selection import KFold +import numpy as np +import memilio.surrogatemodel.ode_secir_simple.network_architectures as architectures +import memilio.surrogatemodel.ode_secir_simple.model as md + +# Function to train and evaluate the model using cross-validation + + +def train_and_evaluate_model(param, inputs, labels, training_parameter, show_results=False): + """ Training and evaluating a model with given architecture using 5-fold cross validation, returning a dictionary with the main training statistics. + + :param param: tuple of parameters describing the model architecture, it should be of the form + (num_days_per_output, num_outputs, num_hidden_layers, neurons_per_layer, name_activation, name_architecture) + :param inputs: training inputs + :param labels: training output labels + :param training_parameter: tuple of parameters used for the training process, it should be of the form + (early_stop, max_epochs, loss, optimizer, metrics), where loss is a loss-function implemented in keras, optimizer is the name of the used optimizer, + metrics is a list of used training metrics, e.g. [tf.keras.metrics.MeanAbsoluteError(), tf.keras.metrics.MeanAbsolutePercentageError()] + :param show_results: Boolean, whether or not the evaluation results are printed. + :returns: a dictionary of training statistics of the form + {"model", "activation","optimizer","mean_train_loss_kfold","mean_val_loss_kfold","training_time", "train_losses", "val_losses"} + + """ + # Unpacking parameters + _, _, _, _, activation, modelname = param + early_stop, max_epochs, loss, optimizer, metrics = training_parameter + early_stopping = tf.keras.callbacks.EarlyStopping(monitor='val_loss', + patience=early_stop, + mode='min') + + # Preparing K-Fold Cross-Validation + kf = KFold(n_splits=5) + train_losses = [] + val_losses = [] + + losses_history_all = [] + val_losses_history_all = [] + + start = time.perf_counter() + for train_idx, val_idx in kf.split(inputs): + # Clearing any information about the previous model + tf.keras.backend.clear_session() + + # Gather training and validation data based on the fold + train_inputs = tf.gather(inputs, indices=train_idx) + train_labels = tf.gather(labels, indices=train_idx) + valid_inputs = tf.gather(inputs, indices=val_idx) + valid_labels = tf.gather(labels, indices=val_idx) + + # Initializing model + model = md.initialize_model(param) + # Compile the model + model.compile(loss=loss, + optimizer=optimizer, + metrics=metrics) + + # Train the model + history = model.fit(train_inputs, train_labels, epochs=max_epochs, + validation_data=(valid_inputs, valid_labels), + callbacks=[early_stopping]) + + train_losses.append(np.min(history.history['loss'])) + val_losses.append(np.min(history.history['val_loss'])) + losses_history_all.append(history.history['loss']) + val_losses_history_all.append(history.history['val_loss']) + + elapsed = time.perf_counter() - start + + # Print out the results + if show_results: + print(f"Best train losses: {train_losses}") + print(f"Best validation losses: {val_losses}") + print("--------------------------------------------") + print(f"K-Fold Train Score: {np.mean(train_losses)}") + print(f"K-Fold Validation Score: {np.mean(val_losses)}") + print(f"Time for training: {elapsed:.4f} seconds") + print(f"Time for training: {elapsed / 60:.4f} minutes") + + # After cross-validation, we can test on the withhold dataset (outside of the loop) ? + return { + "model": modelname, + "activation": activation, + "optimizer": optimizer, + "mean_train_loss_kfold": np.mean(train_losses), + "mean_val_loss_kfold": np.mean(val_losses), + "training_time": elapsed/60, + "train_losses": [losses_history_all], + "val_losses": [val_losses_history_all] + } + + +def perform_grid_search(model_parameters, inputs, labels, training_parameters, filename_df, path=None): + """ Performing grid search for a given set of model parameters + + The results are stored in directory 'secir_simple_grid_search', each row has the form + ['model', 'optimizer', 'number_of_hidden_layers', 'number_of_neurons', 'activation', + 'mean_test_MAPE', 'kfold_train', 'kfold_val', + 'kfold_test', 'training_time', 'train_losses', 'val_losses'] + + :param model_parameters: List of tuples of model parameters, each entry should be a tuple of the form + (num_days_per_output, num_outputs, num_hidden_layers, neurons_per_layer, name_activation, name_architecture) + :param inputs: training input data + :param labels: training label data + :param training_parameters: List of tuples of parameters used for the training process, each should be of the form + (early_stop, max_epochs, loss, optimizer, metrics), where loss is a loss-function implemented in keras, optimizer is the name of the used optimizer, + metrics is a list of used training metrics, e.g. [tf.keras.metrics.MeanAbsoluteError(), tf.keras.metrics.MeanAbsolutePercentageError()] + :param filename_df: String, giving name of the file, where the data is stored, actual filename is given by filename_df + ".pickle" + :param path: String representing the path, where dataframe should be stored + """ + # Create a DataFrame to store the results + df_results = pd.DataFrame(columns=['model', 'optimizer', 'number_of_hidden_layers', 'number_of_neurons', 'activation', + 'mean_test_MAPE', 'kfold_train', 'kfold_val', + 'kfold_test', 'training_time', 'train_losses', 'val_losses']) + + # Iterate the different model architectures and save the training results + for param in model_parameters: + for training_parameter in training_parameters: + _, _, layer, neuron_number, activation, modelname = param + results = train_and_evaluate_model( + param, inputs, labels, training_parameter) + df_results.loc[len(df_results.index)] = [ + # Placeholder for test score + modelname, results["optimizer"], layer, neuron_number, activation, np.nan, + results["mean_train_loss_kfold"], + results["mean_val_loss_kfold"], + np.nan, # Placeholder for final test score + results["training_time"], + results["train_losses"], + results["val_losses"] + ] + + # Save the results in file + folder_name = 'secir_simple_grid_search' + + if path is None: + path = os.path.dirname(os.path.realpath(__file__)) + file_path = os.path.join(os.path.dirname(os.path.realpath(path)), + folder_name) + else: + file_path = os.path.join(path, folder_name) + + if not os.path.isdir(file_path): + os.mkdir(file_path) + + file_path = os.path.join(file_path, filename_df + ".pickle") + + with open(file_path, "wb") as f: + pickle.dump(df_results, f) diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/hyperparameter_tuning.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/hyperparameter_tuning.py new file mode 100644 index 0000000000..0123c354c3 --- /dev/null +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/hyperparameter_tuning.py @@ -0,0 +1,98 @@ +############################################################################# +# Copyright (C) 2020-2025 MEmilio +# +# Authors: Manuel Heger, Henrik Zunker +# +# Contact: Martin J. Kuehn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################# + +import os +import tensorflow as tf +import pickle +import pandas as pd +import time +from sklearn.model_selection import KFold +import numpy as np +import memilio.surrogatemodel.ode_secir_simple.grid_search as grid_search +import memilio.surrogatemodel.ode_secir_simple.model as md + +# Setting random seed +np.random.seed(42) + +label_width = 30 +num_outputs = 8 + +# Name of the file containing training data +filename = "data_secir_simple_30days_10k.pickle" + +# Saving grid search results +filename_df = "dataframe_optimizer_paper.csv" + +# General grid search parameters for the training process: +early_stop = [100] +max_epochs = [200] +losses = [tf.keras.losses.MeanAbsolutePercentageError()] +optimizers = ['Adam', 'AdamW', 'Nadam', 'SGD', 'Adagrad', 'RMSProp'] +metrics = [[tf.keras.metrics.MeanAbsoluteError( +), tf.keras.metrics.MeanAbsolutePercentageError()]] + +# Define grid search parameters for the architecture +hidden_layers = [1, 2] +neurons_in_hidden_layer = [32] +activation_function = ['relu', 'elu', 'softmax', 'sigmoid', 'linear', 'tanh'] +models = ["LSTM"] + +# Collecting parameters +training_parameters = [(early, epochs, loss, optimizer, metric) + for early in early_stop for epochs in max_epochs for loss in losses + for optimizer in optimizers for metric in metrics] + +model_parameters = [(label_width, num_outputs, layer, neuron_number, activation, modelname) + for layer in hidden_layers for neuron_number in neurons_in_hidden_layer + for activation in activation_function for modelname in models] + +# Loading data +path = os.path.dirname(os.path.realpath(__file__)) +path_data = os.path.join(os.path.dirname(os.path.realpath( + os.path.dirname(os.path.realpath(path)))), 'data') + +if not os.path.isfile(os.path.join(path_data, filename)): + raise ValueError(f"No dataset found in path: {path_data}") + +with open(os.path.join(path_data, filename), 'rb') as file: + data = pickle.load(file) + +# Split the data: 80% will be used for grid search with cross-validation, and the remaining 20% is withheld for testing +data_splitted = md.split_data(data['inputs'], data['labels'], 0.8, 0, 0.2) +inputs_grid_search = data_splitted['train_inputs'] +labels_grid_search = data_splitted['train_labels'] +inputs_withhold = data_splitted['test_inputs'] +labels_withhold = data_splitted['test_labels'] + + +start_hyper = time.perf_counter() + +# Performing grid search +grid_search.perform_grid_search( + model_parameters, inputs_grid_search, labels_grid_search, training_parameters, filename_df) + +elapsed_hyper = time.perf_counter() - start_hyper + +print( + "Time for hyperparameter testing: {:.4f} minutes".format( + elapsed_hyper / 60)) +print( + "Time for hyperparameter testing: {:.4f} hours".format( + elapsed_hyper / 60 / 60)) diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/model.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/model.py index 2aa5556ba0..4033a64120 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/model.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/model.py @@ -26,7 +26,133 @@ import tensorflow as tf from memilio.simulation.osecir import InfectionState -from memilio.surrogatemodel.ode_secir_simple import network_architectures +import memilio.surrogatemodel.ode_secir_simple.network_architectures as architectures + + +def initialize_model(parameters): + """ Initialize model from given list of parameters + + :param parameters: tuple of parameters describing the model architecture, it should be of the form + (num_days_per_output, num_outputs, num_hidden_layers, neurons_per_layer, name_activation, name_architecture) + :returns: tensor flow keras model with the given architecture + """ + label_width, num_outputs, layer, neuron_number, activation, modelname = parameters + + if modelname == "Dense": + return architectures.mlp_multi_input_multi_output( + label_width=label_width, + num_outputs=num_outputs, + num_hidden_layers=layer, + num_neurons_per_layer=neuron_number, + activation=activation + ) + elif modelname == "LSTM": + return architectures.lstm_multi_input_multi_output( + label_width=label_width, + num_outputs=num_outputs, + num_hidden_layers=layer, + num_neurons_per_layer=neuron_number, + activation=activation + ) + elif modelname == "CNN": + return architectures.cnn_multi_input_multi_output( + label_width=label_width, + num_outputs=num_outputs, + num_hidden_layers=layer, + num_neurons_per_layer=neuron_number, + activation=activation + ) + else: + raise ValueError( + "name_architecture must be one of 'Dense', 'LSTM' or 'CNN'") + + +def network_fit(model, inputs, labels, training_parameter, plot=True): + """ Training and evaluation of a given model under given training parameters. + + :param model: Keras sequential model. + :param inputs: Input data used for training + :param labels: Output data (labels) used for training + :param training_parameter: tuple of parameters used for the training process, it should be of the form + (early_stop, max_epochs, loss, optimizer, metrics), where loss is a loss-function implemented in keras, optimizer is the name of the used optimizer, + metrics is a list of used training metrics, e.g. [tf.keras.metrics.MeanAbsoluteError(), tf.keras.metrics.MeanAbsolutePercentageError()] + :param plot: Default: True] + :returns: trained model with training history (output of model.fit()) + """ + early_stop, max_epochs, loss, optimizer, metrics = training_parameter + data_splitted = split_data(inputs, labels) + + train_inputs = data_splitted['train_inputs'] + train_labels = data_splitted['train_labels'] + valid_inputs = data_splitted['valid_inputs'] + valid_labels = data_splitted['valid_labels'] + test_inputs = data_splitted['test_inputs'] + test_labels = data_splitted['test_labels'] + + early_stopping = tf.keras.callbacks.EarlyStopping(monitor='val_loss', + patience=early_stop, + mode='min') + + model.compile( + loss=loss, + optimizer=optimizer, + metrics=metrics) + + history = model.fit(train_inputs, train_labels, epochs=max_epochs, + validation_data=(valid_inputs, valid_labels), + callbacks=[early_stopping]) + + if (plot): + plot_losses(history) + plot_compartment_prediction_model( + test_inputs, test_labels, model=model, + plot_compartment='InfectedSymptoms', max_subplots=3) + df = get_test_statistic(test_inputs, test_labels, model) + print(df) + + return history + + +def split_data(inputs, labels, split_train=0.7, + split_valid=0.2, split_test=0.1): + """ Split data set in training, validation and testing data sets. + + :param inputs: input dataset + :param labels: label dataset + :param split_train: Share of training data sets. (Default value = 0.7) + :param split_valid: Share of validation data sets. (Default value = 0.2) + :param split_test: Share of testing data sets. (Default value = 0.1) + :returns: dictionary of the form {'train_inputs', 'train_labels', 'valid_inputs', + 'valid_labels', 'test_inputs', 'test_labels'} + """ + + if split_train + split_valid + split_test > 1 + 1e-10: + raise ValueError( + "Summed data set shares are greater than 1. Please adjust the values.") + elif inputs.shape[0] != labels.shape[0] or inputs.shape[2] != labels.shape[2]: + raise ValueError( + "Number of batches or features different for input and labels") + + n = inputs.shape[0] + n_train = int(n * split_train) + n_valid = int(n * split_valid) + n_test = n - n_train - n_valid + + inputs_train, inputs_valid, inputs_test = tf.split( + inputs, [n_train, n_valid, n_test], 0) + labels_train, labels_valid, labels_test = tf.split( + labels, [n_train, n_valid, n_test], 0) + + data = { + 'train_inputs': inputs_train, + 'train_labels': labels_train, + 'valid_inputs': inputs_valid, + 'valid_labels': labels_valid, + 'test_inputs': inputs_test, + 'test_labels': labels_test + } + + return data def plot_compartment_prediction_model( @@ -35,14 +161,15 @@ def plot_compartment_prediction_model( """ Plot prediction of the model and label for one compartment. If model is none, we just plot the inputs and labels for the selected compartment without any predictions. + Plots are saved in directory ./plots/evaluation_secir_simple_{compartment}.png :param inputs: test inputs for model prediction. :param labels: test labels. - :param model: trained model. (Default value = None) - :param plot_col: string name of compartment to be plotted. - :param max_subplots: Number of the simulation runs to be plotted and compared against. (Default value = 8) - :param plot_compartment: (Default value = 'InfectedSymptoms') - + :param model: Default: None] trained model. + :param plot_col: string name of compartment to be plotted, + :param max_subplots: Default: 8] Number of the simulation runs to be plotted and compared against. + :param plot_compartment: possible compartment names "Susceptible", "Exposed", + "InfectedNoSymptoms", "InfectedSevere", "InfectedCritical", "Recovered", "Dead" """ input_width = inputs.shape[1] @@ -88,57 +215,8 @@ def plot_compartment_prediction_model( plt.savefig('plots/evaluation_secir_simple_' + plot_compartment + '.png') -def network_fit(path, model, max_epochs=30, early_stop=500, plot=True): - """ Training and evaluation of a given model with mean squared error loss and Adam optimizer using the mean absolute error as a metric. - - :param path: path of the dataset. - :param model: Keras sequential model. - :param max_epochs: int maximum number of epochs in training. (Default value = 30) - :param early_stop: Integer that forces an early stop of training if the given number of epochs does not give a significant reduction of validation loss. (Default value = 500) - :param plot: (Default value = True) - - """ - - if not os.path.isfile(os.path.join(path, 'data_secir_simple.pickle')): - ValueError("no dataset found in path: " + path) - - file = open(os.path.join(path, 'data_secir_simple.pickle'), 'rb') - - data = pickle.load(file) - data_splitted = split_data(data['inputs'], data['labels']) - - train_inputs = data_splitted['train_inputs'] - train_labels = data_splitted['train_labels'] - valid_inputs = data_splitted['valid_inputs'] - valid_labels = data_splitted['valid_labels'] - test_inputs = data_splitted['test_inputs'] - test_labels = data_splitted['test_labels'] - - early_stopping = tf.keras.callbacks.EarlyStopping(monitor='val_loss', - patience=early_stop, - mode='min') - - model.compile( - loss=tf.keras.losses.MeanSquaredError(), - optimizer=tf.keras.optimizers.Adam(), - metrics=[tf.keras.metrics.MeanAbsoluteError()]) - - history = model.fit(train_inputs, train_labels, epochs=max_epochs, - validation_data=(valid_inputs, valid_labels), - callbacks=[early_stopping]) - - if (plot): - plot_losses(history) - plot_compartment_prediction_model( - test_inputs, test_labels, model=model, - plot_compartment='InfectedSymptoms', max_subplots=3) - df = get_test_statistic(test_inputs, test_labels, model) - print(df) - return history - - def plot_losses(history): - """ Plots the losses of the model training. + """ Plotting the losses of the model training. :param history: model training history. @@ -161,7 +239,7 @@ def get_test_statistic(test_inputs, test_labels, model): :param test_inputs: inputs from test data. :param test_labels: labels (output) from test data. :param model: trained model. - + :returns: Dataframe with the corresponding MAPEs of the different compartments """ pred = model(test_inputs) @@ -173,70 +251,61 @@ def get_test_statistic(test_inputs, test_labels, model): # reshape [batch, time, features] -> [features, time * batch] relative_err_transformed = relative_err.transpose(2, 0, 1).reshape(8, -1) relative_err_means_percentage = relative_err_transformed.mean(axis=1) * 100 + compartments = [str(compartment).split('.')[1] + for compartment in InfectionState.values()] + + # Excluding combined compartments + compartments.pop(5) + compartments.pop(3) mean_percentage = pd.DataFrame( data=relative_err_means_percentage, - index=[str(compartment).split('.')[1] - for compartment in InfectionState.values()], + index=[compartments], columns=['Percentage Error']) return mean_percentage -def split_data(inputs, labels, split_train=0.7, - split_valid=0.2, split_test=0.1): - """ Split data set in training, validation and testing data sets. - - :param inputs: input dataset - :param labels: label dataset - :param split_train: Share of training data sets. (Default value = 0.7) - :param split_valid: Share of validation data sets. (Default value = 0.2) - :param split_test: Share of testing data sets. (Default value = 0.1) - - """ - - if split_train + split_valid + split_test > 1 + 1e-10: - raise ValueError( - "Summed data set shares are greater than 1. Please adjust the values.") - elif inputs.shape[0] != labels.shape[0] or inputs.shape[2] != labels.shape[2]: - raise ValueError( - "Number of batches or features different for input and labels") - - n = inputs.shape[0] - n_train = int(n * split_train) - n_valid = int(n * split_valid) - n_test = n - n_train - n_valid - - inputs_train, inputs_valid, inputs_test = tf.split( - inputs, [n_train, n_valid, n_test], 0) - labels_train, labels_valid, labels_test = tf.split( - labels, [n_train, n_valid, n_test], 0) - - data = { - 'train_inputs': inputs_train, - 'train_labels': labels_train, - 'valid_inputs': inputs_valid, - 'valid_labels': labels_valid, - 'test_inputs': inputs_test, - 'test_labels': labels_test - } - - return data - - if __name__ == "__main__": + label_width = 30 + num_outputs = 8 + # General parameters for the training process: + early_stop = [100] + max_epochs = [200] + losses = [tf.keras.losses.MeanAbsolutePercentageError()] + optimizers = ['AdamW'] + metrics = [[tf.keras.metrics.MeanAbsoluteError(), + tf.keras.metrics.MeanAbsolutePercentageError()]] + + # Define parameters for the architecture + hidden_layers = [3] + neurons_in_hidden_layer = [32] + activation_function = ['relu'] + models = ["LSTM"] + + # Collecting parameters + training_parameters = [(early, epochs, loss, optimizer, metric) + for early in early_stop for epochs in max_epochs for loss in losses + for optimizer in optimizers for metric in metrics] + + model_parameters = [(label_width, num_outputs, layer, neuron_number, activation, modelname) + for layer in hidden_layers for neuron_number in neurons_in_hidden_layer + for activation in activation_function for modelname in models] + + # getting training data + filename = "data_secir_simple_30days_10k.pickle" path = os.path.dirname(os.path.realpath(__file__)) path_data = os.path.join(os.path.dirname(os.path.realpath( os.path.dirname(os.path.realpath(path)))), 'data') - max_epochs = 400 - model = "LSTM" - if model == "Dense": - model = network_architectures.mlp_multi_input_single_output() - elif model == "LSTM": - model = network_architectures.lstm_multi_input_multi_output(30) - elif model == "CNN": - model = network_architectures.cnn_multi_input_multi_output(30) + if not os.path.isfile(os.path.join(path_data, filename)): + raise ValueError(f"No dataset found in path: {path_data}") + + with open(os.path.join(path_data, filename), 'rb') as file: + data = pickle.load(file) + input_data = data['inputs'] + label_data = data['labels'] + model = initialize_model(model_parameters[0]) model_output = network_fit( - path_data, model=model, - max_epochs=max_epochs) + model, input_data, label_data, training_parameters[0]) + plot_compartment_prediction_model(input_data, label_data, model) diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/network_architectures.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/network_architectures.py index a85bb227a6..8d9278052a 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/network_architectures.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel/ode_secir_simple/network_architectures.py @@ -20,39 +20,129 @@ import tensorflow as tf -def mlp_multi_input_single_output(num_outputs=8): +def mlp_multi_input_single_output(num_outputs=8, num_hidden_layers=3, num_neurons_per_layer=32, activation='relu'): """ Simple MLP Network which takes the compartments for multiple time steps as input and returns the 8 compartments for one single time step. Reshaping adds an extra dimension to the output, so the shape of the output is 1x8. This makes the shape comparable to that of the multi-output models. :param num_outputs: Default: 8 Number of compartments. Default value is reached when aggregating the confirmed compartments. + :param num_hidden_layers: Number of hidden dense layers in the MLP architecture + :param num_neurons_per_layer: Number of neurons per hidden layer + :param activation: name of the used activation function + :returns: tensorflow keras model with the desired MLP architecture + """ + + # Catching unallowed inputs + if num_outputs < 1: + raise ValueError( + "Output dimension must be at least 1, here %d" % (num_outputs)) + if num_hidden_layers < 0: + raise ValueError( + "Number of layers must be at least 0, here %d" % (num_hidden_layers)) + if num_neurons_per_layer < 1: + raise ValueError("Number of neurons per layer must be at least 1, here %d" % ( + num_neurons_per_layer)) + + # Setting up basic model + model = tf.keras.Sequential() + model.add(tf.keras.layers.Flatten()) + + # Adding new hidden layers + for _ in range(num_hidden_layers): + model.add(tf.keras.layers.Dense( + units=num_neurons_per_layer, activation=activation)) + + # Adding output layer and reshaping output + model.add(tf.keras.layers.Flatten()) + model.add(tf.keras.layers.Dense(units=num_outputs)) + + return model + + +def mlp_multi_input_multi_output(label_width, num_outputs=8, num_hidden_layers=3, num_neurons_per_layer=32, activation='relu'): + """ Simple MLP Network which takes the compartments for multiple time steps as input and returns the 8 compartments for multiple time steps. + Reshaping adds an extra dimension to the output, so the shape of the output is 1xnum_outputs. This makes the shape comparable to that of the multi-output models. + + :param label_width: Number of time steps in the output. + :param num_outputs: Default: 8 Number of compartments. Default value is reached when aggregating the confirmed compartments. + :param num_hidden_layers: Number of hidden dense layers in the MLP architecture + :param num_neurons_per_layer: Number of neurons per hidden layer + :param activation: name of the used activation function + :returns: tensorflow keras model with the desired MLP architecture """ - model = tf.keras.Sequential([ - tf.keras.layers.Flatten(), - tf.keras.layers.Dense(units=32, activation='relu'), - tf.keras.layers.Dense(units=32, activation='relu'), - tf.keras.layers.Dense(units=num_outputs), - tf.keras.layers.Reshape([1, -1]), ]) + # Catching unallowed inputs + if label_width < 1: + raise ValueError("label width has to be a positive integer") + if num_outputs < 1: + raise ValueError( + "Output dimension must be at least 1, here %d" % (num_outputs)) + if num_hidden_layers < 0: + raise ValueError( + "Number of layers must be at least 0, here %d" % (num_hidden_layers)) + if num_neurons_per_layer < 1: + raise ValueError("Number of neurons per layer must be at least 1, here %d" % ( + num_neurons_per_layer)) + + # Setting up basic model + model = tf.keras.Sequential() + model.add(tf.keras.layers.Flatten()) + + # Adding new hidden layers + for _ in range(num_hidden_layers): + model.add(tf.keras.layers.Dense( + units=num_neurons_per_layer, activation=activation)) + + # Adding output layer and reshaping output + model.add(tf.keras.layers.Dense(label_width*num_outputs)) + model.add(tf.keras.layers.Reshape([label_width, num_outputs])) + return model -def lstm_network_multi_input_single_output(num_outputs=8): +def lstm_network_multi_input_single_output(num_outputs=8, internal_dimension=32, num_hidden_layers=1, num_neurons_per_layer=32, activation="relu"): """ LSTM Network which uses multiple time steps as input and returns the 8 compartments for one single time step in the future. Input and output have shape [number of expert model simulations, time points in simulation, number of individuals in infection states]. :param num_outputs: Default: 8 Number of compartments. Default value is reached when aggregating the confirmed compartments. - + :param internal_dimension: Output dimension of the LSTM-layer. + :param num_hidden_layers: Number of hidden layers in the dense network. + :param num_neurons_per_layer: Number of neurons per hidden layer. + :param activation: Name of the used activation function. + :returns: tensorflow keras model with the desired LSTM architecture """ + # Catching unallowed inputs + if num_outputs < 1: + raise ValueError( + "Output dimension must be at least 1, here %d" % (num_outputs)) + if internal_dimension < 1: + raise ValueError( + "Internal dimension must be at least 1, here %d" % (internal_dimension)) + if num_hidden_layers < 0: + raise ValueError( + "Number of layers must be at least 0, here %d" % (num_hidden_layers)) + if num_neurons_per_layer < 1: + raise ValueError("Number of neurons per layer must be at least 1, here %d" % ( + num_neurons_per_layer)) + + # Defining final model starting with LSTM-layer model = tf.keras.models.Sequential([ - tf.keras.layers.LSTM(32, return_sequences=True), - tf.keras.layers.Dense(units=num_outputs) + tf.keras.layers.LSTM(internal_dimension, return_sequences=True) ]) + + # Adding hidden layers in the dense network + for _ in range(num_hidden_layers): + model.add(tf.keras.layers.Dense( + units=num_neurons_per_layer, activation=activation)) + + # Adding output layer and reshaping + model.add(tf.keras.layers.Flatten()) + model.add(tf.keras.layers.Dense(units=num_outputs)) return model -def cnn_multi_input_multi_output(label_width, conv_size=3, num_outputs=8): +def cnn_multi_input_multi_output(label_width, conv_size=3, num_outputs=8, num_filters=256, num_hidden_layers=1, num_neurons_per_layer=256, activation="relu"): """ CNN Network which uses multiple time steps as input and returns the 8 compartments for multiple time step in the future. Input and output have shape [number of expert model simulations, time points in simulation, number of individuals in infection states]. @@ -60,34 +150,96 @@ def cnn_multi_input_multi_output(label_width, conv_size=3, num_outputs=8): We also use the parameter in combination with a lambda layer to transform the input to shape [batch, CONV_WIDTH, features]. :param label_width: Number of time steps in the output. - :param conv_size: Default: 3 Convolution kernel width which is 3 per default. - :param num_outputs: Default: 8 Number of compartments. Default value is reached when aggregating the confirmed compartments. - + :param conv_size: Default: 3] Convolution kernel width which is 3 per default. + :param num_outputs: Default: 8] Number of compartments. Default value is reached when aggregating the confirmed compartments. + :param num_filters: Number of different filters used in the Conv1D-Layer + :param num_hidden_layers: Number of layers in the dense network following the convolution layer + :param num_neurons_per_layer: Number of neurons in each of the hidden layers (except the output layer) + :param activation: activation function used in the hidden MLP-layers. + :returns: tensorflow keras model with the desired CNN architecture """ + # Catching unallowed inputs + if label_width < 1: + raise ValueError("label width has to be a positive integer") + if conv_size < 2: + raise ValueError( + "Size of the convolution kernel has to be larger than 1, here %d" % (conv_size)) + if num_outputs < 1: + raise ValueError( + "Output dimension must be at least 1, here %d" % (num_outputs)) + if num_filters < 1: + raise ValueError( + "Number of filters must be at least 1, here %d" % (num_filters) + ) + if num_hidden_layers < 0: + raise ValueError( + "Number of hidden layers must be at least 0, here %d" % (num_hidden_layers)) + if num_neurons_per_layer < 1: + raise ValueError("Number of neurons per layer must be at least 1, here %d" % ( + num_neurons_per_layer)) + + # Defining convolutional layer model = tf.keras.Sequential([ tf.keras.layers.Lambda(lambda x: x[:, -conv_size:, :]), - tf.keras.layers.Conv1D(256, activation='relu', - kernel_size=(conv_size)), - tf.keras.layers.Dense(label_width*num_outputs, - kernel_initializer=tf.initializers.zeros()), - tf.keras.layers.Reshape([label_width, num_outputs]) + tf.keras.layers.Conv1D(num_filters, activation='relu', + kernel_size=(conv_size)) ]) + # Constructing the hidden layers + for _ in range(num_hidden_layers): + model.add(tf.keras.layers.Dense( + units=num_neurons_per_layer, activation=activation)) + + # Adding the output layer + model.add(tf.keras.layers.Dense(label_width*num_outputs, + kernel_initializer=tf.initializers.zeros())) + # Adding reshaping layer + model.add(tf.keras.layers.Reshape([label_width, num_outputs])) + return model -def lstm_multi_input_multi_output(label_width, num_outputs=8): - """ LSTM Network which uses multiple time steps as input and returns the 8 compartments for one single time step in the future. +def lstm_multi_input_multi_output(label_width, num_outputs=8, internal_dimension=32, num_hidden_layers=1, num_neurons_per_layer=32, activation="relu"): + """ LSTM Network which uses multiple time steps as input and returns the compartments for multiple time step in the future. Input and output have shape [number of expert model simulations, time points in simulation, number of individuals in infection states]. :param label_width: Number of time steps in the output. :param num_outputs: Default: 8 Number of compartments. Default value is reached when aggregating the confirmed compartments. - + :param internal_dimension: Output dimension of the LSTM-layer. + :param num_hidden_layers: Number of hidden layers in the dense network + :param num_neurons_per_layer: Number of neurons per hidden layer + :param activation: Name of the used activation function + :returns: tensorflow keras model with the desired LSTM architecture """ + # Catching unallowed inputs + if label_width < 1: + raise ValueError("label width has to be a positive integer") + if num_outputs < 1: + raise ValueError( + "Output dimension must be at least 1, here %d" % (num_outputs)) + if internal_dimension < 1: + raise ValueError( + "Internal dimension must be at least 1, here %d" % (internal_dimension)) + if num_hidden_layers < 0: + raise ValueError( + "Number of hidden layers must be at least 0, here %d" % (num_hidden_layers)) + if num_neurons_per_layer < 1: + raise ValueError("Number of neurons per layer must be at least 1, here %d" % ( + num_neurons_per_layer)) + + # Defining the model with one LSTM layer model = tf.keras.Sequential([ - tf.keras.layers.LSTM(32, return_sequences=False), - tf.keras.layers.Dense(label_width*num_outputs, - kernel_initializer=tf.initializers.zeros()), - tf.keras.layers.Reshape([label_width, num_outputs])]) + tf.keras.layers.LSTM(internal_dimension, return_sequences=False) + ]) + + # Constructing the hidden layers + for _ in range(num_hidden_layers): + model.add(tf.keras.layers.Dense( + units=num_neurons_per_layer, activation=activation)) + + # Adding output and reshape layer + model.add(tf.keras.layers.Dense(label_width*num_outputs, + kernel_initializer=tf.initializers.zeros())) + model.add(tf.keras.layers.Reshape([label_width, num_outputs])) return model diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py index cb2f28be09..937e9a4efd 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_groups.py @@ -26,6 +26,7 @@ import unittest import numpy as np +import tensorflow as tf import logging # suppress all autograph warnings from tensorflow @@ -49,8 +50,8 @@ def setUp(self): def test_simulation_run(self, mock_baseline, mock_minimum): """ - :param mock_baseline: - :param mock_minimum: + :param mock_baseline: + :param mock_minimum: """ days_1 = 10 @@ -104,9 +105,9 @@ def test_data_generation_runs( self, mock_population, mock_baseline, mock_minimum): """ - :param mock_population: - :param mock_baseline: - :param mock_minimum: + :param mock_population: + :param mock_baseline: + :param mock_minimum: """ @@ -150,9 +151,9 @@ def test_data_generation_save( self, mock_population, mock_baseline, mock_minimum): """ - :param mock_population: - :param mock_baseline: - :param mock_minimum: + :param mock_population: + :param mock_baseline: + :param mock_minimum: """ @@ -167,45 +168,309 @@ def test_data_generation_save( self.assertEqual(os.listdir(self.path), ['data_secir_groups.pickle']) - def test_structures_networks(self): - """ """ +# Testing network_architectures.py + def test_mlp_multi_single(self): + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_single_output( + 0, 12, 12, 12 + ) + error_message = "Number of age groups has to be positive, here 0" + self.assertEqual(str(error.exception), error_message) - model_mlp_multi_input_single_output = network_architectures.mlp_multi_input_single_output() - self.assertEqual(len(model_mlp_multi_input_single_output.layers), 5) - input_zero = np.zeros((1, 5, 8)) - output_zeros = model_mlp_multi_input_single_output(input_zero) - self.assertEqual(output_zeros.shape[0], 1) - self.assertEqual(output_zeros.shape[1], 1) - self.assertEqual(output_zeros.shape[2], 48) + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_single_output( + 12, 0, 12, 12 + ) + error_message = "Output dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) - label_width = 30 - model_mlp_multi_input_multi_output = network_architectures.mlp_multi_input_multi_output( - label_width) - self.assertEqual(len(model_mlp_multi_input_multi_output.layers), 5) - input_zero = np.zeros((1, 1, 8)) - output_zeros = model_mlp_multi_input_multi_output(input_zero) - self.assertEqual(output_zeros.shape[0], 1) - self.assertEqual(output_zeros.shape[1], 30) - self.assertEqual(output_zeros.shape[2], 48) + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_single_output( + 12, 12, -1, 12 + ) + error_message = "Number of layers must be at least 0, here -1" + self.assertEqual(str(error.exception), error_message) - label_width = 5 - model_cnn = network_architectures.cnn_multi_input_multi_output( - label_width, num_age_groups=1) - self.assertEqual(len(model_cnn.layers), 7) - input_zero = np.zeros((1, label_width, 41)) - output_zeros = model_cnn(input_zero) - self.assertEqual(output_zeros.shape[0], 1) - self.assertEqual(output_zeros.shape[1], 5) - self.assertEqual(output_zeros.shape[2], 8) + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_single_output( + 12, 12, 12, 0 + ) + error_message = "Number of neurons per layer must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) - model_lstm_multi = network_architectures.lstm_multi_input_multi_output( - label_width) - self.assertEqual(len(model_lstm_multi.layers), 3) + model = network_architectures.mlp_multi_input_single_output( + 7, 3, 2, 12 + ) + self.assertEqual(len(model.layers), 5) input_zero = np.zeros((1, 5, 8)) - output_zeros = model_lstm_multi(input_zero) + output_zeros = model(input_zero) self.assertEqual(output_zeros.shape[0], 1) - self.assertEqual(output_zeros.shape[1], 5) - self.assertEqual(output_zeros.shape[2], 48) + self.assertEqual(output_zeros.shape[1], 21) + + def test_mlp_multi_multi(self): + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_multi_output( + 0, 12, 12, 12, 12 + ) + error_message = "label width has to be a positive integer" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_multi_output( + 12, 0, 12, 12, 12 + ) + error_message = "Number of age groups must be positive, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_multi_output( + 12, 12, 0, 12, 12 + ) + error_message = "Output dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_multi_output( + 12, 12, 12, -1, 12 + ) + error_message = "Number of layers must be at least 0, here -1" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_multi_output( + 12, 12, 12, 12, 0 + ) + error_message = "Number of neurons per layer must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + model = network_architectures.mlp_multi_input_multi_output( + 12, 3, 7, 2, 16) + self.assertEqual(len(model.layers), 5) + input_zero = np.zeros((5, 3, 7)) + output_zeros = model(input_zero) + self.assertEqual(output_zeros.shape[0], 5) + self.assertEqual(output_zeros.shape[1], 12) + self.assertEqual(output_zeros.shape[2], 21) + + def test_cnn_multi_multi(self): + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + 0, 1, 1, 2, 1, 1, 1 + ) + error_message = "label width has to be a positive integer" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + 1, 0, 1, 2, 1, 1, 1 + ) + error_message = "Number of age groups has to be positive, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + 1, 1, 0, 2, 1, 1, 1 + ) + error_message = "Output dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + 1, 1, 1, 0, 1, 1, 1 + ) + error_message = "Size of the convolution kernel has to be larger than 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + 1, 1, 1, 2, 0, 1, 1 + ) + error_message = "Number of filters must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + 1, 1, 1, 2, 1, -1, 1 + ) + error_message = "Number of hidden layers must be at least 0, here -1" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + 1, 1, 1, 2, 1, 1, 0 + ) + error_message = "Number of neurons per layer must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + model = network_architectures.cnn_multi_input_multi_output( + 21, 4, 3, 3, 256, 2) + self.assertEqual(len(model.layers), 6) + input_zero = np.zeros((12, 5, 7)) + output_zeros = model(input_zero) + # Number of time series + self.assertEqual(output_zeros.shape[0], 12) + # length of one time series + self.assertEqual(output_zeros.shape[1], 21) + # Dimension of one time step + self.assertEqual(output_zeros.shape[2], 12) + + def test_lstm_multi_multi(self): + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 0, 1, 1, 1, 1, 1 + ) + error_message = "label width has to be a positive integer" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 1, 0, 1, 1, 1, 1 + ) + error_message = "Number of age groups has to be positive, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 1, 1, 0, 1, 1, 1 + ) + error_message = "Output dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 1, 1, 1, 0, 1, 1 + ) + error_message = "Internal dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 1, 1, 1, 1, -1, 1 + ) + error_message = "Number of hidden layers must be at least 0, here -1" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 1, 1, 1, 1, 1, 0 + ) + error_message = "Number of neurons per layer must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + model = network_architectures.lstm_multi_input_multi_output( + 21, 4, 3, 12, 3, 12) + self.assertEqual(len(model.layers), 6) + input_zero = np.zeros((12, 5, 7)) + output_zeros = model(input_zero) + # Number of time series + self.assertEqual(output_zeros.shape[0], 12) + # length of one time series + self.assertEqual(output_zeros.shape[1], 21) + # Dimension of one time step + self.assertEqual(output_zeros.shape[2], 12) + +# Testing model.py + def test_calc_split_index(self): + with self.assertRaises(ValueError) as error: + model.calc_split_index( + 10, 0.9, 0.1, 0.1 + ) + error_message = "Summed data set shares are greater than 1. Please adjust the values." + self.assertEqual(str(error.exception), error_message) + split_index = model.calc_split_index(10, 0.7, 0.1, 0.2) + self.assertEqual(split_index, [7, 1, 2]) + + def test_prepare_data_classic(self): + data = { + "inputs": np.zeros((10, 5, 1)), + "labels": np.zeros((10, 2)), + "contact_matrix": [np.zeros((2, 2)) for _ in np.arange(10)], + "damping_day": [[1] for _ in np.arange(10)] + } + + data_new = model.prepare_data_classic(data) + res = tf.cast([0, 0, 0, 0, 0, 0, 0, 0, 0, 1], tf.float32) + self.assertTrue(all(res == data_new["train_inputs"][0])) + self.assertEqual(data_new["train_inputs"].shape, (7, 10)) + self.assertEqual(data_new["valid_labels"].shape, (2, 2)) + + def test_prod_time_series(self): + obj = [1, 0, 2] + ts = model.prod_time_series(obj, 3, 2) + bl = (ts == tf.reshape( + tf.stack([[1, 1], [0, 0], [2, 2]]), [3, 2, 1])) + self.assertTrue(tf.math.reduce_all(bl)) + + """ + def test_prepare_data_timeseries(self): + data = { + "inputs": np.zeros((10, 5, 1)), + "labels": np.zeros((10, 2)), + "contact_matrix": [np.zeros((2, 2)) for _ in np.arange(10)], + "damping_day": [[1] for _ in np.arange(10)] + } + data_new = model.prepare_data_timeseries(data) + res = tf.cast([[0, 0, 0, 0, 0, 1] for _ in np.arange(5)], tf.float16) + self.assertTrue(tf.math.reduce_all(res == data_new["train_inputs"][0])) + """ + + def test_initialize_model(self): + # Helper function to normalize the .getconfig() output + def normalize_config(config): + config.pop('name', None) + for layer in config["layers"]: + layer["config"].pop("name", None) + return config + + label_width = 8 + number_age_groups = 6 + num_outputs = 30 + hidden_layers = 2 + neurons_in_hidden_layer = 32 + activation_function = "relu" + models = ["LSTM", "Dense", "CNN", "Heinz"] + + # Generating valid models + model_lstm = network_architectures.lstm_multi_input_multi_output( + label_width, number_age_groups, num_outputs, 32, hidden_layers, + neurons_in_hidden_layer, activation_function + ) + model_cnn = network_architectures.cnn_multi_input_multi_output( + label_width=label_width, num_age_groups=number_age_groups, num_outputs=num_outputs, + num_hidden_layers=hidden_layers, num_neurons_per_layer=neurons_in_hidden_layer, activation=activation_function + ) + model_mlp = network_architectures.mlp_multi_input_multi_output( + label_width=label_width, num_age_groups=number_age_groups, num_outputs=num_outputs, + num_hidden_layers=hidden_layers, num_neurons_per_layer=neurons_in_hidden_layer, activation=activation_function + ) + + model_parameters = [(label_width, number_age_groups, num_outputs, hidden_layers, + neurons_in_hidden_layer, activation_function, modelname) for modelname in models] + + for parameter in model_parameters: + _, _, _, _, _, _, modelname = parameter + if modelname == "Heinz": + with self.assertRaises(ValueError) as error: + md = model.initialize_model(parameter) + error_message = "name_architecture must be one of 'Dense', 'LSTM' or 'CNN'" + self.assertEqual(str(error.exception), error_message) + else: + md = model.initialize_model(parameter) + if modelname == "Dense": + self.assertDictEqual( + normalize_config(md.get_config()), + normalize_config(model_mlp.get_config()) + ) + if modelname == "LSTM": + self.assertDictEqual( + normalize_config(md.get_config()), + normalize_config(model_lstm.get_config()) + ) + if modelname == "CNN": + self.assertDictEqual( + normalize_config(md.get_config()), + normalize_config(model_cnn.get_config()) + ) @patch('memilio.surrogatemodel.ode_secir_groups.data_generation.getMinimumMatrix', return_value=0.1 * np.ones((6, 6))) @@ -213,43 +478,37 @@ def test_structures_networks(self): return_value=0.6 * np.ones((6, 6))) @patch('memilio.surrogatemodel.ode_secir_groups.data_generation.get_population', return_value=[[5256.0, 10551, 32368.5, - 43637.833333333336, 22874.066666666666, 8473.6]]) + 43637.833333333336, 22874.066666666666, 8473.6]]) def test_network_fit(self, mock_population, mock_baseline, mock_minimum): """ - :param mock_population: - :param mock_baseline: - :param mock_minimum: + :param mock_population: + :param mock_baseline: + :param mock_minimum: """ - max_epochs = 5 - # models with single output - model_mlp_multi_input_single_output = network_architectures.mlp_multi_input_single_output() - + max_epochs = 1 + early_stop = 100 + loss = tf.keras.losses.MeanAbsolutePercentageError() + optimizer = "Adam" + metric = [tf.keras.metrics.MeanAbsoluteError(), + tf.keras.metrics.MeanAbsolutePercentageError()] + + training_parameter = ( + early_stop, max_epochs, loss, optimizer, metric) + + md = network_architectures.mlp_multi_input_multi_output( + label_width=10 + ) # no existing dataset with self.assertRaises(FileNotFoundError) as error: model.network_fit( - self.path, model=model_mlp_multi_input_single_output, - modeltype='classic', max_epochs=max_epochs) + path=self.path, model=md, + modeltype='classic', training_parameter=training_parameter, filename="data_secir_groups.pickle") error_message = "[Errno 2] No such file or directory in the fake filesystem: '" + \ self.path + "data_secir_groups.pickle'" self.assertEqual(str(error.exception), error_message) - # generate dataset with single output - input_width = 5 - label_width = 1 - num_runs = 5 - data_generation.generate_data(num_runs, self.path, "", input_width, - label_width) - - # test different network_architectures - mlp_output = model.network_fit( - self.path, model=model_mlp_multi_input_single_output, - modeltype='classic', max_epochs=max_epochs, plot=False) - - self.assertEqual( - len(mlp_output.history['val_loss']), max_epochs) - # generate dataset with multiple output input_width = 5 label_width = 10 @@ -259,10 +518,12 @@ def test_network_fit(self, mock_population, mock_baseline, mock_minimum): # models with multiple outputs model_mlp_multi_input_multi_output = model.network_fit( - self.path, model=network_architectures.mlp_multi_input_multi_output( label_width), - modeltype='classic', max_epochs=max_epochs, plot=False) + training_parameter=training_parameter, + path=self.path, + filename="data_secir_groups.pickle", + modeltype='classic', plot=False) self.assertEqual( model_mlp_multi_input_multi_output.model.output_shape[1], label_width) @@ -271,20 +532,26 @@ def test_network_fit(self, mock_population, mock_baseline, mock_minimum): max_epochs) model_lstm_multi_output = model.network_fit( - self.path, model=network_architectures.lstm_multi_input_multi_output( label_width), - modeltype='timeseries', max_epochs=max_epochs, plot=False) + modeltype='timeseries', + training_parameter=training_parameter, + path=self.path, + filename="data_secir_groups.pickle", + plot=False) self.assertEqual( model_lstm_multi_output.model.output_shape[1], label_width) self.assertEqual( len(model_lstm_multi_output.history['val_loss']), max_epochs) cnn_output = model.network_fit( - self.path, model=network_architectures.cnn_multi_input_multi_output( label_width), - modeltype='timeseries', max_epochs=max_epochs, plot=False) + modeltype='timeseries', + training_parameter=training_parameter, + path=self.path, + filename="data_secir_groups.pickle", + plot=False) self.assertEqual( cnn_output.model.output_shape[1], label_width) self.assertEqual( diff --git a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_simple.py b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_simple.py index a2aff3ea29..8bd3dd5a13 100644 --- a/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_simple.py +++ b/pycode/memilio-surrogatemodel/memilio/surrogatemodel_test/test_surrogatemodel_ode_secir_simple.py @@ -1,7 +1,7 @@ ############################################################################# # Copyright (C) 2020-2025 MEmilio # -# Authors: +# Authors: Manuel Heger # # Contact: Martin J. Kuehn # @@ -19,12 +19,17 @@ ############################################################################# from pyfakefs import fake_filesystem_unittest -from memilio.surrogatemodel.ode_secir_simple import (data_generation, model, - network_architectures) +from memilio.surrogatemodel.ode_secir_simple import (data_generation, + network_architectures, grid_search) +from memilio.surrogatemodel.ode_secir_simple import model as md import os +import pickle +import pandas as pd import unittest +from unittest.mock import patch import numpy as np +import tensorflow as tf import logging # suppress all autograph warnings from tensorflow @@ -42,7 +47,7 @@ def setUp(self): self.setUpPyfakefs() def test_simulation_run(self): - """ """ + """ Checking if simulated data for one sample using run_secir_simple_simulation has the desired length""" days_1 = 10 days_2 = 30 @@ -57,7 +62,7 @@ def test_simulation_run(self): self.assertEqual(len(simulation_3), days_3+1) def test_data_generation_runs(self): - """ """ + """ Testing if data_generation produces data with desired dimensions """ input_width_1 = 1 input_width_2 = 3 @@ -89,7 +94,7 @@ def test_data_generation_runs(self): self.assertEqual(len(data_2['labels'][0][0]), 8) def test_data_generation_save(self): - """ """ + """ Checking if saving results as pickle-file is possible """ input_width = 2 label_width = 3 @@ -100,110 +105,405 @@ def test_data_generation_save(self): self.assertEqual(len(os.listdir(self.path)), 1) self.assertEqual(os.listdir(self.path), - ['data_secir_simple.pickle']) + ['data_secir_simple_3days_1.pickle']) + + def test_mlp_multi_single(self): + """" testing initialization of MLP architectures works correctly """ + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_single_output( + 0, 1, 1 + ) + error_message = "Output dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) - def test_structures_networks(self): - """ """ + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_single_output( + 1, -1, 1 + ) + error_message = "Number of layers must be at least 0, here -1" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_single_output( + 1, 1, 0 + ) + error_message = "Number of neurons per layer must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) - model_mlp_multi_input_single_output = network_architectures.mlp_multi_input_single_output() - self.assertEqual(len(model_mlp_multi_input_single_output.layers), 5) + model = network_architectures.mlp_multi_input_single_output(3, 4, 84) + # Layers: Flatten layer + hidden layers + flatten layer + output layer + self.assertEqual(len(model.layers), 7) input_zero = np.zeros((1, 5, 8)) - output_zeros = model_mlp_multi_input_single_output(input_zero) + output_zeros = model(input_zero) + # Output has shape (1, num_outputs) self.assertEqual(output_zeros.shape[0], 1) - self.assertEqual(output_zeros.shape[1], 1) - self.assertEqual(output_zeros.shape[2], 8) + self.assertEqual(output_zeros.shape[1], 3) + + def test_mlp_multi_multi(self): + """" testing initialization of MLP architectures works correctly """ + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_multi_output( + 0, 12, 12, 12 + ) + error_message = "label width has to be a positive integer" + self.assertEqual(str(error.exception), error_message) - model_lstm_single = network_architectures.lstm_network_multi_input_single_output() - self.assertEqual(len(model_lstm_single.layers), 2) - input_zero = np.zeros((1, 1, 8)) - output_zeros = model_lstm_single(input_zero) - self.assertEqual(output_zeros.shape[0], 1) - self.assertEqual(output_zeros.shape[1], 1) - self.assertEqual(output_zeros.shape[2], 8) + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_multi_output( + 12, 0, 12, 12 + ) + error_message = "Output dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) - label_width = 5 - model_cnn = network_architectures.cnn_multi_input_multi_output( - label_width) - self.assertEqual(len(model_cnn.layers), 4) - input_zero = np.zeros((1, 5, 8)) - output_zeros = model_cnn(input_zero) + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_multi_output( + 12, 12, -1, 12 + ) + error_message = "Number of layers must be at least 0, here -1" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.mlp_multi_input_multi_output( + 12, 12, 12, 0 + ) + error_message = "Number of neurons per layer must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + model = network_architectures.mlp_multi_input_multi_output( + 12, 3, 4, 84) + # Layers: flatten layers + hidden layers + output layer + reshape layer + self.assertEqual(len(model.layers), 7) + input_zero = np.zeros((2, 5, 8)) + output_zeros = model(input_zero) + # Output has shape (Number of samples, label width, number of outputs) + self.assertEqual(output_zeros.shape[0], 2) + self.assertEqual(output_zeros.shape[1], 12) + self.assertEqual(output_zeros.shape[2], 3) + + def test_lstm_multi_single(self): + """" testing initialization of LSTM architectures works correctly """ + with self.assertRaises(ValueError) as error: + network_architectures.lstm_network_multi_input_single_output( + 0, 12, 12, 12 + ) + error_message = "Output dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_network_multi_input_single_output( + 12, 0, 12, 12 + ) + error_message = "Internal dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_network_multi_input_single_output( + 12, 12, -1, 12 + ) + error_message = "Number of layers must be at least 0, here -1" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_network_multi_input_single_output( + 12, 12, 12, 0 + ) + error_message = "Number of neurons per layer must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + model = network_architectures.lstm_network_multi_input_single_output( + num_outputs=5, num_hidden_layers=4) + # model: LSTM-layer + hidden layers + output layer + reshape layer, here 1+4+1+1 + self.assertEqual(len(model.layers), 7) + input_zero = np.zeros((1, 5, 3)) + output_zeros = model(input_zero) + # Dimensionality output (1, Number of outputs), here (1,5) self.assertEqual(output_zeros.shape[0], 1) self.assertEqual(output_zeros.shape[1], 5) - self.assertEqual(output_zeros.shape[2], 8) - model_lstm_multi = network_architectures.lstm_multi_input_multi_output( - label_width) - self.assertEqual(len(model_lstm_multi.layers), 3) + def test_lstm_multi_multi(self): + """" testing initialization of LSTM architectures works correctly """ + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 0, 12, 12, 12, 12 + ) + error_message = "label width has to be a positive integer" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 12, 0, 12, 12, 12 + ) + error_message = "Output dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 12, 12, 0, 12, 12 + ) + error_message = "Internal dimension must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 12, 12, 12, -1, 12 + ) + error_message = "Number of hidden layers must be at least 0, here -1" + self.assertEqual(str(error.exception), error_message) + + with self.assertRaises(ValueError) as error: + network_architectures.lstm_multi_input_multi_output( + 12, 12, 12, 12, 0 + ) + error_message = "Number of neurons per layer must be at least 1, here 0" + self.assertEqual(str(error.exception), error_message) + + model = network_architectures.lstm_multi_input_multi_output( + 3, 84, 4, 12) + # model: LSTM-layer + hidden layers + output layer + reshape layer, here 1+12+1+1 + self.assertEqual(len(model.layers), 15) input_zero = np.zeros((1, 5, 8)) - output_zeros = model_lstm_multi(input_zero) + output_zeros = model(input_zero) + # Dimensionality output (Number of samples, label width, Number of outputs per day), here (1,3,84) self.assertEqual(output_zeros.shape[0], 1) - self.assertEqual(output_zeros.shape[1], 5) - self.assertEqual(output_zeros.shape[2], 8) + self.assertEqual(output_zeros.shape[1], 3) + self.assertEqual(output_zeros.shape[2], 84) + + def test_cnn_multi_multi(self): + """" testing initialization of CNN architectures works correctly """ + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + label_width=0, conv_size=3, num_outputs=8, num_filters=12, + num_hidden_layers=12, num_neurons_per_layer=12 + ) + error_message = "label width has to be a positive integer" + self.assertEqual(str(error.exception), error_message) - def test_network_fit(self): - """ """ + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + label_width=2, conv_size=1, num_outputs=8, num_filters=12, + num_hidden_layers=12, num_neurons_per_layer=12 + ) + error_message = "Size of the convolution kernel has to be larger than 1, here 1" + self.assertEqual(str(error.exception), error_message) - max_epochs = 5 + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + label_width=3, conv_size=3, num_outputs=-23, num_filters=12, + num_hidden_layers=12, num_neurons_per_layer=12 + ) + error_message = "Output dimension must be at least 1, here -23" + self.assertEqual(str(error.exception), error_message) - # models with single output - model_mlp_multi_input_single_output = network_architectures.mlp_multi_input_single_output() - model_lstm_multi_input_single_output = network_architectures.lstm_network_multi_input_single_output() + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + label_width=3, conv_size=3, num_outputs=8, num_filters=-2, + num_hidden_layers=12, num_neurons_per_layer=12 + ) + error_message = "Number of filters must be at least 1, here -2" + self.assertEqual(str(error.exception), error_message) - # no existing dataset - with self.assertRaises(FileNotFoundError) as error: - model.network_fit( - self.path, model=model_mlp_multi_input_single_output, - max_epochs=max_epochs) - error_message = "[Errno 2] No such file or directory in the fake filesystem: '" + \ - self.path + "data_secir_simple.pickle'" + with self.assertRaises(ValueError) as error: + network_architectures.cnn_multi_input_multi_output( + label_width=3, conv_size=3, num_outputs=8, num_filters=2, + num_hidden_layers=-23, num_neurons_per_layer=12 + ) + error_message = "Number of hidden layers must be at least 0, here -23" self.assertEqual(str(error.exception), error_message) - # generate dataset with single output - input_width = 5 - label_width = 1 - num_runs = 15 - data_generation.generate_data(num_runs, self.path, input_width, - label_width) + model = network_architectures.cnn_multi_input_multi_output( + label_width=3, conv_size=3, num_outputs=42, num_filters=2, + num_hidden_layers=23, num_neurons_per_layer=12 + ) + # lambda layer + conv layer + hidden layers + output layer + reshape layer -> 1+1+23+1+1 = 27 + self.assertEqual(len(model.layers), 27) + input_zero = np.zeros((1, 5, 8)) + output_zeros = model(input_zero) + # Dimensionality output (Number of samples, label width, Number of outputs per day), here (1,3,42) + self.assertEqual(output_zeros.shape[0], 1) + self.assertEqual(output_zeros.shape[1], 3) + self.assertEqual(output_zeros.shape[2], 42) + + def test_initialize_model(self): + """" Testing wheter initialization via initialize_model works correctly """ + # Deleting any automatically added numbering like ("Dense_01" -> "Dense") + def normalize_config(config): + config.pop('name', None) + for layer in config["layers"]: + layer["config"].pop("name", None) + return config + + label_width = 8 + num_outputs = 32 + hidden_layers = [2] + neurons_in_hidden_layer = [32] + activation_function = ['relu'] + models = ["LSTM", "Dense", "CNN", "MLP"] - # test different network_architectures - mlp_output = model.network_fit( - self.path, model=model_mlp_multi_input_single_output, - max_epochs=max_epochs, plot=False) - self.assertEqual( - len(mlp_output.history['val_loss']), max_epochs) - lstm_single_output = model.network_fit( - self.path, model=model_lstm_multi_input_single_output, - max_epochs=max_epochs, plot=False) - self.assertEqual( - len(lstm_single_output.history['val_loss']), max_epochs) + model_lstm = network_architectures.lstm_multi_input_multi_output( + label_width=label_width, num_outputs=num_outputs, num_hidden_layers=2, + num_neurons_per_layer=32, activation='relu' + ) + model_cnn = network_architectures.cnn_multi_input_multi_output( + label_width=label_width, num_outputs=num_outputs, num_hidden_layers=2, + num_neurons_per_layer=32, activation='relu' + ) + model_mlp = network_architectures.mlp_multi_input_multi_output( + label_width=label_width, num_outputs=num_outputs, num_hidden_layers=2, + num_neurons_per_layer=32, activation='relu' + ) + + model_parameters = [(label_width, num_outputs, layer, neuron_number, activation, modelname) + for layer in hidden_layers for neuron_number in neurons_in_hidden_layer + for activation in activation_function for modelname in models] + for parameter in model_parameters: + _, _, _, _, _, modelname = parameter + if modelname == "MLP": + with self.assertRaises(ValueError) as error: + model = md.initialize_model(parameter) + error_message = "name_architecture must be one of 'Dense', 'LSTM' or 'CNN'" + self.assertEqual(str(error.exception), error_message) + else: + model = md.initialize_model(parameter) + if modelname == "Dense": + self.assertDictEqual( + normalize_config(model.get_config()), + normalize_config(model_mlp.get_config())) + if modelname == "LSTM": + self.assertDictEqual( + normalize_config(model.get_config()), + normalize_config(model_lstm.get_config())) + if modelname == "CNN": + self.assertDictEqual( + normalize_config(model.get_config()), + normalize_config(model_cnn.get_config())) - # generate dataset with multiple output - input_width = 5 + def test_network_fit(self): + """ Test for training of network """ + max_epochs = 1 + + early_stop = 100 + loss = tf.keras.losses.MeanAbsolutePercentageError() + optimizer = 'Adam' + metric = [tf.keras.metrics.MeanAbsoluteError(), + tf.keras.metrics.MeanAbsolutePercentageError()] + + training_parameter = ( + early_stop, max_epochs, loss, optimizer, metric) + + # generate dataset with multiple outputs label_width = 10 - num_runs = 25 - data_generation.generate_data(num_runs, self.path, input_width, - label_width) - # models with multiple outputs + inputs = tf.ones([5, 5, 8]) + labels = tf.ones([5, 10, 8]) + model_cnn = network_architectures.cnn_multi_input_multi_output( - label_width) + label_width=label_width, num_filters=1, num_neurons_per_layer=1, + num_hidden_layers=0) model_lstm = network_architectures.lstm_multi_input_multi_output( - label_width) + label_width=label_width, internal_dimension=1, num_neurons_per_layer=1, + num_hidden_layers=0) - cnn_output = model.network_fit( - self.path, model=model_cnn, max_epochs=max_epochs, plot=False) + cnn_output = md.network_fit( + model=model_cnn, inputs=inputs, labels=labels, + training_parameter=training_parameter, plot=False) self.assertEqual( cnn_output.model.output_shape[1], label_width) self.assertEqual( len(cnn_output.history['val_loss']), max_epochs) - lstm_output = model.network_fit( - self.path, model=model_lstm, max_epochs=max_epochs, plot=False) + lstm_output = md.network_fit( + model=model_lstm, inputs=inputs, labels=labels, + training_parameter=training_parameter, plot=False) self.assertEqual( lstm_output.model.output_shape[1], label_width) self.assertEqual( len(lstm_output.history['val_loss']), max_epochs) + def test_train_and_evaluate_model(self): + """ Testing evaluation procedure """ + inputs = tf.ones([5, 5, 8]) + labels = tf.ones([5, 10, 8]) + max_epochs = 1 + early_stop = 100 + loss = tf.keras.losses.MeanAbsolutePercentageError() + optimizer = 'Adam' + metric = [tf.keras.metrics.MeanAbsoluteError(), + tf.keras.metrics.MeanAbsolutePercentageError()] + model_parameter = ( + 10, 8, 1, 5, "relu", "Dense" + ) + + training_parameter = ( + early_stop, max_epochs, loss, optimizer, metric) + + res = grid_search.train_and_evaluate_model( + model_parameter, inputs, labels, training_parameter, False) + self.assertEqual(len(res["val_losses"][0]), 5) + self.assertEqual(len(res["val_losses"][0][0]), max_epochs) + self.assertEqual(res["optimizer"], "Adam") + + @patch('memilio.surrogatemodel.ode_secir_simple.grid_search.train_and_evaluate_model') + def test_perform_grid_search(self, mock_train_and_evaluate_model): + """ Testing gridsearch along a sequence of possible parameter values""" + mock_train_and_evaluate_model.return_value = { + "model": "._.", + "activation": "relu", + "optimizer": "Eva", + "mean_train_loss_kfold": 42, + "mean_val_loss_kfold": 12, + "training_time": 1, + "train_losses": [3], + "val_losses": [6] + } + filename_df = "dataframe_optimizer" + os.mkdir(self.path) + + # General grid search parameters for the training process: + early_stop = [100] + max_epochs = [1] + losses = [tf.keras.losses.MeanAbsolutePercentageError()] + optimizers = ['Adam', 'Adam'] + metrics = [[tf.keras.metrics.MeanAbsoluteError(), + tf.keras.metrics.MeanAbsolutePercentageError()]] + + # Define grid search parameters for the architecture + label_width = 3 + num_outputs = 2 + hidden_layers = [1, 2] + neurons_in_hidden_layer = [32] + activation_function = ['relu'] + models = ["Dense"] + + # Collecting parameters + training_parameters = [(early, epochs, loss, optimizer, metric) + for early in early_stop for epochs in max_epochs for loss in losses + for optimizer in optimizers for metric in metrics] + + model_parameters = [(label_width, num_outputs, layer, neuron_number, activation, modelname) + for layer in hidden_layers for neuron_number in neurons_in_hidden_layer + for activation in activation_function for modelname in models] + + # generate dataset with multiple output + inputs = tf.ones([5, 1, 2]) + labels = tf.ones([5, 3, 2]) + + grid_search.perform_grid_search( + model_parameters, inputs, labels, training_parameters, + filename_df, self.path + ) + + # testing saving the results + self.assertEqual(len(os.listdir(self.path)), 1) + + self.assertEqual(os.listdir(os.path.join(self.path, 'secir_simple_grid_search')), + [filename_df+".pickle"]) + + # testing size of the output + path_name = os.path.join(self.path, 'secir_simple_grid_search') + with open(os.path.join(path_name, filename_df + ".pickle"), "rb") as f: + d = pickle.load(f) + df = pd.DataFrame(d) + self.assertEqual(len(df['model']), 4) + if __name__ == '__main__': unittest.main()