%matplotlib inline
%reload_ext autoreload
%autoreload 2

bin_preds[source]

bin_preds()

Bins predictions over specified range

get_shape[source]

get_shape()

Extracts normalised shape of class from binned predictions. Empty bins are filled with a small quantity to avoid zeros.

get_paper_syst_shapes[source]

get_paper_syst_shapes()

Pass background data through trained model in order to get up/down shape variations.

get_likelihood_width[source]

get_likelihood_width(nll:ndarray, mu_scan:ndarray, val:float=0.5)

Compute width of likelihood at 95% confidence-level

interp_shape[source]

interp_shape(alpha:Tensor, f_nom:Tensor, f_up:Tensor, f_dw:Tensor)

Use quadratic interpolation between up/down systematic shapes and nominal in order to estimate shapes at arbitrary nuisance values. Linear extrapolation for absolute nuisances values greater than 1 (outside up/down shape range). Does not account for co-dependence of nuisances. Adapted from https://github.com/pablodecm/paper-inferno/blob/master/code/template_model.py under BSD 3-clause licence Copyright (c) 2018, Pablo de Castro, Tommaso Dorigo

def parallel_calc_nll(s_true:float, b_true:float, s_exp:Tensor, f_s:Tensor, alpha:Tensor,
             f_b_nom:Tensor, f_b_up:Tensor, f_b_dw:Tensor) -> Tensor:
    r'''Unused
    Compute multiple negative log-likelihood for specified parameters. Unused due to difficulty of batch-wise hessians in PyTorch.'''
    f_b = interp_shape(alpha, f_b_nom, f_b_up, f_b_dw)
    t_exp = (s_exp[:,None]*f_s[None,])+(b_true*f_b)
    asimov = (s_true*f_s)+(b_true*f_b_nom)
    p = torch.distributions.Poisson(t_exp, False)
    return -p.log_prob(asimov).sum(1)

def calc_diag_grad_hesse(nll:Tensor, alpha:Tensor) -> Tuple[Tensor,Tensor]:
    r'''Unused
    Compute batch-wise gradient and hessian, but only the diagonal elements.'''
    grad = autograd.grad(nll, alpha, torch.ones_like(nll, device=nll.device), create_graph=True)[0]
    hesse = autograd.grad(grad, alpha, torch.ones_like(alpha, device=nll.device), create_graph=True, retain_graph=True)[0]
    alpha.grad=None
    return grad, hesse

def calc_diag_profile(f_s:Tensor, f_b_nom:Tensor, f_b_up:Tensor, f_b_dw:Tensor, n:int,
                      mu_scan:Tensor, true_mu:int, n_steps:int=100, lr:float=0.1,  verbose:bool=True) -> Tensor:
    r'''Unused
    Compute profile likelihood for range of mu values, but only optimise using diagonal hessian elements.'''
    alpha = torch.zeros((len(mu_scan),f_b_up.shape[0]), requires_grad=True, device=f_b_nom.device)
    f_b_nom = f_b_nom.unsqueeze(0)
    get_nll = partialler(parallel_calc_nll, s_true=true_mu, b_true=n-true_mu, s_exp=mu_scan,
                         f_s=f_s, f_b_nom=f_b_nom, f_b_up=f_b_up, f_b_dw=f_b_dw)
    for i in range(n_steps):  # Newton optimise nuisances
        nll = get_nll(alpha=alpha)
        grad, hesse = calc_diag_grad_hesse(nll, alpha)
        step = torch.clamp(lr*grad.detach()/(hesse+1e-7), -100, 100)
        alpha = alpha-step
    return get_nll(alpha=alpha), alpha

calc_nll[source]

calc_nll(s_true:float, b_true:float, mu:Tensor, f_s_nom:Tensor, f_b_nom:Tensor, shape_alpha:Optional[Tensor]=None, s_norm_alpha:Optional[Tensor]=None, b_norm_alpha:Optional[Tensor]=None, f_s_up:Optional[Tensor]=None, f_s_dw:Optional[Tensor]=None, f_b_up:Optional[Tensor]=None, f_b_dw:Optional[Tensor]=None, s_norm_aux:Optional[Distribution]=None, b_norm_aux:Optional[Distribution]=None, shape_aux:Optional[List[Distribution]]=None)

Compute negative log-likelihood for specified parameters.

jacobian[source]

jacobian(y:Tensor, x:Tensor, create_graph=False)

Compute full jacobian matrix for single tensor. Call twice for hessian. Copied from https://gist.github.com/apaszke/226abdf867c4e9d6698bd198f3b45fb7 credits: Adam Paszke TODO: Fix this to work batch-wise (maybe https://gist.github.com/sbarratt/37356c46ad1350d4c30aefbd488a4faa)

calc_grad_hesse[source]

calc_grad_hesse(nll:Tensor, alpha:Tensor, create_graph:bool=False)

Compute full hessian and jacobian for single tensor

calc_profile[source]

calc_profile(f_s_nom:Tensor, f_b_nom:Tensor, n_obs:int, mu_scan:Tensor, mu_true:int, f_s_up:Optional[Tensor]=None, f_s_dw:Optional[Tensor]=None, f_b_up:Optional[Tensor]=None, f_b_dw:Optional[Tensor]=None, shape_aux:Optional[List[Distribution]]=None, s_norm_aux:Optional[List[Distribution]]=None, b_norm_aux:Optional[List[Distribution]]=None, nonaux_b_norm:bool=False, n_steps:int=100, lr:float=0.1, verbose:bool=True)

Compute profile likelihoods for range of mu values, optimising on full hessian. Ideally mu-values should be computed in parallel, but batch-wise hessian in PyTorch is difficult.