implemented new CTDAS structure for the CTE-CO2 setup, including XCO2...

implemented new CTDAS structure for the CTE-CO2 setup, including XCO2 routines, plus some remaining changes between py2 -> py3

da/analysis/expand_fluxes.py

@@ -113,6 +113,10 @@ def save_weekly_avg_1x1_data(dacycle, statevector):
         fire = np.array(file.get_variable(dacycle.dasystem['background.co2.fires.flux']))
         fossil = np.array(file.get_variable(dacycle.dasystem['background.co2.fossil.flux']))
         #mapped_parameters   = np.array(file.get_variable(dacycle.dasystem['final.param.mean.1x1']))
+        if dacycle.dasystem['background.co2.bio.lt.flux'] in list(file.variables.keys()):
+            lt_flux = True
+            bio_lt  = np.array(file.get_variable(dacycle.dasystem['background.co2.bio.lt.flux']))
+        else: lt_flux = False
         if dacycle.dasystem['background.co2.biosam.flux'] in list(file.variables.keys()):
             sam = True
             biosam = np.array(file.get_variable(dacycle.dasystem['background.co2.biosam.flux']))
@@ -162,10 +166,13 @@ def save_weekly_avg_1x1_data(dacycle, statevector):
 #
 # if prior, do not multiply fluxes with parameters, otherwise do
 #
-            print(gridensemble.shape, bio.shape, gridmean.shape)
-            biomapped = bio * gridmean 
-            oceanmapped = ocean * gridmean 
-            biovarmapped = bio * gridensemble
+            if lt_flux:
+                biomapped    = bio_lt + bio * gridmean
+                biovarmapped = bio_lt + bio * gridensemble
+            else:
+                biomapped    = bio * gridmean 
+                biovarmapped = bio * gridensemble
+            oceanmapped    = ocean * gridmean 
             oceanvarmapped = ocean * gridensemble
 
 #

da/analysis/expand_molefractions.py

@@ -19,7 +19,6 @@ import shutil
 import logging
 import netCDF4
 import numpy as np
-from string import join
 from datetime import datetime, timedelta
 sys.path.append('../../')
 from da.tools.general import date2num, num2date
@@ -71,6 +70,10 @@ def write_mole_fractions(dacycle):
 
     ncf_in = io.ct_read(infile, 'read')
 
+    if len(ncf_in.dimensions['obs']) == 0:
+        ncf_in.close()
+        return None
+
     obs_num = ncf_in.get_variable('obs_num')
     obs_val = ncf_in.get_variable('observed')
     simulated = ncf_in.get_variable('modelsamples')
@@ -148,11 +151,10 @@ def write_mole_fractions(dacycle):
             ncf_out = io.CT_CDF(copy_file, 'write')
 
             # Modify the attributes of the file to reflect added data from CTDAS properly
-
-	    try:
-	       host=os.environ['HOSTNAME']
-	    except:
-	       host='unknown'
+            try:
+                host=os.environ['HOSTNAME']
+            except:
+                host='unknown'
      
 
             ncf_out.Caution = '==================================================================================='

da/analysis/summarize_obs.py

@@ -15,7 +15,7 @@ import sys
 sys.path.append('../../')
 import os
 import numpy as np
-import string
+#import string
 import datetime as dt
 import logging
 import re
@@ -36,9 +36,9 @@ def nice_lat(cls,format='html'):
 
     #return string.strip('%2d %2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h))
     if format == 'python':
