J. Kissel

Picking up where I last left off, in LHO:69898, I had divided up the ten 4-hour segments of IFO thermalization times during ER15 into slices of average X&Y arm cavity power, or into power bins of 1.0 W. This aLOG covers the attempts to create a representative model of the collection of transfer functions during each power bin using methodology developed during O3 -- fitting the data using a Gaussian Process Regression with the same radial basis function kernel we use to fit for other unknown frequency dependence in the sensing and actuation functions.

Executive Summary: Gaussian Process Regression's cons continue outweigh its pros as a transfer function fitter with the kernel we typically use. I'm able to finagle the hyper parameters of the kernel model such that the regression output tracks the median behavior of the transfer function, but it does *not* reflect the large, measured, variability in the transfer function (which we traditionally try to cover with the Gaussian uncertainty). This may have to be "good enough," but I'm going to have to farm this out to other specialists if we need better.

Step 1: Reorganize the *real* transfer function data, and add uncertainty
Attachment 1: sensingFunction_syserror_vs_powerslices_originaldata_vs_binneddata.pdf
    As of May 24 2023, and LHO:69898, I had found *a* way to group the data into power slices and sort the 3D data arrays. However, in doing so, I had NOT saved any of the coherence-based uncertainty on the measurement points. So, as of this aLOG, I have now added the measurement, coherence-based uncertainty to each of the data points in all the collections of data. While doing so, I had some more time with 3-D python array manipulation, so I reorganized things as well. 

This collection of plots, shows all the collection of power bin's transfer functions and their associated power, now with quantified uncertainty. On the left, the bode plot of the transfer function organized in one way, and on the right, the same data is just organized differently -- but serves as a sanity check that my new organization preserves all the right data and uncertainty assignment.

Step 2: Gather some statistical intuition on the data by finding 68th, 95th, 99.7th quantiles as one way to quantify the uncertainty of the parent distribution at each frequency point. 
Attachment 2: sensingFunction_syserror_vs_powerslices_binneddata_freqhist.pdf
    One of the fundamental assumptions of a Gaussian Process regression is that the family of curves that you're given it is sampled from a -- you guessed it -- Gaussian parent distribution of curves. So... is this collection of data Gaussian?

This plot shows histograms of each frequency point's magnitude and phase within a power bin. Overplotted on top of the histogram is 
   - the (red) median, or 50th percentile, then 
   - the (darkest pink) 0.1587th,0.8413th "1-sigma" quantiles, 
   - the (medium pink) 0.9772nd,0.0228th "2-sigma" quantiles, and 
   - the (light pink) 0.9986th,0.0013rd "3-sigma" quantiles.
Of course, each bin's worth of data is diverse in number of transfer function values. Some bins only have 3 traces, and others have 170. As such, I've let python's histogram function determine the number of bins by setting histogram_bin_edges='auto'. But still, this display gives *me* a good impression of the data, and the first impressions that a Guassian Process regression is just not the right tool:
   - the distributions of transfer function values at each frequency are *not* Gaussian, and
   - the "uncertainty" of the distribution, quantified by the 1-, 2-, and 3- sigma quantiles is not symmetric. 

Alas, we proceed.

Step 3: Add "fake" data the the transfer function collections, convert from mag/phase, to complex, then to real/imaginary parts, then flatten in prep for input into GPR.
Attachment 3: sensingFunction_syserror_vs_powerslices_binneddata_vs_gprinputdata.pdf
    As we've done in O3, we want a model of these transfer functions to serve as a systematic error (with associated uncertainty) that can be "stapled" in to the overall calibration's modeled response function systematic error (with associated uncertainty). In doing so, we want this model to add systematic error only where where think the response function is changing as a function of thermalization. However, Guassian Process Regression gives you large uncertainty where there are no measured frequency points to "anchor" the fit. Further, due to the querks of the skikit learn's python implementation of GPR, (or at least the way that Craig Cahillane and then Evan Goetz figured out how to get it to "work"), you must flatten the collection of input data on to a single vector (i.e. the nTFs by mFreq complex 2D array must be reshaped into a (nTFs * mFreq) x 1 vector).

This collection of plots shows the sanity check that the addition of fake data points, and the flattening works.

The fake data points have a transfer function value of 1.0 + 0.0j i.e. a magnitude of 1.0 [(ct/m) / (ct/m)], and a phase of 0, the sensing function residual [meas/model] if model perfectly agrees with measurement. I've assigned these fake data points an uncertainty of 1e-4, or 0.01%. The *number* of fake data points required depends on the radial basis function hyper parameter "length scale"

