MadMiner physics tutorial (part 3B)

Johann Brehmer, Felix Kling, Irina Espejo, and Kyle Cranmer 2018-2019

In part 3A of this tutorial we will finally train a neural network to estimate likelihood ratios. 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.


Make sure you’ve run the first tutorial before executing this notebook!

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 import ScoreEstimator
# MadMiner output
    format='%(asctime)-5.5s %(name)-20.20s %(levelname)-7.7s %(message)s',

# Output of all other modules (e.g. matplotlib)
for key in logging.Logger.manager.loggerDict:
    if "madminer" not in key:

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.

sampler = SampleAugmenter('data/lhe_data_shuffled.h5')
# sampler = SampleAugmenter('data/delphes_data_shuffled.h5')
19:25 madminer.analysis    INFO    Loading data from data/lhe_data_shuffled.h5
19:25 madminer.analysis    INFO    Found 2 parameters
19:25 madminer.analysis    INFO    Did not find nuisance parameters
19:25 madminer.analysis    INFO    Found 6 benchmarks, of which 6 physical
19:25 madminer.analysis    INFO    Found 3 observables
19:25 madminer.analysis    INFO    Found 89991 events
19:25 madminer.analysis    INFO      49991 signal events sampled from benchmark sm
19:25 madminer.analysis    INFO      10000 signal events sampled from benchmark w
19:25 madminer.analysis    INFO      10000 signal events sampled from benchmark neg_w
19:25 madminer.analysis    INFO      10000 signal events sampled from benchmark ww
19:25 madminer.analysis    INFO      10000 signal events sampled from benchmark neg_ww
19:25 madminer.analysis    INFO    Found morphing setup with 6 components
19:25 madminer.analysis    INFO    Did not find nuisance morphing setup

The relevant SampleAugmenter function for local score estimators is extract_samples_train_local(). As in part 3a of the tutorial, for the argument theta you can use the helper functions sampling.benchmark(), sampling.benchmarks(), sampling.morphing_point(), sampling.morphing_points(), and sampling.random_morphing_points().

