More generic hessian correction allows for different RT models

kafka/inference/broadbandSAIL_tools.py

@@ -0,0 +1,97 @@
+import logging
+import gdal
+import numpy as np
+import os
+
+from .utils import block_diag
+from .kf_tools import propagate_single_parameter
+
+
+def tip_prior_values():
+    """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
+    where TLAI = exp(-0.5*LAIe).
+
+    Returns
+    -------
+    The mean prior vector, covariance and inverse covariance matrices."""
+    # broadly TLAI 0->7 for 1sigma
+    sigma = np.array([0.12, 0.7, 0.0959, 0.15, 1.5, 0.2, 0.5])
+    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)
+    little_p[5, 2] = 0.8862*0.0959*0.2
+    little_p[2, 5] = 0.8862*0.0959*0.2
+
+    inv_p = np.linalg.inv(little_p)
+    return x0, little_p, inv_p
+
+
+class JRCPrior(object):
+    """Dummpy 2.7/3.6 prior class following the same interface as 3.6 only
+    version."""
+
+    def __init__(self, parameter_list, state_mask):
+        """It makes sense to have the list of parameters and state mask
+        defined at this point, as they won't change during processing."""
+        self.mean, self.covar, self.inv_covar = self._tip_prior()
+        self.parameter_list = parameter_list
+        if isinstance(state_mask, (np.ndarray, np.generic) ):
+            self.state_mask = state_mask
+        else:
+            self.state_mask = self._read_mask(state_mask)
+
+    def _read_mask(self, fname):
+        """Tries to read the mask as a GDAL dataset"""
+        if not os.path.exists(fname):
+            raise IOError("State mask is neither an array or a file that exists!")
+        g = gdal.Open(fname)
+        if g is None:
+            raise IOError("{:s} can't be opened with GDAL!".format(fname))
+        mask = g.ReadAsArray()
+        return mask
+
+    def _tip_prior(self):
+        """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
+        where TLAI = exp(-0.5*LAIe).
+
+        Returns
+        -------
+        The mean prior vector, covariance and inverse covariance matrices."""
+        mean, c_prior, c_inv_prior = tip_prior_values()
+        self.mean = mean
+        self.covar = c_prior
+        self.inv_covar = c_inv_prior
+        return mean, c_prior, c_inv_prior
+
+    def process_prior ( self, time, inv_cov=True):
+        # Presumably, self._inference_prior has some method to retrieve
+        # a bunch of files for a given date...
+        n_pixels = self.state_mask.sum()
+        x0 = np.array([self.mean for i in range(n_pixels)]).flatten()
+        if inv_cov:
+            inv_covar_list = [self.inv_covar for m in range(n_pixels)]
+            inv_covar = block_diag(inv_covar_list, dtype=np.float32)
+            return x0, inv_covar
+        else:
+            covar_list = [self.covar for m in range(n_pixels)]
+            covar = block_diag(covar_list, dtype=np.float32)
+            return x0, covar
+
+
+def propagate_LAI_broadbandSAIL(x_analysis, P_analysis,
+                                     P_analysis_inverse,
+                                     M_matrix, Q_matrix,
+                                     prior=None, state_propagator=None, date=None):
+    """Propagate a single parameter and
+      set the rest of the parameter propagations to the prior.
+      This should be used with a prior for the remaining parameters"""
+    nparameters = 7
+    lai_position = 6
+    x_prior, c_prior, c_inv_prior = tip_prior_values()
+    return propagate_single_parameter(x_analysis, P_analysis,
+                                      P_analysis_inverse,
+                                      M_matrix, Q_matrix,
+                                      nparameters, lai_position,
+                                      x_prior, c_inv_prior)

kafka/inference/kf_tools.py

@@ -16,59 +16,47 @@ class NoHessianMethod(Exception):
     def __init__(self, message):
         self.message = message
 
