%matplotlib inline
%reload_ext autoreload
%autoreload 2
from pytorch_inferno.model_wrapper import ModelWrapper
from pytorch_inferno.callback import *
from pytorch_inferno.data import get_paper_data
from pytorch_inferno.plotting import *
from pytorch_inferno.inference import *
from pytorch_inferno.utils import *
from fastcore.all import partialler
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import itertools
from typing import *
from collections import OrderedDict
import torch.nn.functional as F
from torch import optim, autograd
import torch
from torch.distributions import Normal
bs = 2000
data, test = get_paper_data(200000, bs=bs, n_test=1000000)
x = torch.randn((1,10))
VariableSoftmax(0.1)(x), VariableSoftmax(0.5)(x), VariableSoftmax(1)(x)
net = nn.Sequential(nn.Linear(3,100), nn.ReLU(),
nn.Linear(100,100),nn.ReLU(),
nn.Linear(100,10), VariableSoftmax(0.1))
init_net(net)
x,y,w = next(iter(data.trn_dl))
preds = net(x)
assert preds.shape == (bs,10)
def to_shape(p:Tensor) -> Tensor:
f = p.sum(0)
f = f/f.sum()
return f
m = y.squeeze()==0
f_s = to_shape(preds[~m])
f_b = to_shape(preds[m])
plt.plot(to_np(f_s))
plt.plot(to_np(f_b))
Class initialised with a true number of signal events $\mu$ and a true number of events $n$. The true number of background events is then computed as $n-mu$. The number of nuisance parameters must also be passed.
This process includes the analytical effect of the nuisances, and their interdependence, since the augmented batch is passed through the model. This does, however, require that the effect of each nuisance on the input features is accurately parameterised.
Running cycle is:
- on_train_begin: set alpha as tensor with gradient and
n_alpha+1elements. Set element zero to $\mu$. - for e in #epochs:
- for b in #batches:
- on_batch_begin: add relevant nuisances to incoming batch of inputs (as implemented in inheriting function)
- on_forwards_end:
- compute the full hessian of the negative log-likelihood w.r.t alpha (shape systematics included since they're added to the input data and passed through the model)
- invert hessian , and set the diagonal element corresponding to $\mu$ as the loss value
- on_batch_end: zero gradient of alpha
- for b in #batches:
net = nn.Sequential(nn.Linear(3,100), nn.ReLU(),
nn.Linear(100,100),nn.ReLU(),
nn.Linear(100,10), VariableSoftmax(0.1))
init_net(net)
model = ModelWrapper(net)
%%time
model.fit(1, data=data, opt=partialler(optim.Adam,lr=1e-3), loss=None,
cbs=[PaperInferno(float_r=True, float_l=True, shape_aux=[Normal(0,2), Normal(0,2)], nonaux_b_norm=False, b_norm_aux=[Normal(0,100)]),
LossTracker(),SaveBest('weights/best_ie.h5'),EarlyStopping(2)])
model.save('weights/Inferno_Test_exact_bm2.h5')
model.load('weights/Inferno_Test_exact_bm2.h5')
preds = model._predict_dl(test, pred_cb=InfernoPred())
df = pd.DataFrame({'pred':preds})
df['gen_target'] = test.dataset.y
df.head()
plot_preds(df, bin_edges=np.linspace(0,10,11))
bin_preds(df)
df.head()
f_s,f_b = get_shape(df,1),get_shape(df,0)
assert is_close(f_s.sum().data.item(),f_b.sum().data.item()) and is_close(f_b.sum().data.item(), 1.0)
asimov = (50*f_s)+(1000*f_b)
asimov, asimov.sum()
n = 1050
mu_scan = torch.linspace(20,80,61)
y = torch.zeros_like(mu_scan)
for i,m in enumerate(mu_scan):
pois = torch.distributions.Poisson((m*f_s)+(1000*f_b), False)
y[i] = -pois.log_prob(asimov).sum()
y
w1 = plot_likelihood(y, mu_scan=mu_scan)
profiler = partialler(calc_profile, f_s_nom=f_s, f_b_nom=f_b, n_obs=1050, mu_scan=mu_scan, mu_true=50, n_steps=100)
nll = profiler()
w2 = plot_likelihood(nll, mu_scan=mu_scan)
assert w1 == w2
bkg = test.dataset.x[test.dataset.y.squeeze() == 0]
assert len(bkg) == 500000
b_shapes = get_paper_syst_shapes(bkg, df, model=model, pred_cb=InfernoPred())
df
plot_preds(df, pred_names=['pred', 'pred_-0.2_3', 'pred_0.2_3', 'pred_0_2.5', 'pred_0_3.5'], bin_edges=np.linspace(0,10,11))
nll = profiler(f_b_up=b_shapes['f_b_up'][0], f_b_dw=b_shapes['f_b_dw'][0])
plot_likelihood(nll, mu_scan)
nll = profiler(f_b_up= b_shapes['f_b_up'][1], f_b_dw=b_shapes['f_b_dw'][1])
plot_likelihood(nll, mu_scan)
nll = profiler(**b_shapes)
plot_likelihood(nll, mu_scan)
b_shapes
f_s
shape_aux = [Normal(0,2), Normal(0,2)]
nll = profiler(shape_aux=shape_aux, **b_shapes)
plot_likelihood(nll, mu_scan)
nll = profiler(shape_aux=shape_aux, b_norm_aux=[Normal(0,100)], **b_shapes)
plot_likelihood(nll, mu_scan)