Step 4: Explore the "standard" O3-era "radial basis function"-inclusive hyper parameter space in order to create a "satisfying" fit.
Attachment 4: sensingFunction_syserror_vs_powerslices_gprfit_lengthscale_0p25_nxtrapnts30.pdf
    The GPR transfer function kernel we continue use to describe unknonwn frequency dependence in "our" transfer function residuals since Evan added additional terms after Craig's O1/O2 work (see "Slides for follow up on action item" from G1901479) is as follows:

    F * RBF(\ell) + C

where
    F = a "ConstantKernel" function that serves to be what I call an "amount of frequency dependence" coefficient that determines how much of the radial basis function is "used"

    RBF(\ell) = a "Radial Basis" function, which is a Gaussian function of the dimensionless hyper-parameter \ell, where the function computes the correlation of neighboring data points with respect to the length scale, \ell. As we've shoved the transfer function data into this kernel function, it's computing the Euclidean distance between the neighboring, logarithmic frequency points of the real and imaginary points of the complex transfer function.
        Working through a bit of algebra, as in G2101319, the physical meaning of the length scale becomes

            delta f = f * (10^\ell - 1)

        i.e. the RBF looks for lower frequency points have tighter (in frequency) correlation.

    C = another "ConstantKernel," but in this case used as an additive constant. This is used as a representative of what an "ideal" residual should look like, i.e. a transfer function of 1.0 + 0.0j at all frequencies. 

As I've been reminded through my exploration, if you start probing for length scales, \ells, that are smaller than the log(f) frequency separation, the GPR will return a large uncertainty on the transfer function in between your data points. Thus, the number of *fake* data points needs to also start to come into play if your length scale gets small.

So, the parameter space to explore is the priors on the constant for the two ConstantKernels as well as the prior on the RBF length scale \ell.
    In the end, I used
        F = value of 0.5, with uniform probability between (0.4, 1.0). 
            I interpret this as allowing for the frequency dependence to be "all the way on" (i.e. 1.0) or a little-bit 
            less than "half" on (0.4).

        \ell = value of 0.25, with uniform probability between (0.125, 0.375). 
            This value I determined entirely empirically, by trying values between [0.125, 0.25, 0.5, 0.75, 1.0, 2.0]

        C = value of 1.0, with uniform probability between (0.9, 1.1).

        nExtraPoints = 30. 
            One can safely run with 15 points on length scales down to 0.5, but below 0.5 the length scales become comparable 
            to the frequency separation in these points, so the frequency space *between* the extra points starts to become 
            prohibitively uncertain, which spoils the "frequencies above the thermalization region shall not impact the greater 
            systematic error budget" rule.

This plot shows the Gaussian Process Regression fit against the Median and 1-, 2-, and 3- sigma quantiles of the 4 *real* data frequency points, as well as the original data.
From here, we arrive at our executive summary: I'm able to finagle the hyper parameters of the kernel model such that the regression output tracks the median behavior of the transfer function, but it does *not* reflect the large, measured, variability in the transfer function (which we traditionally try to cover with the Gaussian uncertainty) -- i.e. some of the way there, and may be good enough. For the fits that are the worst, the lowest power bins, I think this is an OK compromise if we recall that the IFO stays in these lowest power bins for the *least* amount of time given the "exponential" nature of the thermalization, and sometimes the IFO is past these arm powers before we even hit "NOMINAL_LOW_NOISE."

In fact, over my years of asking for One Transfer Function Fitter to Rule Them All, I've found that this exact problem is the hardest part about transfer function fitting in LIGO -- finding a fit or model that appropriately describes the *uncertainty*, or in our application here, the *variability* in the transfer function. I guess it should be a surprise that I wasn't able to figure out an entire field of computational metrology (or loosely, parameter estimation) in 1.5 weeks.

So, at this point, 
    (a) We need to install *something* in the overall systematic error budget, and this is a good start,
    (b) There's talk of reducing the power and reverting the SRCL offset, perhaps soon rendering this data set "no longer representative,"
    (c) I need the help of others to do better in terms of capturing the variability (and maybe that means exploring different kernels, or switching to zpk-based fitting like IIRrational or Ethan's ZPK Fitting), 
    (d) We could also try to do all of this same process for the *response function,* rather than the sensing function (which *may* be quicker than 2 months worth of work because I already have examples of how to manipulate the data which was a huge time suck for me)
    (e) I need to move on with my life and make progress elsewhere

So, I hand it off to my fellow calibrators, transfer function fitters, and the world.