-        return string.strip('%3d$^\circ$%2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h))
+        return ('%3d$^\circ$%2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h)).strip()
     if format == 'html':
-        return string.strip('%3d&deg%2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h))
+        return ('%3d&deg%2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h)).strip()
 
 def nice_lon(cls,format='html'):
     #
@@ -53,16 +53,16 @@ def nice_lon(cls,format='html'):
 
     #return string.strip('%3d %2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h))
     if format == 'python':
-        return string.strip('%3d$^\circ$%2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h))
+        return ('%3d$^\circ$%2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h)).strip()
     if format == 'html':
-        return string.strip('%3d&deg%2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h))
+        return ('%3d&deg%2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h)).strip()
 
 def nice_alt(cls):
     #
     # Reformat elevation or altitude
     #
     #return string.strip('%10.1f masl' % round(cls, -1))
-    return string.strip('%i masl' %cls)
+    return ('%i masl' %cls).strip()
 
 
 def summarize_obs(analysisdir, printfmt='html'):
@@ -282,7 +282,6 @@ def make_map(analysisdir): #makes a map of amount of assimilated observations pe
     import netCDF4 as cdf
     import matplotlib.pyplot as plt 
     import matplotlib
-    from maptools import *
     from matplotlib.font_manager import FontProperties
 
     sumdir = os.path.join(analysisdir, 'summary')

da/analysis/time_avg_fluxes.py

@@ -29,7 +29,7 @@ from dateutil.relativedelta import relativedelta
 import datetime
 import subprocess
 
-def time_avg(dacycle,avg='transcom'):
+def time_avg(dacycle,avg='transcom',monthly=False):
     """ Function to create a set of averaged files in a folder, needed to make longer term means """
     
     if avg not in ['transcom','transcom_extended','olson','olson_extended','country','flux1x1']:
@@ -40,9 +40,9 @@ def time_avg(dacycle,avg='transcom'):
     if not os.path.exists(analysisdir):
         raise IOError('analysis dir requested (%s) does not exist, exiting...'%analysisdir)
 
-    daily_avg(dacycle,avg)
-
-    monthly_avg(dacycle,avg)
+    if not monthly:
+        daily_avg(dacycle,avg)
+        monthly_avg(dacycle,avg)
 
     yearly_avg(dacycle,avg)
 

da/analysis/tools_regions.py

@@ -66,12 +66,12 @@ def state_to_grid(values, regionmap, reverse=False, avg=False, mapname=None):
     if reverse:
         """ project 1x1 degree map onto ecoregions """
 
-        result = np.zeros(nregions, float)
+        result = np.zeros(int(nregions), float)
         for k, v in regionselect.items():
-            if avg: 
-                result[k - 1] = values.ravel().take(v).mean()
-            else : 
-                result[k - 1] = values.ravel().take(v).sum()
+            if avg:
+                result[int(k) - 1] = values.ravel().take(v).mean()
+            else :
+                result[int(k) - 1] = values.ravel().take(v).sum()
         return result
 
     else:

da/analysis/tools_transcom.py

@@ -18,7 +18,6 @@ sys.path.append('../../')
 rootdir = os.getcwd().split('da/')[0]
 analysisdir = os.path.join(rootdir, 'da/analysis')
 
-from string import join, split
 from numpy import array, identity, zeros, arange, dot
 import da.tools.io4 as io
 
@@ -49,7 +48,7 @@ temp = open(os.path.join(analysisdir, 't3_region_names'), 'r').readlines()
 for line in temp:
     items = line.split()
     if items:
-        num, abbr, name = (items[0], items[1], join(items[2:]),)
+        num, abbr, name = (items[0], items[1], "".join(items[2:]))
         transnams.append(name.strip('"'))
         transshort.append(abbr)
         if abbr.startswith('T'):
@@ -66,7 +65,7 @@ temp = open(os.path.join(analysisdir, 'olson19_region_names'), 'r').readlines()
 for line in temp:
     items = line.split()
     if items:
-        num, abbr, name = (items[0], items[1], join(items[2:]),)
+        num, abbr, name = (items[0], items[1], "".join(items[2:]))
         olsonnams.append(name.strip('"'))
         olsonshort.append(abbr)
 
@@ -192,7 +191,7 @@ translocs = [ (-177, 0), \
             (68, 2,)]
 
 #olsonshort=[str(name.split()[1:2]).join('  ') for name in olsonnams]
-old_olsonshort = [join(split(name, ' ')[1:2], ' ') for name in olsonnams]
+old_olsonshort = ["".join((name.split())[1:2]) for name in olsonnams]
 
 olsonlabs = ['Conifer Forest', 'Broadleaf Forest', 'Mixed Forest', 'Grass/Shrub', 'Tropical Forest', 'Scrub/Woods', 'Semitundra', 'Fields/Woods/\nSavanna', \
      'Northern Taiga', 'Forest/Field', 'Wetland', 'Deserts', 'Shrub/Tree/\nSuc ', 'Crops', 'Conifer\n Snowy/Coastal', \

da/carbondioxide/dasystem.py → da/dasystems/dasystem_cteco2.py

@@ -25,7 +25,7 @@ import logging
 
 ################### Begin Class CO2DaSystem ###################
 
-from da.baseclasses.dasystem import DaSystem
+from da.dasystems.dasystem_baseclass import DaSystem
 
 class CO2DaSystem(DaSystem):
     """ Information on the data assimilation system used. This is normally an rc-file with settings.

da/tm5/observationoperator.py → da/obsoperators/observationoperator_tm5_cteco2.py

@@ -38,14 +38,14 @@ import logging
 import shutil
 import datetime
 import subprocess
-from string import join
+#from string import join
 import glob
 sys.path.append(os.getcwd())
 sys.path.append("../../")
 
 import da.tools.rc as rc
 from da.tools.general import create_dirs, to_datetime
-from da.baseclasses.observationoperator import ObservationOperator
+from da.obsoperators.observationoperator_baseclass import ObservationOperator
 
 identifier = 'TM5'
 version = 'release 3.0'
@@ -192,7 +192,7 @@ class TM5ObservationOperator(ObservationOperator):
         else: cmd = ['python', 'setup_tm5', '--%s' % self.dacycle.daplatform.give_queue_type(), rcfilename]
 
         logging.info('Starting the external TM5 setup script')
-        logging.info('using command ... %s' % join(cmd))
+        logging.info('using command ... %s' % str.join(str(cmd),''))
 
         retcode = subprocess.call(cmd)
         os.chdir(self.dacycle['dir.da_submit'])
@@ -204,8 +204,8 @@ class TM5ObservationOperator(ObservationOperator):
         else:
             logging.info('Compilation successful, continuing')
 
-    def prepare_run(self):
-        """ 
+    def prepare_run(self, samples):
+        """
         Prepare a forward model TM5 run, this consists of:
 
           - reading the working copy TM5 rc-file, 
@@ -226,10 +226,14 @@ class TM5ObservationOperator(ObservationOperator):
                     'jobstep.length': 'inf',
                     'ct.params.input.dir': self.dacycle['dir.input'],
                     'ct.params.input.file': os.path.join(self.dacycle['dir.input'], 'parameters'),
-                    'output.flask.infile': self.dacycle['ObsOperator.inputfile'] ,
-                    'output.flask': 'True'
                     }
 
+        # add sample type specific settings to rc-file
+        for sample in samples:
+            new_items['output.'+sample.get_samples_type()] = 'True'
+            new_items['output.'+sample.get_samples_type()+'.infile'] = self.dacycle['ObsOperator.inputfile.'+sample.get_samples_type()]
+        del sample
+
         if self.dacycle['transition']:
             new_items[self.istartkey] = self.transitionvalue
             logging.debug('Resetting TM5 to perform transition of od meteo from 25 to 34 levels') 
@@ -247,15 +251,17 @@ class TM5ObservationOperator(ObservationOperator):
         if self.dacycle['time.sample.window'] != 0:  # If this is a restart from a previous time step within the filter lag, the TM5 model should do a restart
             new_items[self.istartkey] = self.restartvalue
             logging.debug('Resetting TM5 to perform restart')
-        
-        # If neither one is true, simply take the istart value from the tm5.rc file that was read
 
+        # If neither one is true, simply take the istart value from the tm5.rc file that was read
         self.modify_rc(new_items)
         self.write_rc(self.rc_filename)
 
-	# Define the name of the file that will contain the modeled output of each observation
+	# For each sample type, define the name of the file that will contain the modeled output of each observation
+        self.simulated_file = [None] * len(samples)
+        for i in range(len(samples)):
+            self.simulated_file[i] = os.path.join(self.outputdir, '%s_output.%s.nc' % (samples[i].get_samples_type(),self.dacycle['time.sample.stamp']))
+        del i
 
-    	self.simulated_file = os.path.join(self.outputdir, 'flask_output.%s.nc' % self.dacycle['time.sample.stamp'])
 
     def load_rc(self, name):
         """ 
@@ -346,7 +352,7 @@ class TM5ObservationOperator(ObservationOperator):
         """
     
         for k, v in list(newvalues.items()):
-            if key in self.tm_settings:
+            if k in self.tm_settings:
                 # keep previous value
                 v_orig = self.tm_settings[k]
                 #replace with new
@@ -374,7 +380,9 @@ class TM5ObservationOperator(ObservationOperator):
         rc.write(tm5rcfilename, self.tm_settings)
         logging.debug("Modified rc file for TM5 written (%s)" % tm5rcfilename)
 
-    def validate_input(self):
+
+
+    def validate_input(self, samples):
         """
         Make sure that parameter files are written to the TM5 inputdir, and that observation lists are present
         """
@@ -387,13 +395,15 @@ class TM5ObservationOperator(ObservationOperator):
 
         datafiles = os.listdir(datadir)
 
-        obsfile = self.dacycle['ObsOperator.inputfile']
-
-        if not os.path.exists(obsfile):
-            msg = "The specified obs input file for the TM5 model to read from does not exist (%s), exiting..." % obsfile  
-            logging.error(msg)
-            if 'forward.savestate.dir' not in self.dacycle:
-                raise IOError(msg)
+        # For every sample type, check if input file is present
+        for sample in samples:
+            obsfile = self.dacycle['ObsOperator.inputfile.'+sample.get_samples_type()]
+            if not os.path.exists(obsfile) and len(sample.datalist) > 0:
+                msg = "The specified obs input file for the TM5 model to read from does not exist (%s), exiting..." % obsfile
+                logging.error(msg)
+                if 'forward.savestate.dir' not in self.dacycle:
+                    raise IOError(msg)
+        del sample
 
         for n in range(int(self.dacycle['da.optimizer.nmembers'])):
             paramfile = 'parameters.%03d.nc' % n
@@ -457,11 +467,15 @@ class TM5ObservationOperator(ObservationOperator):
             shutil.copy(fpath, fpath.replace(sourcedir, targetdir))
         logging.debug("All restart data have been copied from the restart/current directory to the TM5 save dir")
 
-    def run_forecast_model(self):
-        self.prepare_run()
-        self.validate_input()
+
+
+    def run_forecast_model(self, samples):
+        self.prepare_run(samples)
+        self.validate_input(samples)
         self.run()
-        self.save_data()
+        self.save_data(samples)
+
+
 
     def run(self):
         """ 
@@ -545,7 +559,9 @@ class TM5ObservationOperator(ObservationOperator):
 
         return code
 
-    def save_data(self):
+
+
+    def save_data(self, samples):
         """ Copy the TM5 recovery data from the outputdir to the TM5 savedir, also add the restart files to a list of names
             that is used by the dacycle object to collect restart data for the filter.
 
@@ -589,7 +605,12 @@ class TM5ObservationOperator(ObservationOperator):
 
         sourcedir = os.path.join(self.tm_settings[self.outputdirkey])
         sd_ed = self.dacycle['time.sample.stamp']
-        filterlist = ['flask_output.%s' % sd_ed, 'flux1x1_%s' % sd_ed]
+
+        filterlist = []
+        for sample in samples:
+            filterlist.append('%s_output.%s' % (sample.get_samples_type(),sd_ed))
+        del sample
+        filterlist.append('flux1x1_%s' % sd_ed)
 
         logging.debug("Creating a new list of TM5 output data to collect")
         logging.debug("           from directory: %s " % sourcedir)

da/baseclasses/optimizer.py → da/optimizers/optimizer_baseclass.py

@@ -86,6 +86,8 @@ class Optimizer(object):
         self.species = np.zeros(self.nobs, str)
         # species type
         self.sitecode = np.zeros(self.nobs, str)
+        # rejection_threshold
+        self.rejection_threshold = np.zeros(self.nobs, float)
 
         # species mask
         self.speciesmask = {}
@@ -104,6 +106,7 @@ class Optimizer(object):
         allflags = []  # collect all model samples for n=1,..,nlag
         allspecies = []  # collect all model samples for n=1,..,nlag
         allsimulated = []  # collect all members model samples for n=1,..,nlag
+        allrej_thres = []  # collect all rejection_thresholds, will be the same for all samples of same source
 
         for n in range(self.nlag):
             samples = statevector.obs_to_assimilate[n]
@@ -111,22 +114,28 @@ class Optimizer(object):
             self.x[n * self.nparams:(n + 1) * self.nparams] = members[0].param_values
             self.X_prime[n * self.nparams:(n + 1) * self.nparams, :] = np.transpose(np.array([m.param_values for m in members]))
 
-            if samples != None:        
-                self.rejection_threshold = samples.rejection_threshold
-
-                allreject.extend(samples.getvalues('may_reject'))
-                alllocalize.extend(samples.getvalues('may_localize'))
-                allflags.extend(samples.getvalues('flag'))
-                allspecies.extend(samples.getvalues('species'))
-                allobs.extend(samples.getvalues('obs'))
-                allsites.extend(samples.getvalues('code'))
-                allmdm.extend(samples.getvalues('mdm'))
-                allids.extend(samples.getvalues('id'))
-
-                simulatedensemble = samples.getvalues('simulated')
-                for s in range(simulatedensemble.shape[0]):
+            # Add observation data for all sample objects
+            if samples != None:
+                if type(samples) != list: samples = [samples]
+                for m in range(len(samples)):
+                    sample = samples[m]
+                    logging.debug('Lag %i, sample %i: rejection_threshold = %i, nobs = %i' %(n, m, sample.rejection_threshold, sample.getlength()))
+
+                    allrej_thres.extend([sample.rejection_threshold] * sample.getlength())
+                    allreject.extend(sample.getvalues('may_reject'))
+                    alllocalize.extend(sample.getvalues('may_localize'))
+                    allflags.extend(sample.getvalues('flag'))
+                    allspecies.extend(sample.getvalues('species'))
+                    allobs.extend(sample.getvalues('obs'))
+                    allsites.extend(sample.getvalues('code'))
+                    allmdm.extend(sample.getvalues('mdm'))
+                    allids.extend(sample.getvalues('id'))
+
+                    simulatedensemble = sample.getvalues('simulated')
+                    for s in range(simulatedensemble.shape[0]):
                         allsimulated.append(simulatedensemble[s])
 
+        self.rejection_threshold[:] = np.array(allrej_thres)
         self.obs[:] = np.array(allobs)
         self.obs_ids[:] = np.array(allids)
         self.HX_prime[:, :] = np.array(allsimulated)
@@ -147,12 +156,12 @@ class Optimizer(object):
         else:
             for i, mdm in enumerate(allmdm):
                 self.R[i, i] = mdm ** 2
-                    
+
     def matrix_to_state(self, statevector):
         for n in range(self.nlag):
             members = statevector.ensemble_members[n]
             for m, mem in enumerate(members):
-                members[m].param_values[:] = self.X_prime[n * self.nparams:(n + 1) * self.nparams, m] + self.x[n * self.nparams:(n + 1) * self.nparams]     
+                members[m].param_values[:] = self.X_prime[n * self.nparams:(n + 1) * self.nparams, m] + self.x[n * self.nparams:(n + 1) * self.nparams]
 
         logging.debug('Returning optimized data to the StateVector, setting "StateVector.isOptimized = True" ')
 
@@ -309,6 +318,12 @@ class Optimizer(object):
 
     def serial_minimum_least_squares(self):
         """ Make minimum least squares solution by looping over obs"""
+
+        # Calculate prior value cost function (observation part)
+        res_prior = np.abs(self.obs-self.Hx)
+        select    = (res_prior < 1E15).nonzero()[0]
+        J_prior   = res_prior.take(select,axis=0)**2/self.R.take(select,axis=0)
+        res_prior = np.mean(res_prior)
         for n in range(self.nobs):
 
             # Screen for flagged observations (for instance site not found, or no sample written from model)
@@ -322,7 +337,7 @@ class Optimizer(object):
             res = self.obs[n] - self.Hx[n]
 
             if self.may_reject[n]:
-                threshold = self.rejection_threshold * np.sqrt(self.R[n])
+                threshold = self.rejection_threshold[n] * 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
@@ -330,10 +345,9 @@ class Optimizer(object):
 
             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, :])
+            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[:]   = PHt / self.HPHR[n]
 
             if self.may_localize[n]:
                 logging.debug('Trying to localize observation %s, %i' % (self.sitecode[n], self.obs_ids[n]))
@@ -347,21 +361,43 @@ class Optimizer(object):
 
             for r in range(self.nmembers):
                 self.X_prime[:, r] = self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r])
+            del r
+
+            # update samples to account for update of statevector based on observation n
+            HXprime_n     = self.HX_prime[n,:].copy()
+            res           = self.obs[n] - self.Hx[n]
+            fac           = 1.0 / (self.nmembers - 1) * np.sum(HXprime_n[np.newaxis,:] * self.HX_prime, axis=1) / self.HPHR[n]
+            self.Hx       = self.Hx + fac*res
+            self.HX_prime = self.HX_prime - alpha* fac[:,np.newaxis]*HXprime_n
+
+
+        del n
+        if 'HXprime_n' in globals(): del HXprime_n
+
+        # calculate posterior value cost function
+        res_post = np.abs(self.obs-self.Hx)
+        select   = (res_post < 1E15).nonzero()[0]
+        J_post   = res_post.take(select,axis=0)**2/self.R.take(select,axis=0)
+        res_post = np.mean(res_post)
+
+        logging.info('Observation part cost function: prior = %s, posterior = %s' % (np.mean(J_prior), np.mean(J_post)))
+        logging.info('Mean residual: prior = %s, posterior = %s' % (res_prior, res_post))
 
 #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, :]
 
-            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, :]
 
             
     def bulk_minimum_least_squares(self):

da/carbondioxide/optimizer.py → da/optimizers/optimizer_cteco2.py

@@ -24,8 +24,10 @@ File created on 28 Jul 2010.
 import os
 import sys
 import logging
+import numpy as np
+import da.tools.io4 as io
 sys.path.append(os.getcwd())
-from da.baseclasses.optimizer import Optimizer
+from da.optimizers.optimizer_baseclass import Optimizer
 
 
 identifier = 'Ensemble Square Root Filter'
@@ -41,6 +43,132 @@ class CO2Optimizer(Optimizer):
         All other methods are inherited from the base class Optimizer.
     """
 
+    def create_matrices(self):
+        """ Create Matrix space needed in optimization routine """
+
+        # mean state  [X]
+        self.x = np.zeros((self.nlag * self.nparams,), float)
+        # deviations from mean state  [X']
+        self.X_prime = np.zeros((self.nlag * self.nparams, self.nmembers,), float)
+        # mean state, transported to observation space [ H(X) ]
+        self.Hx = np.zeros((self.nobs,), float)
+        # deviations from mean state, transported to observation space [ H(X') ]
+        self.HX_prime = np.zeros((self.nobs, self.nmembers), float)
+        # observations
+        self.obs = np.zeros((self.nobs,), float)
+        # observation ids
+        self.obs_ids = np.zeros((self.nobs,), float)
+        # covariance of observations
+        # Total covariance of fluxes and obs in units of obs [H P H^t + R]
+        if self.algorithm == 'Serial':
+            self.R = np.zeros((self.nobs,), float)
+            self.HPHR = np.zeros((self.nobs,), float)
+        else:
+            self.R = np.zeros((self.nobs, self.nobs,), float)
+            self.HPHR = np.zeros((self.nobs, self.nobs,), float)
+        # localization of obs
+        self.may_localize = np.zeros(self.nobs, bool)
+        self.loc_L = np.zeros(self.nobs, int)
+        # rejection of obs
+        self.may_reject = np.zeros(self.nobs, bool)
+        # flags of obs
+        self.flags = np.zeros(self.nobs, int)
+        # species type
+        self.species = np.zeros(self.nobs, str)
+        # species type
+        self.sitecode = np.zeros(self.nobs, str)
+        # rejection_threshold
+        self.rejection_threshold = np.zeros(self.nobs, float)
+        # lat and lon
+        self.latitude = np.zeros(self.nobs, float)
+        self.longitude = np.zeros(self.nobs, float)
+        self.date = np.zeros((self.nobs, 6), float)
+
+        # species mask
+        self.speciesmask = {}
+