From 2afbd7622ab4aa9c3ae0c4c0949c873d5d99e087 Mon Sep 17 00:00:00 2001 From: ivar <amvdw95@gmail.com> Date: Thu, 13 Aug 2015 14:24:22 +0000 Subject: [PATCH] added new optimizer with Rejection flagging bug fixed --- da/stilt/optimizer.py | 167 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100755 da/stilt/optimizer.py diff --git a/da/stilt/optimizer.py b/da/stilt/optimizer.py new file mode 100755 index 0000000..1fe8a30 --- /dev/null +++ b/da/stilt/optimizer.py @@ -0,0 +1,167 @@ +#!/usr/bin/env python +# optimizer.py + +""" +Author : peters + +Revision History: +File created on 28 Jul 2010. + +""" + +import os +import sys +import logging +sys.path.append(os.getcwd()) +from da.baseclasses.optimizer import Optimizer + + +identifier = 'Ensemble Square Root Filter' +version = '0.0' + +################### Begin Class CO2Optimizer ################### + +class CO2Optimizer(Optimizer): + """ + This creates an instance of a CarbonTracker optimization object. The base class it derives from is the optimizer object. + Additionally, this CO2Optimizer implements a special localization option following the CT2007 method. + + All other methods are inherited from the base class Optimizer. + """ + + def set_localization(self, loctype='None'): + """ determine which localization to use """ + + if loctype == 'CT2007': + self.localization = True + self.localizetype = 'CT2007' + #T-test values for two-tailed student's T-test using 95% confidence interval for some options of nmembers + if self.nmembers == 50: + self.tvalue = 2.0086 + elif self.nmembers == 100: + self.tvalue = 1.9840 + elif self.nmembers == 150: + self.tvalue = 1.97591 + elif self.nmembers == 200: + self.tvalue = 1.9719 + else: self.tvalue = 0 + else: + self.localization = False + self.localizetype = 'None' + + logging.info("Current localization option is set to %s" % self.localizetype) + if self.localization == True: + if self.tvalue == 0: + logging.error("Critical tvalue for localization not set for %i ensemble members"%(self.nmembers)) + sys.exit(2) + else: logging.info("Used critical tvalue %0.05f is based on 95%% probability and %i ensemble members in a two-tailed student's T-test"%(self.tvalue,self.nmembers)) + + def localize(self, n): + """ localize the Kalman Gain matrix """ + import numpy as np + + if not self.localization: + logging.debug('Not localized observation %i' % self.obs_ids[n]) + return + if self.localizetype == 'CT2007': + count_localized = 0 + for r in range(self.nlag * self.nparams): + corr = np.corrcoef(self.HX_prime[n, :], self.X_prime[r, :].squeeze())[0, 1] + prob = corr / np.sqrt((1.000000001 - corr ** 2) / (self.nmembers - 2)) + if abs(prob) < self.tvalue: + self.KG[r] = 0.0 + count_localized = count_localized + 1 + logging.debug('Localized observation %i, %i%% of values set to 0' % (self.obs_ids[n],count_localized*100/(self.nlag * self.nparams))) + + def set_algorithm(self, algorithm='Serial'): + """ determine which minimum least squares algorithm to use """ + + if algorithm == 'Serial': + self.algorithm = 'Serial' + else: + self.algorithm = 'Bulk' + + logging.info("Current minimum least squares algorithm is set to %s" % self.algorithm) + + + + + def serial_minimum_least_squares(self): + """ Make minimum least squares solution by looping over obs""" + import numpy as np + for n in range(self.nobs): + + res = self.obs[n] - self.Hx[n] + + if self.may_reject[n]: + threshold = self.rejection_threshold * np.sqrt(self.R[n]) + if np.abs(res) > threshold: + logging.debug('Rejecting observation (%s,%i) because residual (%f) exceeds threshold (%f)' % (self.sitecode[n], self.obs_ids[n], res, threshold)) + self.flags[n] = 2 + continue + + + for n in range(self.nobs): + + # Screen for flagged observations (for instance site not found, or no sample written from model) + + if self.flags[n] != 0: + logging.debug('Skipping observation (%s,%i) because of flag value %d' % (self.sitecode[n], self.obs_ids[n], self.flags[n])) + continue + + + # Screen for outliers greather than 3x model-data mismatch, only apply if obs may be rejected + + #res = self.obs[n] - self.Hx[n] + + #if self.may_reject[n]: + # threshold = self.rejection_threshold * np.sqrt(self.R[n]) + # if np.abs(res) > threshold: + # logging.debug('Rejecting observation (%s,%i) because residual (%f) exceeds threshold (%f)' % (self.sitecode[n], self.obs_ids[n], res, threshold)) + # self.flags[n] = 2 + # continue + + + logging.debug('Proceeding to assimilate observation %s, %i' % (self.sitecode[n], self.obs_ids[n])) + + PHt = 1. / (self.nmembers - 1) * np.dot(self.X_prime, self.HX_prime[n, :]) + self.HPHR[n] = 1. / (self.nmembers - 1) * (self.HX_prime[n, :] * self.HX_prime[n, :]).sum() + self.R[n] + + self.KG[:] = PHt / self.HPHR[n] + + if self.may_localize[n]: + logging.debug('Trying to localize observation %s, %i' % (self.sitecode[n], self.obs_ids[n])) + self.localize(n) + else: + logging.debug('Not allowed to localize observation %s, %i' % (self.sitecode[n], self.obs_ids[n])) + + alpha = np.double(1.0) / (np.double(1.0) + np.sqrt((self.R[n]) / self.HPHR[n])) + + self.x[:] = self.x + self.KG[:] * res + + for r in range(self.nmembers): + self.X_prime[:, r] = self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r]) + + + +#WP !!!! Very important to first do all obervations from n=1 through the end, and only then update 1,...,n. The current observation +#WP should always be updated last because it features in the loop of the adjustments !!!! + + for m in range(n + 1, self.nobs): + res = self.obs[n] - self.Hx[n] + fac = 1.0 / (self.nmembers - 1) * (self.HX_prime[n, :] * self.HX_prime[m, :]).sum() / self.HPHR[n] + self.Hx[m] = self.Hx[m] + fac * res + self.HX_prime[m, :] = self.HX_prime[m, :] - alpha * fac * self.HX_prime[n, :] + + for m in range(1, n + 1): + res = self.obs[n] - self.Hx[n] + fac = 1.0 / (self.nmembers - 1) * (self.HX_prime[n, :] * self.HX_prime[m, :]).sum() / self.HPHR[n] + self.Hx[m] = self.Hx[m] + fac * res + self.HX_prime[m, :] = self.HX_prime[m, :] - alpha * fac * self.HX_prime[n, :] + + + +################### End Class CO2Optimizer ################### + +if __name__ == "__main__": + pass -- GitLab