Skip to content
Snippets Groups Projects
Commit 0e7df5f6 authored by weihe's avatar weihe
Browse files

revised additive scheme

parent 27811f35
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 59 additions and 43 deletions
da/stilt/optimizer.py
+ 59
− 43
View file @ 0e7df5f6
"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
updates of the code. See also: http://www.carbontracker.eu.
This program is free software: you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software Foundation,
version 3. This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this
program. If not, see <http://www.gnu.org/licenses/>."""
#!/usr/bin/env python
# optimizer.py
"""
Author : peters
Author : peters
Revision History:
File created on 28 Jul 2010.
......@@ -55,12 +43,12 @@ class CO2Optimizer(Optimizer):
elif self.nmembers == 150:
self.tvalue = 1.97591
elif self.nmembers == 200:
self.tvalue = 1.9719
else: self.tvalue = 0
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:
......@@ -72,9 +60,9 @@ class CO2Optimizer(Optimizer):
""" localize the Kalman Gain matrix """
import numpy as np
if not self.localization:
if not self.localization:
logging.debug('Not localized observation %i' % self.obs_ids[n])
return
return
if self.localizetype == 'CT2007':
count_localized = 0
for r in range(self.nlag * self.nparams):
......@@ -92,7 +80,7 @@ class CO2Optimizer(Optimizer):
self.algorithm = 'Serial'
else:
self.algorithm = 'Bulk'
logging.info("Current minimum least squares algorithm is set to %s" % self.algorithm)
......@@ -101,17 +89,20 @@ class CO2Optimizer(Optimizer):
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):
#for n in range(self.nobs):
res = self.obs[n] - self.Hx[n]
# 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
# 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
#test=np.zeros((self.Hx.shape[0]),)
#test[:]=self.Hx[:]
#logging.info("Total HX array: %s"%test)
for n in range(self.nobs):
......@@ -121,24 +112,34 @@ class CO2Optimizer(Optimizer):
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]
self.KG[:,n] = PHt / self.HPHR[n]
#if 'surface' in self.sitecode[n]:
# self.KG[-4:,n]=0.
# self.KG[3078:3082,n]=0.
# logging.debug('BC KG value set to zero for %s' %(self.sitecode[n]))
#if 'aircraft' in self.sitecode[n]:
# self.KG[0:3078,n]=0.
# self.KG[3082:-4,n]=0.
# logging.debug('Flux KG values set to zero for %s' %(self.sitecode[n]))
if 'surface' in self.sitecode[n]:
self.KG[-1]=0.
self.KG[3078]=0
logging.debug('BC KG value set to zero for %s' %(self.sitecode[n]))
if 'aircraft' in self.sitecode[n]:
self.KG[0:3078]=0.
self.KG[3079:-1]=0.
logging.debug('Flux KG values set to zero for %s' %(self.sitecode[n]))
if self.may_localize[n]:
logging.debug('Trying to localize observation %s, %i' % (self.sitecode[n], self.obs_ids[n]))
......@@ -149,15 +150,16 @@ class CO2Optimizer(Optimizer):
alpha = np.double(1.0) / (np.double(1.0) + np.sqrt((self.R[n]) / self.HPHR[n]))
self.x[:] = self.x + self.KG[:] * res
logging.debug('Residual %s'%res)
logging.debug('New self.KG BC1 %s' %self.KG[3078])
logging.debug('New self.KG BC2 %s' %self.KG[-1])
logging.debug('New self.x BC1 %s' %self.x[3078])
logging.debug('New self.x BC2 %s' %self.x[-1])
self.x[:] = self.x + self.KG[:,n] * res
#logging.debug('Residual %s'%res)
#logging.debug('obs = %s, Hx = %s,Hx(CAR) = %s, Hxp(CAR) = %s, Hx_prime_std(CAR) = %s'%(self.obs[n],self.Hx[n],self.Hx[19],test[19],self.HX_prime[19,:].std()))
#logging.debug('New self.KG BC1 %s' %self.KG[3078,n])
#logging.debug('New self.KG BC2 %s' %self.KG[-1,n])
#logging.debug('New self.x BC1 %s' %self.x[3078])
#logging.debug('New self.x BC2 %s' %self.x[-1])
for r in range(self.nmembers):
self.X_prime[:, r] = self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r])
self.X_prime[:, r] = self.X_prime[:, r] - alpha * self.KG[:,n] * (self.HX_prime[n, r])
......@@ -165,17 +167,31 @@ class CO2Optimizer(Optimizer):
#WP should always be updated last because it features in the loop of the adjustments !!!!
for m in range(n + 1, self.nobs):
#if 'aircraft' in self.sitecode[m] and not 'aircraft' in self.sitecode[n]:
# continue
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
#if n==0 and m==19:
# logging.debug('self.HX_prime[n, :]= %s'%self.HX_prime[n, :])
# logging.debug('self.HX_prime[m, :]= %s'%self.HX_prime[m, :])
self.HX_prime[m, :] = self.HX_prime[m, :] - alpha * fac * self.HX_prime[n, :]
#logging.debug('m = %s, Corrcoef = %s, fac = %s'%(m, np.corrcoef(self.HX_prime[n, :],self.HX_prime[m, :])[0,1],fac))
for m in range(1, n + 1):
#if 'aircraft' in self.sitecode[m] and not 'aircraft' in self.sitecode[n]:
# continue
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, :]
# logging.debug('m = %s, Corrcoef = %s, fac = %s'%(m, np.corrcoef(self.HX_prime[n, :],self.HX_prime[m, :])[0,1],fac))
################### End Class CO2Optimizer ###################
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment