From e5b81d1484aa68c224157d6b45bdd614df4d8d83 Mon Sep 17 00:00:00 2001
From: Wouter Peters <wouter.peters@wur.nl>
Date: Mon, 23 Apr 2012 15:35:45 +0000
Subject: [PATCH] species mask is introduced in code

---
 da/baseclasses/optimizer.py   |  9 ++++++++-
 da/baseclasses/statevector.py | 33 ++++++++++++++++++++++++++++++++-
 2 files changed, 40 insertions(+), 2 deletions(-)

diff --git a/da/baseclasses/optimizer.py b/da/baseclasses/optimizer.py
index 57cc7297..8f84fc62 100755
--- a/da/baseclasses/optimizer.py
+++ b/da/baseclasses/optimizer.py
@@ -72,6 +72,11 @@ class Optimizer(object):
         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 mask
+        self.speciesmask    = {}
 
         # Total covariance of fluxes and obs in units of obs [H P H^t + R]
         self.HPHR            =  np.zeros( (self.nobs,self.nobs,), float)
@@ -91,6 +96,7 @@ class Optimizer(object):
         allreject=[]  # collect all model samples for n=1,..,nlag
         alllocalize=[]  # collect all model samples for n=1,..,nlag
         allflags=[]  # collect all model samples for n=1,..,nlag
+        allspecies=[]  # collect all model samples for n=1,..,nlag
         allmemsamples=None  # collect all members model samples for n=1,..,nlag
 
         for n in range(self.nlag):
@@ -106,7 +112,7 @@ class Optimizer(object):
                 allreject.extend(                                 members[0].ModelSample.Data.getvalues('may_reject') )
                 alllocalize.extend(                               members[0].ModelSample.Data.getvalues('may_localize') )
                 allflags.extend(                                  members[0].ModelSample.Data.getvalues('flag') )
-
+                allspecies.extend(                                members[0].ModelSample.Data.getvalues('species') )
                 allobs.extend(                                    members[0].ModelSample.Data.getvalues('obs') )
                 allsamples.extend(                                members[0].ModelSample.Data.getvalues('simulated') )
                 allmdm.extend(                                    members[0].ModelSample.Data.getvalues('mdm')  )
@@ -129,6 +135,7 @@ class Optimizer(object):
         self.may_reject[:]                                      = np.array(allreject)
         self.may_localize[:]                                    = np.array(alllocalize)
         self.flags[:]                                           = np.array(allflags)
+        self.species[:]                                         = np.array(allspecies)
 
 
         self.X_prime                                            = self.X_prime - self.x[:,np.newaxis] # make into a deviation matrix
diff --git a/da/baseclasses/statevector.py b/da/baseclasses/statevector.py
index a8ac4d0e..2872d883 100755
--- a/da/baseclasses/statevector.py
+++ b/da/baseclasses/statevector.py
@@ -220,9 +220,41 @@ class StateVector(object):
 
         msg = "A matrix to map states to TransCom regions and vice versa was created" ; logging.debug(msg)
 
+        # Create a mask for species/unknowns
+
+        dummy = self.MakeSpeciesMask()
+
     def __str__(self):
         return "This is a base class derived state vector object"
 
+    def MakeSpeciesMask(self):
+
+        """
+
+        This method creates a dictionary with as key the name of a tracer, and as values an array of 0.0/1.0 values 
+        specifying which StateVector elements are constrained by this tracer. This mask can be used in 
+        the optimization to ensure that certain types of osbervations only update certain unknowns.
+
+        An example would be that the tracer '14CO2' can be allowed to only map onto fossil fuel emissions in the state
+
+        The form of the mask is:
+
+        {'co2': np.ones(self.nparams), 'co2c14', np.zeros(self.nparams)  }
+
+        so that 'co2' maps onto all parameters, and 'co2c14' on none at all. These arrays are used in the Class 
+        optimizer when state updates are actually performed
+
+        """
+
+        self.speciesdict = {'co2': np.ones(self.nparams)}
+
+        msg = "A species mask was created, only the following species are recognized in this system:" ; logging.debug(msg)
+        for k in self.speciesdict.keys(): 
+            msg = "   ->    %s"%k ; logging.debug(msg)
+
+        return None
+
+
     def MakeNewEnsemble(self,lag, covariancematrix = None):
         """ 
         :param lag: an integer indicating the time step in the lag order
@@ -409,7 +441,6 @@ class StateVector(object):
         EnsembleMembers = f.GetVariable('statevectorensemble_'+qual)
         dummy           = f.close()
 
-    
         for n in range(self.nlag):
             if not self.EnsembleMembers[n] == []:
                 self.EnsembleMembers[n] = []
-- 
GitLab