Interpolator: Use more advanced interpolation methods

[Jammy2211]

The interpolator now works for the CTI use case, here is a simple boiled-down example which is a simplified version of (https://github.com/Jammy2211/autocti_workspace_test/blob/main/temporal/model.py):

"""
Temporal: Individual Fits
=========================

In this script, we will fit multiple charge injection imaging to calibrate CTI, where:

 - The CTI model consists of one parallel `TrapInstantCapture` species.
 - The `CCD` volume filling is a simple parameterization with just a `well_fill_power` parameter.
 - The `ImagingCI` is simulated with uniform charge injection lines and no cosmic rays.

The multiple datasets are representative of CTI calibration data taken over the course of a space mission. Therefore,
the model-fitting aims to determine the increase in the density of traps with time.

This script fits each dataset one-by-one and uses the results post-analysis to determine the density evolution
parameters. Other scripts perform more detailed analyses that use more advanced statistical methods to provide a
better estimate.
"""
import numpy as np
from os import path
import autofit as af
import autocti as ac

"""
We now load every image, noise-map and pre-CTI charge injection image as instances of the `ImagingCI` object.

We load and fit each dataset, accquried at different times, one-by-one. We do this in a for loop to avoid loading 
everything into memory.
"""
instance_list = []

time_list = range(0, 10)

for time in time_list:

    dataset_time = f"time_{time}"
    dataset_time_path = path.join(dataset_path, dataset_time)

    imaging_ci_list = [
        ac.ImagingCI.from_fits(
            image_path=path.join(dataset_time_path, f"data_{int(norm)}.fits"),
            noise_map_path=path.join(
                dataset_time_path, f"norm_{int(norm)}", "noise_map.fits"
            ),
        )
        for layout, norm in zip(layout_list, norm_list)
    ]

    """
    __Model__
    
    We now compose our CTI model, which represents the trap species and CCD volume filling behaviour used to fit the 
    charge  injection data. In this example we fit a CTI model with:
    
     - One parallel `TrapInstantCapture`'s which capture electrons during clocking instantly in the parallel direction
     [2 parameters].
    
     - A simple `CCD` volume filling parametrization with fixed notch depth and capacity [1 parameter].
    
    The number of free parameters and therefore the dimensionality of non-linear parameter space is N=3.
    """
    parallel_trap_0 = af.Model(ac.TrapInstantCapture)
    parallel_traps = [parallel_trap_0]
    parallel_ccd = af.Model(ac.CCDPhase)
    parallel_ccd.well_notch_depth = 0.0
    parallel_ccd.full_well_depth = 200000.0

    model = af.Collection(
        cti=af.Model(
            ac.CTI2D,
            parallel_trap_list=[parallel_trap_0],
            parallel_ccd=parallel_ccd,
        ),
        time=time, # <- This is key to the interpolator
    )

    """
    __Search__
    
    The CTI model is fitted to the data using a `NonLinearSearch`. In this example, we use the
    nested sampling algorithm Dynesty (https://dynesty.readthedocs.io/en/latest/).
    """
    search = af.DynestyStatic(
        path_prefix=path.join(dataset_label, dataset_time),
        name="parallel[x1]",
        nlive=50,
    )

    """
    __Analysis__
    
    The `AnalysisImagingCI` object defines the `log_likelihood_function` used by the non-linear search to fit the 
    model to the `ImagingCI` dataset.
    """
    analysis_list = [
        ac.AnalysisImagingCI(dataset=imaging_ci)
        for imaging_ci in imaging_ci_list
    ]
    analysis = sum(analysis_list)
    analysis.n_cores = 1

    """
    __Model-Fit__
    
    We can now begin the model-fit by passing the model and analysis object to the search, which performs a non-linear
    search to find which models fit the data with the highest likelihood.
    """
    result_list = search.fit(model=model, analysis=analysis)

    instance_list.append(result_list[0].instance)

"""
We now use the interpolator to create a new instance of the model, where the trap densities
are at time=1.5.
"""

interpolator = af.LinearInterpolator(instances=instance_list)
instance = interpolator[interpolator.time == 1.5]

print(instance.cti.parallel_trap_list[0].density)

Interpolation Method:

For Euclid, we do not want to assume y = mx + c but instead need a more general interpolation tool.

This would either be some sort of spline fit, or an equation with more parameters.

We also need the interpolation to factor in the inferred errors on each model fit, which I believe we currently do not do.

@rjmassey do you have any suggestions for how to do this, or what algorithm to you?

[rjmassey]

TLDR: don't know until we get the data. But could we have a hierarchical model in which beta has a universal value; ratios of trap densities vary slowly over time (as splines?); and the total trap density is y=mt+c within specific ranges of t, but c jumps between them.

There are two big uncertainties:

1. Either all the data shows a smooth evolution for the entire mission duration, and can be modelled in a single function, or it can only be modelled in chunks. I expect Euclid to be only piecewise smooth, because of e.g. Solar events (but I expected the same for Hubble, and that's been smooth for 2 decades).

Assuming it will be only piecewise smooth, we still have a choice: we could either look at data +/-n days from the date in question - that's essentially a spline fit, as you sugest. However, we will need to choose n - and whatever we choose, that'll give a bad fit near any discontinuity. I would instead suggest a model with lots of (y=mt+c) models, each within a range t1<t<t2, where t1 and t2 are times of discontinuity events found from an initial nonlinear search, i.e. t2 for one period is the same as t1 for the next period. Is that feasible?

2. The other issue is that we expect some parameters to be constant (within measurement uncertainty). For example, the values of beta and the relative density of trap species were constant for Hubble. It is probably be worth treating these differently from the parameters that we expect to grow over time, such as the sum of trap desnities (which can be modelled as the piecewise y=mt+c or otherwise). There could be several ways of dealing with this, but my intuition is that it could be as simple as retaining the covariance between parameters, so we can include them in the fit, to get a single value of beta for each continuous period. The best approach would be to use a hierarchical model in which beta has a universal value; ratios of trap densities vay slowly over time; and the total trap density varies as y=mt+c within specific ranges of t, but c jumps between them.

[rhayes777]

I've added a CubicSpline interpolator here:

#687

From what RIchard said it sounds like it might make sense to instead use a sliding window around the time we're interested in and only for those attributes which we know vary over time. If the relationship tends to be linear locally then the existing LinearInterpolator may be sufficient.

[rjmassey]

Sorry for the delay in getting back to you. Thanks for the CubicSpline.
There are two big items in my wishlist for this time evolution fitter.

1). Please could the likelihood calculation use the 2D covariance matrix between measurements, rather than just the sigma on individual measurements?

• The covariance between parameters is (according to James) already returned by the fit to each epoch. This should be fairly easily incorporated to the time series fitter - the difference is neatly explained in one equation at https://stats.stackexchange.com/questions/369167/prove-positivity-of-chi-squared-statistic-with-general-covariance-matrix and, for a concrete an example, scipy.optimize.curve_fit accepts sigma as either a 1D vector or 2D matrix (the 1D vector is the diagonal elements of the 2D matrix).
• To implement this in pyAutoFit, rather than separately fitting beta(t) and rho(t) and tau(t) etc, we'd first need to simultaneously fit [beta,rho,tau,...](t), using errors [sigma_beta,sigma_rho,sigma_tau,...]. That won't change anything by itself, but it will make it possible to switch the error to the covariance matrix. We'd then need to construct the covariance for multiple epochs by pasting individual epochs' coavriance matrices near the diagonal, with zeros everywhere else. The hardest part will be to invert this matrix (equivalent to 1/sigma), but it is a very sparse matrix.

2). Rather than a sliding window, I would prefer to use a model in which all the data (from a time window) is fitted, but jumps are allowed in certain parameters. Having talked to people from the Gaia mission, I expect most parameters to very smoothly (linear/polynomial/cubicspline), but rho_i parameters to vary linearly (or even be constant) between jumps at discrete times t_1, t_2, t_3,... t_n.

