sf6 code to be updated

da/sf6/statevector.py

+#!/usr/bin/env python
+# ct_statevector_tools.py
+
+"""
+Author : peters 
+
+Revision History:
+File created on 28 Jul 2010.
+
+"""
+
+import os
+import sys
+sys.path.append(os.getcwd())
+
+import logging
+from da.baseclasses.statevector import StateVector
+import numpy as np
+import da.tools.io4 as io
+
+identifier = 'SF6 Statevector '
+version = '0.0'
+
+################### Begin Class SF6StateVector ###################
+
+class SF6StateVector(StateVector):
+    """ This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """
+
+    def get_covariance(self, date):
+        """ Make a new ensemble from specified matrices, the attribute lag refers to the position in the state vector. 
+            Note that lag=1 means an index of 0 in python, hence the notation lag-1 in the indexing below.
+            The argument is thus referring to the lagged state vector as [1,2,3,4,5,..., nlag]
+        """    
+
+        import da.tools.io4 as io
+        try:
+            import matplotlib.pyplot as plt
+        except:
+            pass
+
+        # Get the needed matrices from the specified covariance files
+
+        fullcov = np.zeros((self.nparams, self.nparams), float)
+
+        fullcov[:,:] = 0.2**2
+
+        return fullcov
+
+    def Initialize(self):
+        """
+        Initialize the object by specifying the dimensions. 
+        There are two major requirements for each statvector that you want to build:
+        
+            (1) is that the statevector can map itself onto a regular grid
+            (2) is that the statevector can map itself (mean+covariance) onto TransCom regions
+
+        An example is given below.
+        """
+
+        self.nlag = int(self.DaCycle['time.nlag'])
+        self.nmembers = int(self.DaCycle['da.optimizer.nmembers'])
+        self.nparams = int(self.DaCycle.DaSystem['nparameters'])
+        self.nobs = 0
+        
+        self.isOptimized = False
+        self.ObsToAssimmilate = ()  # empty containter to hold observations to assimilate later on
+
+        # These list objects hold the data for each time step of lag in the system. Note that the ensembles for each time step consist 
+        # of lists of EnsembleMember objects, we define member 0 as the mean of the distribution and n=1,...,nmembers as the spread.
+
+        self.EnsembleMembers = range(self.nlag)
+
+        for n in range(self.nlag):
+            self.EnsembleMembers[n] = []
+
+
+        # This specifies the file to read with the gridded mask at 1x1 degrees. Each gridbox holds a number that specifies the parametermember
+        #  that maps onto it. From this map, a dictionary is created that allows a reverse look-up so that we can map parameters to a grid.
+
+        mapfile = os.path.join(self.DaCycle.DaSystem['regionsfile'])
+        ncf = io.CT_Read(mapfile, 'read')
+        self.tcmap = ncf.GetVariable('transcom_regions')
+        ncf.close()
+
+        self.gridmap = np.ones((180,360),'float')
+
+        logging.debug("A TransCom  map on 1x1 degree was read from file %s" % self.DaCycle.DaSystem['regionsfile'])
+        logging.debug("A parameter map on 1x1 degree was created")
+
+        # Create a dictionary for state <-> gridded map conversions
+
+        nparams = self.gridmap.max()
+        self.griddict = {}
+        for r in range(1, nparams + 1):
+            sel = (self.gridmap.flat == r).nonzero()
+            if len(sel[0]) > 0: 
+                self.griddict[r] = sel
+
+        logging.debug("A dictionary to map grids to states and vice versa was created")
+
+        # Create a matrix for state <-> TransCom conversions
+
+        self.tcmatrix = np.zeros((self.nparams, 23), 'float') 
+
+        logging.debug("A matrix to map states to TransCom regions and vice versa was created")
+
+        # Create a mask for species/unknowns
+
+        self.make_species_mask()
+
+    def make_species_mask(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 = {'sf6': np.ones(self.nparams)}
+        logging.debug("A species mask was created, only the following species are recognized in this system:")
+        for k in self.speciesdict.keys(): 
+            logging.debug("   ->    %s" % k)
+
+    def propagate(self):
+        """ 
+        :rtype: None
+
+        Propagate the parameter values in the StateVector to the next cycle. This means a shift by one cycle 
+        step for all states that will
+        be optimized once more, and the creation of a new ensemble for the time step that just 
+        comes in for the first time (step=nlag). 
+        In the future, this routine can incorporate a formal propagation of the statevector.
+
+        """
+        from datetime import timedelta
+
+        ensemble_deviations = np.array([p.ParameterValues for p in self.EnsembleMembers[0] ])
+        if ensemble_deviations.std(axis=0) < 0.05 :
+            for m in self.EnsembleMembers[0]:
+                m  = 3.0 * m  # inflate deviations by a factor of 3.0
+
+            logging.info('The state vector covariance was inflated to 1-sigma of %5.2f ' % (ensemble_deviations.std(axis=0)*3.0))
+
+        logging.info('The state vector remains the same in the SF6 run')
+        logging.info('The state vector has been propagated by one cycle')
+
+
+################### End Class SF6StateVector ###################
+
+
+if __name__ == "__main__":
+    from da.tools.initexit import StartLogger 
+    from da.tools.pipeline import start_job
+
+    sys.path.append(os.getcwd())
+
+    opts = ['-v']
+    args = {'rc':'da.rc', 'logfile':'da_initexit.log', 'jobrcfilename':'test.rc'}
+
+    StartLogger()
+
+    DaCycle = start_job(opts, args)
+
+    DaCycle.Initialize()
+
+    StateVector = CtStateVector()
+
+    StateVector.Initialize()
+
+    for n in range(dims[0]):
+        cov = StateVector.get_covariance()
+        dummy = StateVector.make_new_ensemble(n + 1, cov)
+
+    StateVector.propagate()
+
+    savedir = DaCycle['dir.output']
+    filename = os.path.join(savedir, 'savestate.nc')
+
+    dummy = StateVector.WriteToFile()
+
+    StateVector.ReadFromFile(filename)
+