Automatic Differentiation
All code in Newtrinos.jl is fully differentiable via AutoDiff, using ForwardDiff.jl under the hood. This means you can compute exact gradients of the joint log-likelihood with respect to all physics and nuisance parameters — which is essential for gradient-based optimizers and samplers.
Configuration
We must first configure the experiments (with default physics):
using Newtrinos, DensityInterface
experiments = (
deepcore = Newtrinos.deepcore.configure(),
dayabay = Newtrinos.dayabay.configure(),
kamland = Newtrinos.kamland.configure(),
minos = Newtrinos.minos.configure(),
orca = Newtrinos.orca.configure(),
);This is enough to generate a joint likelihood and collect the nominal parameters for everything:
likelihood = Newtrinos.generate_likelihood(experiments)
params = Newtrinos.get_params(experiments)
priors = Newtrinos.get_params(experiments)Computing gradients
In a neutrino fit the gradient is useful for:
- Gradient-based optimizers (e.g. L-BFGS) – finding the best-fit point much faster than grid scans
- HMC/NUTS samplers – needed for efficient Bayesian posterior sampling via BAT.jl
- Sensitivity studies – large $|\partial \text{llh} / \partial \theta |$ means the data is highly sensitive to that parameter
We can calculate gradients directly using ForwardDiff
using ForwardDiff
# Wrap the likelihood as a plain function of parameters
f(params) = logdensityof(likelihood, params)
# Evaluate the gradient at nominal parameter values
grad = ForwardDiff.gradient(f, params)The result is a NamedTuple with the same fields as p, giving you ∂llh/∂θ for every parameter:
grad.δCP # ∂llh/∂δCP
grad.θ₂₃ # ∂llh/∂θ₂₃
grad.Δm²₃₁ # ∂llh/∂Δm²₃₁The gradient calculation can be automatically activated in likelihood scans by setting gradient_map=true. Then, at each grid point the ForwardDiff.gradient is evaluated and the per-parameter gradient is included in the NewtrinosResult.values.
using DataStructures
vars_to_scan=OrderedDict(:θ₂₃ => 21,)
result=Newtrinos.scan(likelihood, priors, vars_to_scan, params, gradient_map=true)
#acces gradient of parameter via result.values.paramter_grad
result.values.θ₂₃_grad