<jupyter_start><jupyter_text>MadMiner physics tutorial (part 3C)<NAME>, <NAME>, <NAME>, and <NAME> 2018-2019 In part 3C of this tutorial we will train a third neural estimator: this time of the likelihood function itself (rather than its ratio). We assume that you have run part 1 and 2A of this tutorial. If, instead of 2A, you have run part 2B, you just have to load a different filename later.  Preparations<jupyter_code>import logging
import numpy as np
import matplotlib
from matplotlib import pyplot as plt

%matplotlib inline

from madminer.sampling import SampleAugmenter
from madminer import sampling
from madminer.ml import LikelihoodEstimator
# MadMiner output
logging.basicConfig(
    format="%(asctime)-5.5s %(name)-20.20s %(levelname)-7.7s %(message)s",
    datefmt="%H:%M",
    level=logging.INFO,
)

# Output of all other modules (e.g. matplotlib)
for key in logging.Logger.manager.loggerDict:
    if "madminer" not in key:
        logging.getLogger(key).setLevel(logging.WARNING)<jupyter_output><empty_output><jupyter_text>1. Make (unweighted) training and test samples with augmented data At this point, we have all the information we need from the simulations. But the data is not quite ready to be used for machine learning. The `madminer.sampling` class `SampleAugmenter` will take care of the remaining book-keeping steps before we can train our estimators:First, it unweights the samples, i.e. for a given parameter vector `theta` (or a distribution `p(theta)`) it picks events `x` such that their distribution follows `p(x|theta)`. The selected samples will all come from the event file we have so far, but their frequency is changed -- some events will appear multiple times, some will disappear.Second, `SampleAugmenter` calculates all the augmented data ("gold") that is the key to our new inference methods. Depending on the specific technique, these are the joint likelihood ratio and / or the joint score. It saves all these pieces of information for the selected events in a set of numpy files that can easily be used in any machine learning framework.<jupyter_code>sampler = SampleAugmenter("data/lhe_data_shuffled.h5")
# sampler = SampleAugmenter('data/delphes_data_shuffled.h5')<jupyter_output>03:42 madminer.analysis    INFO    Loading data from data/lhe_data_shuffled.h5
03:42 madminer.analysis    INFO    Found 2 parameters
03:42 madminer.analysis    INFO    Did not find nuisance parameters
03:42 madminer.analysis    INFO    Found 6 benchmarks, of which 6 physical
03:42 madminer.analysis    INFO    Found 3 observables
03:42 madminer.analysis    INFO    Found 89991 events
03:42 madminer.analysis    INFO      49991 signal events sampled from benchmark sm
03:42 madminer.analysis    INFO      10000 signal events sampled from benchmark w
03:42 madminer.analysis    INFO      10000 signal events sampled from benchmark neg_w
03:42 madminer.analysis    INFO      10000 signal events sampled from benchmark ww
03:42 madminer.analysis    INFO      10000 signal events sampled from benchmark neg_ww
03:42 madminer.analysis    INFO    Found morphing setup with 6 components
03:42 madminer.analysis    INFO    Did not find nuisance morphing setup<jupyter_text>The `SampleAugmenter` class defines five different high-level functions to generate train or test samples:- `sample_train_plain()`, which only saves observations x, for instance for histograms or ABC;- `sample_train_local()` for methods like SALLY and SALLINO, which will be demonstrated in the second part of the tutorial;- `sample_train_density()` for neural density estimation techniques like MAF or SCANDAL;- `sample_train_ratio()` for techniques like CARL, ROLR, CASCAL, and RASCAL, when only theta0 is parameterized;- `sample_train_more_ratios()` for the same techniques, but with both theta0 and theta1 parameterized;- `sample_test()` for the evaluation of any method.For the arguments `theta`, `theta0`, or `theta1`, you can (and should!) use the helper functions `benchmark()`, `benchmarks()`, `morphing_point()`, `morphing_points()`, and `random_morphing_points()`, all defined in the `madminer.sampling` module.Here we'll train a likelihood estimator with the SCANDAL method, so we focus on the `extract_samples_train_density()` function. We'll sample the numerator hypothesis in the likelihood ratio with 1000 points drawn from a Gaussian prior, and fix the denominator hypothesis to the SM.Note the keyword `sample_only_from_closest_benchmark=True`, which makes sure that for each parameter point we only use the events that were originally (in MG) generated from the closest benchmark. This reduces the statistical fluctuations in the outcome quite a bit.<jupyter_code>x, theta, t_xz, _ = sampler.sample_train_density(
    theta=sampling.random_morphing_points(10000, [("gaussian", 0.0, 0.5), ("gaussian", 0.0, 0.5)]),
    n_samples=500000,
    folder="./data/samples",
    filename="train_density",
    sample_only_from_closest_benchmark=True,
)<jupyter_output>03:42 madminer.sampling    INFO    Extracting training sample for non-local score-based methods. Sampling and score evaluation according to (u'random_morphing_points', (10000, [(u'gaussian', 0.0, 0.5), (u'gaussian', 0.0, 0.5)]))
03:42 madminer.sampling    INFO    Starting sampling serially
03:42 madminer.sampling    WARNING Large statistical uncertainty on the total cross section when sampling from theta = [0.58710176 1.2685223 ]: (0.001242 +/- 0.000441) pb (35.5351101412644 %). Skipping these warnings in the future...
03:43 madminer.sampling    INFO    Sampling from parameter point 500 / 10000
03:43 madminer.sampling    INFO    Sampling from parameter point 1000 / 10000
03:44 madminer.sampling    INFO    Sampling from parameter point 1500 / 10000
03:45 madminer.sampling    INFO    Sampling from parameter point 2000 / 10000
03:45 madminer.sampling    INFO    Sampling from parameter point 2500 / 10000
03:47 madminer.sampling    INFO    Sampling from parameter point 3000 / 10000
03:48 ma[...]<jupyter_text>For the evaluation we'll need a test sample:<jupyter_code>_ = sampler.sample_test(
    theta=sampling.benchmark("sm"),
    n_samples=1000,
    folder="./data/samples",
    filename="test",
)<jupyter_output>05:18 madminer.sampling    INFO    Extracting evaluation sample. Sampling according to sm
05:18 madminer.sampling    INFO    Starting sampling serially
05:18 madminer.sampling    INFO    Sampling from parameter point 1 / 1
05:18 madminer.sampling    INFO    Effective number of samples: mean 10015.0, with individual thetas ranging from 10015.0 to 10015.0<jupyter_text>2. Train likelihood estimator It's now time to build the neural network that estimates the likelihood ratio. The central object for this is the `madminer.ml.ParameterizedRatioEstimator` class. It defines functions that train, save, load, and evaluate the estimators.In the initialization, the keywords `n_hidden` and `activation` define the architecture of the (fully connected) neural network:<jupyter_code>estimator = LikelihoodEstimator(n_mades=3, n_hidden=(60, 60), activation="relu")<jupyter_output><empty_output><jupyter_text>To train this model we will minimize the SCANDAL loss function described in ["Mining gold from implicit models to improve likelihood-free inference"](https://arxiv.org/abs/1805.12244).<jupyter_code>estimator.train(
    method="scandal",
    theta="data/samples/theta_train_density.npy",
    x="data/samples/x_train_density.npy",
    t_xz="data/samples/t_xz_train_density.npy",
    alpha=10.0,
    n_epochs=10,
    initial_lr=0.01,
)

