Correction term for Hessian uncertainty

kafka/kf_tools.py

@@ -11,6 +11,39 @@
 from utils import block_diag
 
 
+def band_selecta(band):
+    if band == 0:
+        return np.array([0, 1, 6, 2])
+    else:
+        return np.array([3, 4, 6, 5])
+
+
+def hessian_correction_pixel(gp, x0, C_obs_inv, innovation, band, nparams):
+    selecta = band_selecta(band)
+    ddH = gp.hessian(np.atleast_2d(x0[selecta]))
+    big_ddH = np.zeros((nparams,nparams))
+    for i, ii in enumerate(selecta):
+        for j, jj in enumerate(selecta):
+            big_ddH[ii,jj] = ddH.squeeze()[i,j]        
+    big_Hessian_corr = big_ddH*C_obs_inv*innovation
+    return big_Hessian_corr
+
+
+def hessian_correction(gp, x0, C_obs_inv, innovation, mask, state_mask, band, nparams):
+    C_obs_inv = C_obs_inv.diagonal()[state_mask.flatten()]
+    mask = mask[state_mask].flatten()
+    little_hess = []
+    for i, (innov, C, m) in enumerate(zip(innovation, C_obs_inv, mask)):
+        x0_pixel =  x0.squeeze()[nparams*i:nparams*(i+1)]
+        if not m:
+            hessian_corr = np.zeros((nparams,nparams))
+        else:
+            hessian_corr = m * hessian_correction_pixel(gp, x0_pixel,C, innov, band, nparams)
+        little_hess.append(hessian_corr)
+    hessian_corr = block_diag(little_hess)
+    return hessian_corr
+
+
 def tip_prior():
     """The JRC-TIP prior in a convenient function which is fun for the whole
     family. Note that the effective LAI is here defined in transformed space
@@ -19,7 +52,7 @@ def tip_prior():
     Returns
     -------
     The mean prior vector, covariance and inverse covariance matrices."""
-    sigma = np.array([0.12, 0.7, 0.0959, 0.15, 1.5, 0.2, 0.5]) # broadly TLAI 0->7 for 1sigma
+    sigma = np.array([0.12, 0.7, 0.0959, 0.15, 1.5, 0.2, 0.35]) # broadly TLAI 0->7 for 1sigma
     x0 = np.array([0.17, 1.0, 0.1, 0.7, 2.0, 0.18, np.exp(-0.5*1.5)])
     # The individual covariance matrix
     little_p = np.diag ( sigma**2).astype(np.float32)
@@ -159,7 +192,7 @@ def propagate_information_filter_LAI(x_analysis, P_analysis,
     lai_post_cov = P_analysis_inverse.diagonal()
     c_inv_prior_mat = []
     for n in xrange(n_pixels):
-        c_inv_prior[6,6] =  lai_post_cov[n]
+        c_inv_prior[6,6] =  40 #lai_post_cov[n]
         c_inv_prior_mat.append(c_inv_prior)
 
     P_forecast_inverse=block_diag(c_inv_prior_mat, dtype=np.float32)

kafka/linear_kf.py

@@ -35,7 +35,7 @@
 from utils import create_linear_observation_operator
 from utils import create_nonlinear_observation_operator
 from utils import iterate_time_grid
-from kf_tools import propagate_information_filter
+from kf_tools import hessian_correction, propagate_information_filter_SLOW
 
 # Set up logging
 
@@ -58,7 +58,7 @@ class LinearKalman (object):
     goal of this class is not to consider complex, time evolving models, but
     rather grotty "0-th" order models!"""
     def __init__(self, observations, output, state_mask,
-                 state_propagation=propagate_information_filter,
+                 state_propagation=propagate_information_filter_SLOW,
                  linear=True, n_params=1, diagnostics=True,
                  bands_per_observation=1):
         """The class creator takes (i) an observations object, (ii) an output
@@ -267,9 +267,16 @@ def assimilate(self, locate_times, x_forecast, P_forecast,
                             x_forecast, P_forecast, P_forecast_inverse,
                             R_mat, the_metadata)
 
+
+                    P_correction = hessian_correction(
+                                the_emulator, x_analysis, R_mat, 
+                                innovations_prime, mask, self.state_mask, band, self.n_params)
+
+                    P_analysis_inverse = P_analysis_inverse - P_correction
+
                     x_forecast = x_analysis*1
                     P_forecast = P_analysis
-                    P_forecast_inverse = P_analysis_inverse
+                    P_forecast_inverse = P_analysis_inverse*1
 
                 if iter_obs_op:
                     # this should be an option...

kafka/observations.py

@@ -51,7 +51,7 @@
 
 from BRDF_descriptors import RetrieveBRDFDescriptors
 
-from kernels import Kernels
+#from kernels import Kernels
 
 import scipy.sparse as sp
 from scipy.ndimage import zoom
@@ -83,7 +83,7 @@ def get_modis_dates(fnames):
 # TODO needs class for MODIS L1b product too
 # These classes should define emulators
 
-
+'''
 class MOD09_ObservationsKernels(object):
     """A generic M*D09 data reader"""
     def __init__(self, dates, filenames):
@@ -143,7 +143,7 @@ def get_band_data(self, the_date, band_no):
         data_object = MOD09_data(refl, mask, uncertainty, K, sza, vza, raa)
 
         return data_object
-
+'''
 
 class SynergyKernels(object):
     """An object to store, process and update linear kernel weights datasets

kafka/solvers.py

@@ -24,7 +24,7 @@
 import scipy.sparse as sp
 import matplotlib.pyplot as plt
 
-from utils import  matrix_squeeze, spsolve2, reconstruct_array
+from utils import spsolve2
 
 # Set up logging
 import logging

tests/test_kf.py

@@ -13,7 +13,7 @@
 sys.path.insert(0, myPath + '/../')
 
 from kafka.kf_tools import propagate_standard_kalman
-from kafka.kf_tools import propagate_information_filter
+from kafka.kf_tools import propagate_information_filter_SLOW
 
 
 def test_propagate_standard_kalman():
@@ -39,7 +39,7 @@ def test_propagate_information_filter():
     Pi = np.linalg.inv(Pd)
     Q_matrix = np.eye(7)*0.1
 
-    x_forecast, P_forecast, P_forecast_inverse = propagate_information_filter(
+    x_forecast, P_forecast, P_forecast_inverse = propagate_information_filter_SLOW(
         x_analysis, None, Pi, M_matrix, Q_matrix)
     assert np.allclose(
         np.array(P_forecast_inverse.todense()).squeeze().diagonal(),