From 2a11c77b93a6a39519351c8e74721714522eb67b Mon Sep 17 00:00:00 2001 From: NPounder Date: Wed, 13 Dec 2017 16:01:40 +0000 Subject: [PATCH 1/4] remove imports of unused code nolonger in repository --- kafka/linear_kf.py | 5 +++-- kafka/observations.py | 6 +++--- kafka/solvers.py | 2 +- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/kafka/linear_kf.py b/kafka/linear_kf.py index 4213803..1e04c2b 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 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 @@ -269,6 +269,7 @@ def assimilate(self, locate_times, x_forecast, P_forecast, x_forecast = x_analysis*1 P_forecast = P_analysis + P_forecast_inverse = P_analysis_inverse if iter_obs_op: 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 From 6bd18247daa623f4ea15b245342d1fded70f121f Mon Sep 17 00:00:00 2001 From: NPounder Date: Thu, 21 Dec 2017 12:25:27 +0000 Subject: [PATCH 2/4] correction to Hessian calculation --- kafka/kf_tools.py | 37 +++++++++++++++++++++++++++++++++++-- kafka/linear_kf.py | 12 +++++++++--- 2 files changed, 44 insertions(+), 5 deletions(-) diff --git a/kafka/kf_tools.py b/kafka/kf_tools.py index 4247678..bf277c6 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,P_inv, innovation, mask, state_mask, band, nparams): + C_obs_inv = P_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 1e04c2b..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_SLOW +from kf_tools import hessian_correction, propagate_information_filter_SLOW # Set up logging @@ -267,10 +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... From d8a12c59987e3f0d2aacfdc1677107f9cac9772e Mon Sep 17 00:00:00 2001 From: NPounder Date: Tue, 2 Jan 2018 16:50:57 +0000 Subject: [PATCH 3/4] update test with correct method name (propagate_information_filter_SLOW) --- tests/test_kf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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(), From 540c2bf2336c76470c3a65ea9cefe5c1c9f15ebd Mon Sep 17 00:00:00 2001 From: Nicola Pounder Date: Tue, 9 Jan 2018 15:58:31 +0000 Subject: [PATCH 4/4] Change variable name to reduce confusion. --- kafka/kf_tools.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/kafka/kf_tools.py b/kafka/kf_tools.py index bf277c6..71ddf91 100644 --- a/kafka/kf_tools.py +++ b/kafka/kf_tools.py @@ -29,8 +29,8 @@ def hessian_correction_pixel(gp, x0, C_obs_inv, innovation, band, nparams): return big_Hessian_corr -def hessian_correction(gp, x0,P_inv, innovation, mask, state_mask, band, nparams): - C_obs_inv = P_inv.diagonal()[state_mask.flatten()] +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)):