estimator.save("models/scandal")<jupyter_output>05:18 madminer.ml          INFO    Starting training
05:18 madminer.ml          INFO      Method:                 scandal
05:18 madminer.ml          INFO      alpha:                  10.0
05:18 madminer.ml          INFO      Batch size:             128
05:18 madminer.ml          INFO      Optimizer:              amsgrad
05:18 madminer.ml          INFO      Epochs:                 10
05:18 madminer.ml          INFO      Learning rate:          0.01 initially, decaying to 0.0001
05:18 madminer.ml          INFO      Validation split:       0.25
05:18 madminer.ml          INFO      Early stopping:         True
05:18 madminer.ml          INFO      Scale inputs:           True
05:18 madminer.ml          INFO      Shuffle labels          False
05:18 madminer.ml          INFO      Samples:                all
05:18 madminer.ml          INFO    Loading training data
05:18 madminer.utils.vario INFO      Loading data/samples/theta_train_density.npy into RAM
05:18 madminer.utils.vario INFO      Loa[...]<jupyter_text>(nll:  3.830, mse_score:  0.064  3. Evaluate likelihood estimator `estimator.evaluate_log_likelihood(theta,x)` estimated the log likelihood for all combination between the given phase-space points `x` and parameters `theta`. That is, if given 100 events `x` and a grid of 25 `theta` points, it will return 25\*100 estimates for the log likelihood, indexed by `[i_theta,i_x]`.<jupyter_code>theta_each = np.linspace(-1.0, 1.0, 25)
theta0, theta1 = np.meshgrid(theta_each, theta_each)
theta_grid = np.vstack((theta0.flatten(), theta1.flatten())).T
np.save("data/samples/theta_grid.npy", theta_grid)
estimator.load("models/scandal")

log_p_hat, _ = estimator.evaluate_log_likelihood(
    theta="data/samples/theta_grid.npy",
    x="data/samples/x_test.npy",
    evaluate_score=False,
)<jupyter_output>05:28 madminer.ml          INFO    Loading model from models/scandal
05:28 madminer.utils.vario INFO      Loading data/samples/theta_grid.npy into RAM
05:28 madminer.utils.vario INFO      Loading data/samples/x_test.npy into RAM
05:28 madminer.ml          INFO    Starting ratio evaluation for 625000 x-theta combinations
05:28 madminer.ml          INFO    Evaluation done<jupyter_text>Let's look at the result:<jupyter_code>bin_size = theta_each[1] - theta_each[0]
edges = np.linspace(
    theta_each[0] - bin_size / 2,
    theta_each[-1] + bin_size / 2,
    len(theta_each) + 1,
)

fig = plt.figure(figsize=(6, 5))
ax = plt.gca()

expected_llr = np.mean(log_p_hat, axis=1)
best_fit = theta_grid[np.argmin(-2.0 * expected_llr)]

cmin = np.min(-2 * expected_llr)
cmax = np.max(-2 * expected_llr)

pcm = ax.pcolormesh(
    edges,
    edges,
    -2.0 * expected_llr.reshape((25, 25)),
    norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
    cmap="viridis_r",
)

cbar = fig.colorbar(pcm, ax=ax, extend="both")

plt.scatter(
    best_fit[0],
    best_fit[1],
    s=80.0,
    color="black",
    marker="*",
)

plt.xlabel(r"$\theta_0$")
plt.ylabel(r"$\theta_1$")
cbar.set_label(r"$\mathbb{E}_x [ -2\, \log \,\hat{r}(x | \theta, \theta_{SM}) ]$")

plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08858985141575554, 'fraction_numerical': 0.04840669096346136, 'mean_word_length': 4.031499764927127, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '420008', 'license': 'MIT', 'n_tokens_mistral': 3661, 'n_tokens_neox': 3223, 'n_words': 1106, 'path': 'examples/tutorial_particle_physics/3c_likelihood.ipynb', 'repo_name': 'diana-hep/madminer', 'size': '10726'}