Uncertainty correction

[NPounder]

Corrections to the Hessian can happen after the assimilation has returned a value and has calculated the approximation to the Hessian, broadly speaking after every band has been processed, or before the following line

KaFKA/kafka/linear_kf.py

Line 271 in 81c97a1

P_forecast = P_analysis

The function to calculate the correction term will calculate for a single band the Hessian from the emulator (gp.hessian(x)) which yields H'' (a 4x4 matrix). The other term is just innovation divided by the per-band variance, which is a scalar. The other ingredients are

• The analysis state x_analysis
• The per pixel uncertainty
• The innovations (check that they're defined in the right order eg H(x) -y )

Easiest thing is to use the cached_obs structure, which is a n_bands list where each element contains:

1. observations
2. R_mat (uncertainty matrix) (needed)
3. Mask (needed)
4. Metadata
5. Emulator (needed)

Innovations come from innovations_prime (but you might just want to either change their definition, or just put a - in front)

[jgomezdans]

This is the code that calculates the hessian for a single pixel

def hessian ( gp, x0, hess_approx, C_obs_inv, Hx, y):
    big_hess_H = np.zeros((2, 7, 7))
    for band in xrange(2):
        if band == 0:
            selecta = np.array([0,1,6,2])

        else:
            selecta = np.array([3, 4, 6, 5])
        xp = np.atleast_2d(x0[selecta])
        little_hess = gp.hessian(xp).squeeze()
        for i, ii in enumerate(selecta):
            for j, jj in enumerate(selecta):
                big_hess_H[band, ii, jj] = little_hess[i, j]
    # Big hess is now filled up
    r = Hx.squeeze() - y
    z = C_obs_inv.dot(r)
    corr_term = (big_hess_H*z[:, None, None]).sum(axis=0)
    return hess_approx - corr_term

[jgomezdans]

There is propagation between two bands, but they should be yield the same result as they should be getting the same input
The solver gets the following parameters for band 0 and then band 1:

observations, mask*, H_matrix,
x_forecast, P_forecast, P_forecast_inverse,
R_mat, the_metadata

Only the mask is the same. The analysis from band 0 goes in as x_forecast of band 1. So if x_forecast and/or P_forecast_inverse are different when passed on to the solver (call to self.solver in method assimilate in LinearKalman) between both versions of the code, we need to know why. I can't remember where the Hessian correction was at the end, I thought it ought to be after the the solution for both bands has coverged, but you're right that we are using the posterior of band0 as prior for band1. In my simplified code, I do both bands together, and this is different. The error is (likely) in the sequence of solution in KaFKA. Currently, it does

for observations in this timestep:
     while not converged:
           for band in bands:
                xa, Pa_inv = solve(xf, Pf_inv, band)
                xf = xa
                 Pf_inv = Pa_inv

This is wrong, it should be

for observations in this timestep:
    for band in bands:
             while not converged:
                  xa, Pa_inv = solve(xf, Pf_inv, band)
              xf = xa
              Pf = Pa_inv

Failing that, do all bands together, but this takes up a lot of space and limits doing lots of pixels together

[NPounder]

Currently the hessian correction is applied before the Pf =Pa_inv, it has to be! but I guess this is in the wrong place. In current KaFKA this propagation is the only way the two bands are combined (as far as I can see) the results from band 0 don't persist outside of the loop...