Inference using Dynesty (NestedSampling)

In this example we look at how to use a random sampler for inference purposes with Finesse. This is a method that can also be used to infer optical model parameters from measured data. Although here we generate simulated data and compare against some known true values.

Note

Running this complete example takes too long to be included in Finesse’s documentation build. Therefore only the first two cells are ran and are guaranteed to run successfully. The other cells containing the dynesty code are included statically and are not guaranteed to run successfully in the future.

import os
import numpy as np

import finesse
from scipy.stats import truncnorm
import numpy as np
import matplotlib.pyplot as plt
import finesse
from finesse.knm import Map
from finesse.utilities.maps import circular_aperture

Define the Finesse model, parameters. We will generate a mock signal: the power as a function of detuning phase (a cavity scan). We will add some realistic measurement noise on a photodiode in the next cell

# create the finesse model with a fixed RoC, and beam offsets xbeta and ybeta

finesse.init_plotting()

model = finesse.Model()
model.parse(
    """
    l l1
    mod mod1 f=9.1M midx=0.1
    m m1 R=0.984 T=0.014 Rc=-1940 xbeta=0.1e-6 ybeta=0.1e-6 ## treat RoC, betas etc... as unknowns for MCMC
    m m2 R=1 T=0 Rc=2245
    link(l1, mod1, m1, 3994, m2)
    cav cavity m2.p1.o
    modes(maxtem=10)


    fd E_arm m2.p1.i l1.f  # Field detector
    """
)

sol = model.run(
    """
    series(
        eigenmodes(cavity, -mod1.f, name="l9"),
        eigenmodes(cavity, 0,       name="c0"),
        eigenmodes(cavity, mod1.f,  name="u9"),
        xaxis(m2.phi, lin, -10, 190, 500, name="scan")
    )
    """
)

signal = np.sum(abs(sol["scan"]["E_arm"])**2,axis=1)

Now add some simulated noise to make the signal in the data look more realistic. The noise on the power measured by the photodiode is modelled as random draws from a normal distribution with: mean zero and standard deviation \(10^{-2.5}\), and truncated to positive values (because measurements of optical power are non-negative). The mock data is then just the sum of the simulated signal and modelled noise.

mean = 0
variance = 1e-5
lower_bound = 0  # Ensuring no negative values

# Calculate the standard deviation (since variance = std_dev^2)
std_dev = np.sqrt(variance)

# Define the bounds for the standard normal distribution
a = (lower_bound - mean) / std_dev
b = np.inf



# Generate samples from the truncated normal distribution
num_samples = 501  # You can adjust the number of samples as needed
noise_samples = truncnorm.rvs(a, b, loc=mean, scale=std_dev, size=num_samples)

# Add signal to noise to create mock data.
noise = noise_samples
data = signal + noise

Inference on radius of curvature of one of the cavity mirrors

Now assume we do not know the radius of curvature inside the resonant cavity, but a priori, we have an idea of its size, and uncertainty. In this example, we will assume a nominal range of RoC’s between

\[-(1940- 5\times 10^{-1})\,m \leq R \leq -(1940+5\times 10^{-1})\,m. \]

From our noisy data (the cavity scan) we would like to infer what we can about the radius of curvature. For simplicity, we’ll assuming that everything else about the cavity (beam offsets, etc…) are known. Concretely, we will use Bayesian inference to perform this measurement. In this framework, we will compute the probability density of different values of the RoC given our data. This is called the “posterior probability density” of the RoC. We will model our a priori knowledge of the radius of curvature as being uniformly, and randomly distributed within this range. In practice, we would want to determine the distribution of measurements empirically. However, this uniform distribution will serve its purpose in this toy model by allowing us define a prior distribution for inference.

Rc_max = -(1940+5e-1)
Rc_min = -(1940-5e-1)

def prior_transform(params):
    #sample in curvature (inverse RoC)
    curvature_min = 1./Rc_max
    curvature_max = 1./Rc_min
    curvature_unscaled = params

    curvature = uniform_rescale(curvature_min, curvature_max, curvature_unscaled)

    return curvature

def uniform_rescale(minimum, maximum, val):

    return val*(maximum-minimum) + minimum

def lnLikelihood(params):

    Rc = 1./params
    #dx = params[1]
    #dy = params[2]


    try:

        model.m1.Rc = Rc

        _sol = model.run(
            """
            series(
                eigenmodes(cavity, -mod1.f, name="l9"),
                eigenmodes(cavity, 0,       name="c0"),
                eigenmodes(cavity, mod1.f,  name="u9"),
                xaxis(m2.phi, lin, -10, 190, 500, name="scan")
            )
            """
        )


        _signal = np.sum(abs(_sol["scan"]["E_arm"])**2,axis=1)
        res = data - _signal
        if np.any(res<0):
            return -np.inf
        else:
            mean=0
            a_trunc = 0
            a = (a_trunc - mean) / std_dev
            lnL = truncnorm.logpdf(res, a, b, loc=mean, scale=std_dev).sum()
            return lnL

    except Exception as e:
        return -np.inf

nlive = 500      # number of live points
bound = 'multi'   # use MutliNest algorithm for bounds
ndims = 1         # two parameters
sample = 'auto'  # random walk sampling
#sample = 'unif'
tol = 0.1         # the stopping criterion
walks = 50
update_interval = 4.1

import dynesty

sampler = dynesty.NestedSampler(lnLikelihood, prior_transform, ndims,
                   bound=bound, sample=sample, nlive=nlive, update_interval=update_interval)
sampler.run_nested(print_progress=True)

res = sampler.results # get results dictionary from sampler

1057it [16:56:13, 57.69s/it, +500 | bound: 2 | nc: 1 | ncall: 3767 | eff(%): 47.658 | loglstar:   -inf < 2516.248 <    inf | logz: 2512.206 +/-    nan | dlogz:  0.001 >  0.509]

from dynesty.utils import resample_equal

# draw posterior samples
weights = np.exp(res['logwt'] - res['logz'][-1])
posterior_samples = resample_equal(res.samples, weights)

print('Number of posterior samples is {}'.format(len(posterior_samples)))

# plot using corner.py
np.save("posterior_samples", posterior_samples)

Number of posterior samples is 1557

posterior_samples = np.load("posterior_samples.npy")
posterior_samples.T[0] = 1/posterior_samples.T[0]

import corner

corner.corner(posterior_samples,labels=['RoC (m)'],quantiles=(0.95,0.05))
plt.show()
../_images/posterior.svg

The resulting histogram posterior and prior can be plotted using:

plt.hist(np.ravel(posterior_samples),density=True)
plt.hist(np.random.uniform(low=Rc_min, high=Rc_max, size=len(np.ravel(posterior_samples))),
         density=True,alpha=0.5)
plt.xlabel("RoC (m)")
plt.ylabel("probability density")

Text(0, 0.5, 'probability density')
../_images/post_prior.svg