ivalaginja
commented
5 years ago This whole thing about getting the errors from the correct position in the Hessian is kind of nasty, so here we go step by step:
There are 22 parameters that make a systematic model. The difference between the 50 systematic models of WFC3 is, which of these 22 are included in the model fit and which are not. All of these 22 parameters used to be grouped together in an array called p0 and it looked like this:
p0 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21 ]
p0 = np.array([rl, flux0, epoch, inclin, MsMpR, ecc, omega, Per, T0, c1, c2, c3, c4, m_fac, HSTP1, HSTP2, HSTP3, HSTP4, xshift1, xshift2, xshift3, xshift4])So the barebones systematic model has all of these parameters in it. The way in which we go through all the 50 different models is that for each model, we construct a single vector that contains a 1 in all the parameters spots that are fixed (frozen) and a 0 for all the parameters that are getting fit (thawed). If we stack all of these vectors together, we get a 22 x 50 grid that we call our systematic grid.
During the fitting that happens in our code, we feed the fitting routine each 22-parameter vector one at a time, and we get results per systematic model. While a big part of the systematic grid is unique to the instrument that produced our data (WFC3 currently), other parts of the grid are variable because they are not inherent to the instrument, but to the observed object. And while we want to go through all the 50 parametrized models that are inherent to WFC3, there are six parameters out of the 22 that we just set such that we decide which of the six are frozen and which are thawed, purely based on the scientific result we need.
All of this is captured in the function wfc3_systematic_model_grid_selection() in the module margmodule.py. Let's illustrate the previous paragraph by looking at p0 above:
1) Parametrized systematic models:
m_fac, HSTP1, HSTP2, HSTP3, HSTP4, xshift1, xshift2, xshift3, xshift4
These 9 parameters, with indices 13-21 (so, the 9 very last ones) define the 50 different systematic models described in Wakeford et al. 2016. While they are different from model to model, their possible combinations are static and never get changed.
2) Input constants:
omega, Per, T0, c1, c2, c3, c4
These parameters with indices 6-12 (right before the parametrized systematic model parameters in 1. ) are constant input parameters and will always have a 1 in their spot.
3) Permanent free parameters:
rl, flux0
These two are at positions 0-1 of the grid and will always carry a 0, because we always want to fit for them. Note: flux0 is the first data point of our y-array. As such, we don't really care about it in the marginalisation, nor do we care about its error, especially since the fit directly returns an error array on the entire flux array, as opposed to just its first point. This means that we never touch the error on flux0 that we would get from the Hessian.
4) Variable parameters:
epoch, inclin, MsMpR, ecc
At indices 2-5, we have to make a decision before we run the code about which of these four are thawed and which are frozen. To implement this in the easiest possible way, we confined the combinations for these four parameters to six choices, one of which we have to declare in our configfile with the variable grid_selection before running the code:
-- fix_time: 1 1 1 1
-- fit_time: 0 1 1 1
-- fit_inclin: 1 0 1 1
-- fit_msmpr: 1 1 0 1
-- fit_ecc: 1 1 1 0
-- fit_all: 0 0 0 1
These do not capture all possible permutations for the four, but that is fine, since they capture all the possible cases we could sensibly be interested in.
The Hessian matrix of the fit we run returns the errors of all thawed parameters. This means that it contains the errors of all parameters that were marked by a 0, and only those. This means that the error on one particular parameter will appear at a different index, depending on what we're currently fitting. This makes it really hard to extract the correct error from the Hessian.
However, thankfully we are not interested in all the errors that the fit calculates. Particularly, we only need the errors on the first six parameters of the p0 array (indices 0-5), which correspond to the variables rl, flux0, epoch, incl, msmpr and ecc. If we take an even closer look, we know that rl and flux0 are always thawed anyway, so we can rely on the fact that the error on rl will always come back on index 0 of the Hessian and the error on flux0 will always come back on index 1 of the Hessian. This leaves us with only having to deal with the correct ordering of the errors on epoch, incl, msmpr and ecc. And since we have the six options mentioned above for our grid_selection, which tells us which of these four parameters even return an error (because they are thawed), this is what we have to deal with. For each grid_selection, the according errors of these four parameters will be in a different spot in the Hessian.
Still thinking about how to go about this. The most straightforward solution is to implement what feels like a thousand if statements, but I am trying to come up with something slightly smarter than that.
As a side note:
The bright side is, because we switched to actively tiling and stacking the individual parts of the systematic grid together as opposed to purely hard coding it, the grid segment that defines the parametrized systematic models can easily be swapped out for a different parametrized grid. I.e., it should be very straightforward to implement this for a new instrument, e.g. STIS, since all parameters up to including index 12 in the p0 array will be treated exactly the same like they are now.
As outlined in issue #24, we want to discard the "Covariance" option for the error option and make sure that the correct errors get read from the Hessian, no matter which parameters are frozen and which are thawed, which we take care of in this PR.
Getting the errors from the Hessian in the correct order will be tricky, since it only returns error values on thawed parameters. So, getting the correct errors will depend on the instrument (we have one systematic model grid per instrument; we're currently only working with WFC3 data, so that part of the grid stays the same) and on the fit selection, which tells us which of the variable parameters are thawed and which are fixed.
In this PR, I will focus on proper flexibility of getting the errors with the fit selection only and will leave adjustments to other instruments for later.