<jupyter_start><jupyter_text>Problem parameters<jupyter_code># Load problem from predefined ones
model = emupints.Problems.LogisticModel
problem = emupints.Problems.load_problem(model)

problem['values'][:5]
toy.HodgkinHuxleyIKModel().suggested_times()
# take all the required variables from the problem
n_parameters = problem['n_parameters']
log_likelihood = problem['log_likelihood']
log_prior = problem['log_prior']
log_posterior = problem['log_posterior']
real_parameters = problem['parameters']
pints_problem = problem['problem']
bounds = problem['bounds']
index_to_param_name = problem['param_names']<jupyter_output><empty_output><jupyter_text>Apply CMA-ES<jupyter_code># CMA-es
score = pints.SumOfSquaresError(pints_problem)
x0 = bounds.sample(1)[0]

found_parameters, found_value = pints.optimise(
    score,
    x0,
    None,
    bounds,
    method=pints.CMAES
)<jupyter_output>Minimising error measure
using Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
Running in sequential mode.
Population size: 6
Iter. Eval. Best      Time m:s
0     6      512167.8   0:00.1
1     12     512167.8   0:00.1
2     18     511346.8   0:00.1
3     24     511342.4   0:00.1
20    126    510879     0:00.1
40    246    473774.9   0:00.1
60    366    120968.9   0:00.2
80    486    120574.6   0:00.2
100   606    120574.6   0:00.2
120   726    120574.6   0:00.2
140   846    120574.6   0:00.3
160   966    120574.6   0:00.3
180   1086   120574.6   0:00.3
200   1206   120574.6   0:00.4
220   1326   120574.6   0:00.4
240   1446   120574.6   0:00.4
260   1566   120574.6   0:00.5
280   1686   120574.6   0:00.5
300   1806   120574.6   0:00.5
320   1926   120574.6   0:00.6
340   2046   120574.6   0:00.6
360   2166   120574.6   0:00.6
380   2286   120574.6   0:00.6
400   2406   120574.6   0:00.7
420   2526   120574.6   0:00.7
440   2646   120574.6   0:00.7
460   2766   120574.6   0:00[...]<jupyter_text>Use found parameters to sample trianing data<jupyter_code># training size
# use bigger training size for NNs because speed stays the same
train_size = 1000
test_size = 400
sampling_std = 0.05 * found_parameters # take 0.05% around found parameters std

def normal_sampling(size, mean, std):
    return mean + np.random.randn(*size) * std

# create data for training and testing
# genera example using normal distribution around parametrs
train_X = normal_sampling((train_size, n_parameters), found_parameters, sampling_std)
train_y = np.apply_along_axis(log_likelihood, 1, train_X)

# if provide validation data to keras model also need to normalize it 
test_X = normal_sampling((test_size, n_parameters), found_parameters, sampling_std)
test_y = np.apply_along_axis(log_likelihood, 1, test_X)<jupyter_output><empty_output><jupyter_text>Create NN emulator<jupyter_code>EPOCHS = 50
BATCH = 16
validation_split = 0.2
# how to scale data
# for input use standardisation, same as GPs
# for output scale linearly in range (-1, 1)
input_scaler = MinMaxScaler(feature_range=(-1, 1))
output_scaler = StandardScaler()
emu = emupints.NNEmulator(log_likelihood, 
                          train_X,
                          train_y,
                          model_size='small',
                          input_scaler=input_scaler, 
                          output_scaler=output_scaler
                         )
# print summary of the model
emu.summary()
# store best
mc_callback = keras.callbacks.ModelCheckpoint("models/model_{epoch:02d}.h5", 
                                              save_weights_only=True,  
                                             )
# store the results during training
val_umae = []

# Callback to calculate unormalized mae
class UnormalizedMAE(keras.callbacks.Callback):

    def on_epoch_end(self, batch, logs={}):
        X_val, y_val = self.validation_data[0], self.validation_data[1]
        
        y_val_predict = np.asarray(self.model.predict(X_val))

        # unormalize predictions
        y_val = output_scaler.inverse_transform(y_val)
        y_val_predict = output_scaler.inverse_transform(y_val_predict)

        # calculate mae
        mae_val = np.mean(np.abs(y_val - y_val_predict))

        # store results
        val_umae.append(mae_val)

        return
#optimizer = tf.train.RMSPropOptimizer(learning_rate)
emu.set_parameters(
    loss='mse',
    optimizer='adam',
    metrics=['mae'] # mean absolute error
)
history = emu.fit(epochs=EPOCHS,
                  batch_size=BATCH,
                  validation_split=0.2,
                  callbacks=[mc_callback,UnormalizedMAE()],
                 )
best_model_idx = np.argmin(val_umae)
emu._model.load_weights("models/model_{:02d}.h5".format(best_model_idx))
plt.figure()
plt.xlabel('Epoch')
plt.ylabel('Mean Abs Error')
plt.plot(history.epoch, np.array(history.history['mean_absolute_error']),
       label='Train mae')
plt.plot(history.epoch, np.array(history.history['val_mean_absolute_error']),
       label = 'Val mae')
