diff --git a/kafka/kf_tools.py b/kafka/kf_tools.py index 4247678..71ddf91 100644 --- a/kafka/kf_tools.py +++ b/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) diff --git a/kafka/linear_kf.py b/kafka/linear_kf.py index 4213803..ed03907 100644 --- a/kafka/linear_kf.py +++ b/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... diff --git a/kafka/observations.py b/kafka/observations.py index 3ffd205..688aded 100644 --- a/kafka/observations.py +++ b/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 diff --git a/kafka/solvers.py b/kafka/solvers.py index ad8c2f1..04a95d2 100644 --- a/kafka/solvers.py +++ b/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 diff --git a/tests/test_kf.py b/tests/test_kf.py index c867dd8..d03d8a6 100644 --- a/tests/test_kf.py +++ b/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(),