x, theta, t_xz, _ = sampler.sample_train_local(
19:25 madminer.sampling    INFO    Extracting training sample for local score regression. Sampling and score evaluation according to sm
19:25 madminer.sampling    INFO    Starting sampling serially
19:25 madminer.sampling    INFO    Sampling from parameter point 1 / 1
19:25 madminer.sampling    INFO    Effective number of samples: mean 29966.000000000004, with individual thetas ranging from 29966.000000000004 to 29966.000000000004

We can use the same data as in part 3a, so you only have to execute this if you haven’t gone through tutorial 3a:

_ = sampler.sample_test(
19:25 madminer.sampling    INFO    Extracting evaluation sample. Sampling according to sm
19:25 madminer.sampling    INFO    Starting sampling serially
19:25 madminer.sampling    INFO    Sampling from parameter point 1 / 1
19:25 madminer.sampling    INFO    Effective number of samples: mean 10097.0, with individual thetas ranging from 10097.0 to 10097.0

2. Train score estimator

It’s now time to build a neural network. Only this time, instead of the likelihood ratio itself, we will estimate the gradient of the log likelihood with respect to the theory parameters – the score. To be precise, the output of the neural network is an estimate of the score at some reference parameter point, for instance the Standard Model. A neural network that estimates this “local” score can be used to calculate the Fisher information at that point. The estimated score can also be used as a machine learning version of Optimal Observables, and likelihoods can be estimated based on density estimation in the estimated score space. This method for likelihood ratio estimation is called SALLY, and there is a closely related version called SALLINO. Both are explained in “Constraining Effective Field Theories With Machine Learning” and “A Guide to Constraining Effective Field Theories With Machine Learning”.

The central object for this is the class:

estimator = ScoreEstimator(n_hidden=(30,30))
19:25          INFO    Starting training
19:25          INFO      Batch size:             128
19:25          INFO      Optimizer:              amsgrad
19:25          INFO      Epochs:                 50
19:25          INFO      Learning rate:          0.001 initially, decaying to 0.0001
19:25          INFO      Validation split:       0.25
19:25          INFO      Early stopping:         True
19:25          INFO      Scale inputs:           True
19:25          INFO      Shuffle labels          False
19:25          INFO      Samples:                all
19:25          INFO    Loading training data
19:25 madminer.utils.vario INFO      Loading data/samples/x_train_score.npy into RAM
19:25 madminer.utils.vario INFO      Loading data/samples/t_xz_train_score.npy into RAM
19:25          INFO    Found 500000 samples with 2 parameters and 3 observables
19:25          INFO    Setting up input rescaling
19:25          INFO    Creating model
19:25          INFO    Training model
19:25 INFO    Training on CPU with single precision
19:26 INFO      Epoch   3: train loss  0.24239 (mse_score:  0.242)
19:26 INFO                 val. loss   0.25060 (mse_score:  0.251)
19:27 INFO      Epoch   6: train loss  0.22593 (mse_score:  0.226)
19:27 INFO                 val. loss   0.23032 (mse_score:  0.230)
19:28 INFO      Epoch   9: train loss  0.21441 (mse_score:  0.214)
19:28 INFO                 val. loss   0.21741 (mse_score:  0.217)
19:29 INFO      Epoch  12: train loss  0.20653 (mse_score:  0.207)
19:29 INFO                 val. loss   0.20709 (mse_score:  0.207)
19:30 INFO      Epoch  15: train loss  0.19986 (mse_score:  0.200)
19:30 INFO                 val. loss   0.19887 (mse_score:  0.199)
19:31 INFO      Epoch  18: train loss  0.19441 (mse_score:  0.194)
19:31 INFO                 val. loss   0.19143 (mse_score:  0.191)
19:32 INFO      Epoch  21: train loss  0.19016 (mse_score:  0.190)
19:32 INFO                 val. loss   0.18850 (mse_score:  0.189)
19:33 INFO      Epoch  24: train loss  0.18647 (mse_score:  0.186)
19:33 INFO                 val. loss   0.18236 (mse_score:  0.182)
19:35 INFO      Epoch  27: train loss  0.18417 (mse_score:  0.184)
19:35 INFO                 val. loss   0.18023 (mse_score:  0.180)
19:36 INFO      Epoch  30: train loss  0.18195 (mse_score:  0.182)
19:36 INFO                 val. loss   0.17799 (mse_score:  0.178)
19:37 INFO      Epoch  33: train loss  0.17990 (mse_score:  0.180)
19:37 INFO                 val. loss   0.17600 (mse_score:  0.176)
19:38 INFO      Epoch  36: train loss  0.17840 (mse_score:  0.178)
19:38 INFO                 val. loss   0.17421 (mse_score:  0.174)
19:39 INFO      Epoch  39: train loss  0.17712 (mse_score:  0.177)
19:39 INFO                 val. loss   0.17329 (mse_score:  0.173)
19:40 INFO      Epoch  42: train loss  0.17623 (mse_score:  0.176)
19:40 INFO                 val. loss   0.17126 (mse_score:  0.171)
19:41 INFO      Epoch  45: train loss  0.17538 (mse_score:  0.175)
19:41 INFO                 val. loss   0.17093 (mse_score:  0.171)
19:42 INFO      Epoch  48: train loss  0.17472 (mse_score:  0.175)
19:42 INFO                 val. loss   0.17005 (mse_score:  0.170)
19:43 INFO    Early stopping after epoch 49, with loss  0.16998 compared to final loss  0.17045
19:43 INFO    Training time spend on:
19:43 INFO                      initialize model:   0.00h
19:43 INFO                                   ALL:   0.31h
19:43 INFO                            check data:   0.00h
19:43 INFO                          make dataset:   0.00h
19:43 INFO                       make dataloader:   0.00h
19:43 INFO                       setup optimizer:   0.00h
19:43 INFO                   initialize training:   0.00h
19:43 INFO                                set lr:   0.00h
19:43 INFO                   load training batch:   0.10h
19:43 INFO                        fwd: move data:   0.00h
19:43 INFO                   fwd: check for nans:   0.03h
19:43 INFO                    fwd: model.forward:   0.05h
19:43 INFO                 fwd: calculate losses:   0.01h
19:43 INFO                 training forward pass:   0.07h
19:43 INFO                   training sum losses:   0.01h
19:43 INFO                        opt: zero grad:   0.00h
19:43 INFO                         opt: backward:   0.04h
19:43 INFO                   opt: clip grad norm:   0.00h
19:43 INFO                             opt: step:   0.02h
19:43 INFO                        optimizer step:   0.06h
19:43 INFO                 load validation batch:   0.04h
19:43 INFO               validation forward pass:   0.02h
19:43 INFO                 validation sum losses:   0.00h
19:43 INFO                        early stopping:   0.00h
19:43 INFO                          report epoch:   0.00h
19:43          INFO    Saving model to models/sally

3. Evaluate score estimator

Let’s evaluate the SM score on the test data


t_hat = estimator.evaluate_score(
19:43          INFO    Loading model from models/sally
19:43 madminer.utils.vario INFO      Loading data/samples/x_test.npy into RAM
19:43          INFO    Starting score evaluation

Let’s have a look at the estimated score and how it is related to the observables:

x = np.load('data/samples/x_test.npy')

fig = plt.figure(figsize=(10,4))

for i in range(2):
    ax = plt.subplot(1,2,i+1)

    sc = plt.scatter(x[:,0], x[:,1], c=t_hat[:,i], s=25., cmap='viridis', vmin=-1., vmax=1.)
    cbar = plt.colorbar(sc)

    cbar.set_label(r'$\hat{t}_' + str(i) + r'(x | \theta_{ref})$')
    plt.xlabel(r'$p_{T,j1}$ [GeV]')
    plt.ylabel(r'$\Delta \phi_{jj}$')