plt.legend()
plt.figure()
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.plot(history.epoch, np.array(history.history['loss']),
       label='Train Loss')
plt.plot(history.epoch, np.array(history.history['val_loss']),
       label = 'Val loss')
plt.legend()
plt.figure()
plt.xlabel('Epoch')
plt.ylabel('Unormalized mae')
plt.plot(history.epoch, val_umae)<jupyter_output><empty_output><jupyter_text>Likelihood metrics<jupyter_code>emu_time = %timeit -r 50 -o emu(real_parameters)
print("Emulator prediction time (micro-seconds): {:5f}".format(emu_time.average * 1000000))
likelihood_time = %timeit -r 50 -o log_likelihood(real_parameters)
print("Likelihood prediction time (micro-seconds): {:5f}".format(likelihood_time.average * 1000000))
pred_y = np.apply_along_axis(emu, 1, test_X)
print("Mean absolute error (mae): {:.5f}".format(emumet.mae(test_y, pred_y)))
print("Mean percentage error (mape): {:.5f}".format(emumet.mape(test_y, pred_y) * 100))<jupyter_output>Mean percentage error (mape): 8.64502<jupyter_text>Comparing MCMC run times<jupyter_code>mcmc_method = pints.AdaptiveCovarianceMCMC

# MCMC parameters
num_chains = 3
max_iters = 10000
# Create a posterior log-likelihood (log(likelihood * prior))
emu_log_posterior = pints.LogPosterior(emu, log_prior)
xs = [
    found_parameters,
    found_parameters * 1.01,
    found_parameters * 0.99
]
xs
emu_mcmc = pints.MCMCSampling(emu_log_posterior,
                              num_chains,
                              xs,
                              method = mcmc_method,
                             )
emu_mcmc.set_log_to_screen(False)
print('Running...')
emu_chains = emu_mcmc.run()
print('Done!')
real_mcmc = pints.MCMCSampling(log_posterior, 
                               num_chains, 
                               xs, 
                               method = mcmc_method,
                               )
real_mcmc.set_max_iterations(max_iters)
real_mcmc.set_log_to_screen(False)
# Run!
print('Running...')
real_chains = real_mcmc.run()
print('Done!')
pints.plot.trace(emu_chains)
plt.show()
pints.plot.trace(real_chains)
plt.show()<jupyter_output><empty_output><jupyter_text>MCMC metrics<jupyter_code>print('Emulator R-hat:')
print(pints.rhat_all_params(emu_chains))
print('Real R-hat:')
print(pints.rhat_all_params(real_chains))
# percentage
chain = real_chains[0]

chain_mape = emumet.chain_mape(chain, emu, log_likelihood)
print("Chain mean absolute percentage error: {:.5f}".format(chain_mape * 100))
emu_mcmc_time = %timeit -o emu_mcmc.run()
print("Emulator MCMC time: {:.5f}".format(emu_mcmc_time.average))
real_mcmc_time = %timeit -o real_mcmc.run()
print("MCMC time: {:.5f}".format(real_mcmc_time.average))
emu_parameters, emu_parameters_std = emumet.estimate_parameters(emu_chains)
print("Emulator estimated parameters: ", emu_parameters)
print("Emulator parametet std: ", emu_parameters_std)
real_parameters, real_parameters_std = emumet.estimate_parameters(real_chains)
print("Estimated parameters: ", real_parameters)
print("parametet std: ", real_parameters_std)<jupyter_output>Estimated parameters:  [1.49854152e-01 4.99286823e+02]
parametet std:  [9.68200776e-06 3.09614573e-02]<jupyter_text>Plotting emulator surfaces around parameters<jupyter_code># use smaller bounds to get an idea of the surface around real parameters
min_bounds = np.min(chain, axis = 0)
max_bounds = np.max(chain, axis = 0)
bounds = pints.RectangularBoundaries(min_bounds, max_bounds)
bounds
# plot surfaces
if n_parameters >= 3:
    fixed_parameters = emutils.fix_parameters(bounds)
    
    fig, ax = emuplt.plot_fixed_param_grid(
        emu,
        fixed_parameters,
        bounds,
        shape = (n_parameters, n_parameters - 1),
        contour = False,
        index_to_param_name=index_to_param_name,
        additional_log_likelihoods = [log_likelihood]
    )

    plt.show(fig)
    
# for row i we have parameter i fixed and some other parameter j, depending on the column, j != i
print(bounds.lower())
print(bounds.upper())
fig, ax = emuplt.plot_surface_fixed_param(
    emu,
    bounds,
    fixed = [(1, 0.29)],
    index_to_param_name=index_to_param_name,
    contour=False,
    precision=2,
    additional_log_likelihoods=[log_likelihood]
)
ax.view_init(30, 45 + 90 *3)
fig.set_size_inches(10.5, 6)
fig<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10899241603466955, 'fraction_numerical': 0.06619718309859154, 'mean_word_length': 3.4788937409024747, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '91254', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3678, 'n_tokens_neox': 3246, 'n_words': 655, 'path': 'studies/NNExperiments/Logistic.ipynb', 'repo_name': 'MaksymPetyak/emupints', 'size': '9230'}