Before starting, understanding the GRIP nested architecture

GRIP is a tool box, designed to be flexible as a requirement to be compatible with any nulling data. It is also simulation-based, by using a model of the instrument (also named as null depth estimator) to reproduce sequence of null depths which statistical distribution is compared to the one of the measurements.

It is structured in three layers: - top: functions calling the fitting algorithm (e.g. likelihood, least squares or MCMC) - middle: functions calling the cost function (e.g. binomial likelihood, $chi^2$, posterior) and the builder of the model (namely grip.create_histogram_model) - bottom: the instrument model (e.g. grip.lbti_model).

The bottom layer is brought by the user and this tutorial explains how to build a custom model in the second part.

From these layers, it is clear that top functions not only to need to pass their own arguments (and keywords) but also those of the middle and the bottom layers. To do so, arguments of the middle and bottom layers are put in a list or equivalent. This list is passed to the middle layer which extracts its own arguments and pass the others to the bottom layer. The challenge is that the bottom layer has an arbitrary number of parameters of any type, solely constrained by the user, and GRIP needs to get an interface to this user-based layer.

In addition, to stay in the trend of the other fitting algorithm and to make a clean interface with third-party packages, only the constants of the layers are embedded into lists while the quantities of the problem (e.g. the histograms, the bins, the parameters to fit, the error bars) are explicit arguments in the different layers.

Let’s have a look of the different layers, starting with the bottom In the tutorial #3, we have fitted LBTI data (bottom layer) with a likelihood as a cost function (middle layer) run by a scipy.optimize.minimize wrapper (top layer).

The code looks like this:

grip.minimize_fit(grip.neg_log_multinomial, grip.create_histogram_model, initial_guess, \
   null_axis, null_pdf, yerr=None, bounds=bounds_fit, diff_step=diffstep_lklh, \
   func_args=cost_fun_args, func_kwargs=cost_fun_kwargs)\

grip.minimize_fit is the wrapper of scipy.optimize.minimize which takes for its own arguments all what is displayed but the last two. In particular, we recognize the cost function grip.neg_log_multinomial and the histogram builder grip.create_histogram_model, both to be called in the middle layer. initial_guess is both a requirement of scipy.optimize.minimize and an argument for the functions in the middle and bottom layers.

func_args=cost_fun_args and func_kwargs=cost_fun_kwargs are arguments useless for grip.minimize_fit but critical for grip.neg_log_multinomial and grip.create_histogram_model.

And cost_fun_args is defined as follows:

cost_fun_args = (wl_scale, grip.lbti_model, instrument_constants, rvu_opd, cdfs_top, rvus_top)

In cost_fun_args, we recognise the arguments needed by grip.neg_log_multinomial. In particular, the bottom layer is part of the tuple:

  • grip.lbti_model: the instrument model, which forms the bottom layer of GRIP

  • instrument_constants: the constants (i.e. non-random values) of the model.

This means that when the middle layer called the bottom layer, it uses these parameters, along with the one generated by the histogram builder.

cost_fun_kwargs is defined as follows:

cost_fun_kwargs = {'use_this_model':out, 'n_samp_per_loop':n_samp_per_loop, 'nloop':nloop, 'verbose':True}

We recognise a keyword for the cost function grip.neg_log_multinomial: use_this_model. And the others are for the histogram builder create_histogram_model(params_to_fit, xbins, wl_scale0, instrument_model, instrument_args, rvu_forfit, cdfs, rvus, **kwargs) so that:

  • params_to_fit=initial_guess=[an, mu_opd, sig_opd]

  • xbins=null_axis

  • wl_scale0=cost_fun_args[0]

  • instrument_model=cost_fun_args[1]=grip.lbti_model

  • instrument_args=cost_fun_args[2]=instrument_constants

  • rvu_forfit=cost_fun_args[3]=rvu_opd

  • cdfs=cost_fun_args[4]=cdfs_top

  • rvus=cost_fun_args[5]=rvus_top

  • kwargs=cost_fun_kwargs[1], cost_fun_kwargs[2] and cost_fun_kwargs[3]

and in the bottom layer, we have lbti_model(na, wavelength, wl_idx, spec_chan_width, phase_bias, opd, IA, IB, thermal_bckg, sigma_eps):

  • na=initial_guess[0]

  • wavelength=cost_fun_args[0]

  • wl_idx= argument created by create_histogram_model

  • spec_chan_width=instrument_constants[0]

  • phase_bias=instrument_constants[1]

  • opd= sequence of random values generated thanks to initial_guess[1]=mu_opd and initial_guess[2]=sig_opd and rvu_forfit=rvu_opd

  • IA, IB, thermal_bckg and sigma_eps are sequence of random values generated thanks to cdfs=cost_fun_args[4]=cdfs_top and rvus=cost_fun_args[5]=rvus_top.

This structure is an answer to this challenging interface issue between user-based functions, third-party packages and data with different structures. If a simpler solution is to be found, it will be implemented as a major update of GRIP and backward compatibility may not be guaranted.