-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_pixel(hessian, C_obs_inv, innovation):
+    hessian_corr = hessian*C_obs_inv*innovation
+    return hessian_corr
 
 
-def hessian_correction(gp, x0, R_mat, innovation, mask, state_mask, band,
+def hessian_correction(hessian, R_mat, innovation, mask, state_mask,
                        nparams):
     """Calculates higher order Hessian correction for the likelihood term.
     Needs the GP, the Observational uncertainty, the mask...."""
-    if not hasattr(gp, "hessian"):
+    if hessian is None:
         # The observation operator does not provide a Hessian method. We just
         # return 0, meaning no Hessian correction.
         return 0.
     C_obs_inv = R_mat.diagonal()[state_mask.flatten()]
     mask = mask[state_mask].flatten()
+
     little_hess = []
-    for i, (innov, C, m) in enumerate(zip(innovation, C_obs_inv, mask)):
+    for i, (innov, C, m, hess) in enumerate(zip(innovation, C_obs_inv, mask, hessian)):
         if not m:
             # Pixel is masked
             hessian_corr = np.zeros((nparams, nparams))
         else:
-            # Get state for current pixel
-            x0_pixel = x0.squeeze()[(nparams*i):(nparams*(i + 1))]
             # Calculate the Hessian correction for this pixel
-            hessian_corr = m * hessian_correction_pixel(gp, x0_pixel, C,
-                                                        innov, band, nparams)
+            hessian_corr = m * hessian_correction_pixel(hess, C, innov)
         little_hess.append(hessian_corr)
-    hessian_corr = block_diag(little_hess)
+
+    hessian_corr = block_diag(hessian)
     return hessian_corr
 
 
-def hessian_correction_multiband(gp, x0, R_mats, innovations, masks, state_mask, n_bands,
-                       nparams):
+def hessian_correction_multiband(hessians, R_mats, innovations,
+                                 masks, state_mask, n_bands, nparams):
     """ Non linear correction for the Hessian of the cost function. This handles
     multiple bands. """
     little_hess_cor = []
-    for R, innovation, mask, band in zip(R_mats, innovations, masks, range(n_bands)):
-        little_hess_cor.append(hessian_correction(gp, x0, R, innovation, mask, state_mask, band,
-                       nparams))
-    hessian_corr = sum(little_hess_cor) #block_diag(little_hess_cor)
+    for R, hessian, innovation, mask in zip(
+            R_mats, hessians, innovations, masks):
+        little_hess_cor.append(hessian_correction(hessian, R, innovation,
+                                                  mask, state_mask, nparams))
+    hessian_corr = sum(little_hess_cor)
     return hessian_corr
 
 
@@ -96,43 +84,6 @@ def blend_prior(prior_mean, prior_cov_inverse, x_forecast, P_forecast_inverse):
     return x_combined, combined_cov_inv
 
 
-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
-    where TLAI = exp(-0.5*LAIe).
-
-    Returns
-    -------
-    The mean prior vector, covariance and inverse covariance matrices."""
-    # broadly TLAI 0->7 for 1sigma
-    sigma = np.array([0.12, 0.7, 0.0959, 0.15, 1.5, 0.2, 0.5])
-    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)
-    little_p[5, 2] = 0.8862*0.0959*0.2
-    little_p[2, 5] = 0.8862*0.0959*0.2
-
-    inv_p = np.linalg.inv(little_p)
-    return x0, little_p, inv_p
-
-def tip_prior_noLAI(prior):
-    n_pixels = prior['n_pixels']
-    mean, prior_cov_inverse = tip_prior(prior)
-
-
-def tip_prior_full(prior):
-    # This is yet to be properly defined. For now it will create the TIP prior and
-    # prior just contains the size of the array - this function will be replaced with
-    # the real code when we know what the priors look like.
-    x_prior, c_prior, c_inv_prior = tip_prior()
-    n_pixels = prior['n_pixels']
-    mean = np.array([x_prior for i in range(n_pixels)]).flatten()
-    c_inv_prior_mat = [c_inv_prior for n in range(n_pixels)]
-    prior_cov_inverse=block_diag(c_inv_prior_mat, dtype=np.float32)
-
-    return mean, prior_cov_inverse
-
-
 def propagate_and_blend_prior(x_analysis, P_analysis, P_analysis_inverse,
                               M_matrix, Q_matrix, 
                               prior=None, state_propagator=None, date=None):
@@ -240,10 +191,11 @@ def propagate_information_filter_SLOW(x_analysis, P_analysis, P_analysis_inverse
     S= P_analysis_inverse.dot(Q_matrix)
     A = (sp.eye(n) + S).tocsc()
     P_forecast_inverse = spl.spsolve(A, P_analysis_inverse)
-    logging.info("DOne with propagation")
+    logging.info("Done with propagation")
 
     return x_forecast, None, P_forecast_inverse
 
+
 def propagate_information_filter_approx_SLOW(x_analysis, P_analysis, P_analysis_inverse,
                                  M_matrix, Q_matrix,
                                       prior=None, state_propagator=None, date=None):
@@ -289,35 +241,40 @@ def propagate_information_filter_approx_SLOW(x_analysis, P_analysis, P_analysis_
     return x_forecast, None, P_forecast_inverse
 
 
-def propagate_information_filter_LAI(x_analysis, P_analysis,
-                                     P_analysis_inverse,
-                                     M_matrix, Q_matrix,
-                                     prior=None, state_propagator=None, date=None):
-
+def propagate_single_parameter(x_analysis, P_analysis, P_analysis_inverse,
+                               M_matrix, Q_matrix, n_param, location,
+                               x_prior, c_inv_prior):
+    """ Propagate a single parameter and
+     set the rest of the parameter propagations to the prior.
+     This should be used with a prior for the remaining parameters"""
     x_forecast = M_matrix.dot(x_analysis)
-    x_prior, c_prior, c_inv_prior = tip_prior()
-    n_pixels = len(x_analysis)//7
-    x0 = np.array([x_prior for i in range(n_pixels)]).flatten()
-    x0[6::7] = x_forecast[6::7] # Update LAI
-    lai_post_cov = P_analysis_inverse.diagonal()[6::7]
-    lai_Q = Q_matrix.diagonal()[6::7]
+    n_pixels = len(x_analysis)//n_param
+    x0 = np.tile(x_prior, n_pixels)
+    x0[location::n_param] = x_forecast[location::n_param]  # Update LAI
+    lai_post_cov = P_analysis_inverse.diagonal()[location::n_param]
+    lai_Q = Q_matrix.diagonal()[location::n_param]
 
     c_inv_prior_mat = []
-    for n in range(n_pixels):
+    for cov, Q in zip(lai_post_cov, lai_Q):
         # inflate uncertainty
-        lai_inv_cov = 1.0/((1.0/lai_post_cov[n])+lai_Q[n])
+        lai_inv_cov = 1.0/((1.0/cov)+Q)
         little_P_forecast_inverse = c_inv_prior.copy()
-        little_P_forecast_inverse[6, 6] = lai_inv_cov
+        little_P_forecast_inverse[location::location] = lai_inv_cov
         c_inv_prior_mat.append(little_P_forecast_inverse)
     P_forecast_inverse=block_diag(c_inv_prior_mat, dtype=np.float32)
-
     return x0, None, P_forecast_inverse
 
+
 def no_propagation(x_analysis, P_analysis,
                    P_analysis_inverse,
                    M_matrix, Q_matrix,
                    prior=None, state_propagator=None, date=None):
-    """No propagation. In this case, we return the original prior. As the
+    """
+    THIS PROPAGATOR SHOULD NOT BE USED ANY MORE. It is better to set
+    the state_propagator to None and to use the Prior exlicitly.
+
+    THIS IS ONLY SUITABLE FOR BROADBAND SAIL uses TIP prior
+    No propagation. In this case, we return the original prior. As the
     information filter behaviour is the standard behaviour in KaFKA, we
     only return the inverse covariance matrix. **NOTE** the input parameters
     are there to comply with the API, but are **UNUSED**.

kafka/inference/narrowbandSAIL_tools.py

@@ -0,0 +1,82 @@
+import logging
+import gdal
+import numpy as np
+import os
+
+from .utils import block_diag
+from .kf_tools import propagate_single_parameter
+
+def sail_prior_values():
+    """
+    :returns
+    -------
+    The mean prior vector, covariance and inverse covariance matrices."""
+    #parameter_list = ['n', 'cab', 'car', 'cbrown', 'cw', 'cm',
+    #                 'lai', 'ala', 'bsoil', 'psoil']
+    mean = np.array([2.1, np.exp(-60. / 100.),
+                     np.exp(-7.0 / 100.), 0.1,
+                     np.exp(-50 * 0.0176), np.exp(-100. * 0.002),
+                     np.exp(-4. / 2.), 70. / 90., 0.5, 0.9])
+    sigma = np.array([0.01, 0.2,
+                              0.01, 0.05,
+                              0.01, 0.01,
+                              0.50, 0.1, 0.1, 0.1])
+
+    covar = np.diag(sigma ** 2).astype(np.float32)
+    inv_covar = np.diag(1. / sigma ** 2).astype(np.float32)
+    return mean, covar, inv_covar
+
+class SAILPrior(object):
+    def __init__(self, parameter_list, state_mask):
+        self.parameter_list = parameter_list
+        if isinstance(state_mask, (np.ndarray, np.generic)):
+            self.state_mask = state_mask
+        else:
+            self.state_mask = self._read_mask(state_mask)
+            mean, c_prior, c_inv_prior = sail_prior_values()
+            self.mean = mean
+            self.covar = c_prior
+            self.inv_covar = c_inv_prior
+
+
+    def _read_mask(self, fname):
+        """Tries to read the mask as a GDAL dataset"""
+        if not os.path.exists(fname):
+            raise IOError("State mask is neither an array or a file that exists!")
+        g = gdal.Open(fname)
+        if g is None:
+            raise IOError("{:s} can't be opened with GDAL!".format(fname))
+        mask = g.ReadAsArray()
+        return mask
+
+    def process_prior ( self, time, inv_cov=True):
+        # Presumably, self._inference_prior has some method to retrieve
+        # a bunch of files for a given date...
+        n_pixels = self.state_mask.sum()
+        x0 = np.array([self.mean for i in range(n_pixels)]).flatten()
+        if inv_cov:
+            inv_covar_list = [self.inv_covar for m in range(n_pixels)]
+            inv_covar = block_diag(inv_covar_list, dtype=np.float32)
+            return x0, inv_covar
+        else:
+            covar_list = [self.covar for m in range(n_pixels)]
+            covar = block_diag(covar_list, dtype=np.float32)
+            return x0, covar
+
+
+def propagate_LAI_narrowbandSAIL(x_analysis, P_analysis,
+                                     P_analysis_inverse,
+                                     M_matrix, Q_matrix,
+                                     prior=None, state_propagator=None, date=None):
+    ''' Propagate a single parameter and
+     set the rest of the parameter propagations to zero
+     This should be used with a prior for the remaining parameters'''
+    nparameters = 10
+    lai_position = 6
+
+    x_prior, c_prior, c_inv_prior = sail_prior_values()
+    return propagate_single_parameter(x_analysis, P_analysis,
+                                      P_analysis_inverse,
+                                      M_matrix, Q_matrix,
+                                      nparameters, lai_position,
+                                      x_prior, c_inv_prior)