• Would it be possible to have a fitting function like that, which uses a lot of data, but has these extra parameters t_1, t_2, t_3,... t_n and rho_i_0, rho_i_1, rho_i_2,... rho_i_n (where perhaps n is specified in advance, or maybe even the t's are specified in advance)?

[Jammy2211]

The covariance matrix is available here:

https://github.com/rhayes777/PyAutoFit/blob/main/autofit/non_linear/samples/pdf.py

    def covariance_matrix(self) -> np.ndarray:
        """
        Compute the covariance matrix of the non-linear search samples, using the method `np.cov()` which is described
        at the following link:

        https://numpy.org/doc/stable/reference/generated/numpy.cov.html

        Follow that link for a description of what the covariance matrix is.

        Returns
        -------
        ndarray
            A covariance matrix of shape [total_parameters, total_parameters] for the model parameters of the
            non-linear search.
        """
        return np.cov(m=self.parameter_lists, rowvar=False, aweights=self.weight_list)

for a concrete an example, scipy.optimize.curve_fit accepts sigma as either a 1D vector or 2D matrix (the 1D vector is the diagonal elements of the 2D matrix).

Yay.

[rjmassey]

After a conversation with James, a sliding time window is definitely not a great idea - we'd use the calibration process many times and end up with lots of fits in overlapping windows. When it came time to use a best-fit parameter, we would be unsure which one to use. Much cleaner to future us if we just generate a single model - the best-fit parameters of which we upate occasionally, but which can always be used in its most recent incarnation.

So an ability to fit only data within a time window is great functionality if it already exists, but I imagine that time window will probably always be "forever".

[rhayes777]

Could be worth having a call about this if you have time

[rjmassey]

Yep, can talk now if you like

[rhayes777]

Could you email a zoom link?

[Jammy2211]

Going to close this for now, as we are not actively working on it, but I am sure CTI will motivate us to come back to this inthe future.