needs dasystem, first test of statevector

da/ctgridded/dasystem.py

+#!/usr/bin/env python
+# control.py
+
+"""
+Author : peters 
+
+Revision History:
+File created on 26 Aug 2010.
+
+"""
+
+import os
+import sys
+import logging
+import datetime
+
+################### Begin Class CtDaSystem ###################
+
+from da.baseclasses.dasystem import DaSystem
+
+class CtDaSystem(DaSystem):
+    """ Information on the data assimilation system used. This is normally an rc-file with settings.
+    """
+
+    def Initialize(self):
+        """ 
+        Initialize the object
+        """
+        import da.tools.io4 as io
+
+        mapfile   = os.path.join(self['datadir'],self['regionsfile'])
+        ncf       = io.CT_Read(mapfile,'read')
+        self.regionmap = ncf.GetVariable('budget_region')
+        dummy     = ncf.close()
+
+
+    def Validate(self):
+        """ 
+        Validate the contents of the rc-file given a dictionary of required keys
+        """
+
+        needed_rc_items = ['obs.input.dir',
+                           'obs.input.fname',
+                           'ocn.covariance',
+                           'nparameters',
+                           'bio.covariance',
+                           'deltaco2.prefix',
+                           'regtype']
+
+
+        for k,v in self.iteritems():
+            if v == 'True' : self[k] = True
+            if v == 'False': self[k] = False
+
+        for key in needed_rc_items:
+
+            if not self.has_key(key):
+                status,msg = ( False,'Missing a required value in rc-file : %s' % key)
+                logging.error(msg)
+                raise IOError,msg
+
+        status,msg = ( True,'DA System Info settings have been validated succesfully' )  ; logging.debug(msg)
+
+        return None
+################### End Class CtDaSystem ###################
+
+
+if __name__ == "__main__":
+    pass

da/ctgridded/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
+import datetime
+from da.baseclasses.statevector import EnsembleMember, StateVector
+import numpy as np
+
+identifier = 'CarbonTracker Gridded Statevector '
+version    = '0.0'
+
+################### Begin Class CtEnsembleMember ###################
+
+class CTEnsembleMember(EnsembleMember):
+
+    def AddCustomFields(self,filehandle, DaCycle):
+        """ Placeholder to add more custom fields to the output of the ensemble member"""
+        import da.tools.io4 as io
+        #import da.tools.io as io
+        from da.ct.tools import StateToGrid
+        import numpy as np
+
+
+        data                    = StateToGrid(self.ParameterValues,DaCycle.DaSystem.regionmap)
+
+        dimgrid                 = filehandle.AddLatLonDim()
+
+        savedict                = io.std_savedict.copy()
+        savedict['name']        = "parametermap"
+        savedict['long_name']   = "parametermap_for_member_%d"%self.membernumber
+        savedict['units']       = "unitless"
+        savedict['dims']        = dimgrid 
+        savedict['values']      = data.tolist()
+        savedict['comment']     = 'These are gridded parameter values to use for member %d'%self.membernumber
+        dummy                   = filehandle.AddData(savedict)
+
+        return None
+
+################### Begin Class CtStateVector ###################
+
+class CtGriddedStateVector(StateVector):
+    """ This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """
+
+    def getid(self):
+        return identifier
+
+    def getversion(self):
+        return version
+
+    def GetCovariance(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
+
+        file_ocn_cov = self.DaCycle.DaSystem['ocn.covariance'] 
+        file_bio_cov = self.DaCycle.DaSystem['bio.covariance'] 
+
+        # replace YYYY.MM in the ocean covariance file string
+
+        file_ocn_cov = file_ocn_cov.replace('2000.01',date.strftime('%Y.%m'))
+
+        for file in [file_ocn_cov,file_bio_cov]:
+
+            if not os.path.exists(file):
+
+                msg = "Cannot find the specified file %s" % file ; logging.error(msg)
+                raise IOError,msg
+            else:
+
+                msg = "Using covariance file: %s" % file ; logging.info(msg)
+
+        f_ocn   = io.CT_Read(file_ocn_cov,'read')
+        f_bio   = io.CT_Read(file_bio_cov,'read')
+
+        cov_ocn = f_ocn.GetVariable('CORMAT')
+        cov_bio = f_bio.GetVariable('qprior')
+
+        dummy   = f_ocn.close()
+        dummy   = f_bio.close()
+
+        dummy   = logging.debug("Succesfully closed files after retrieving prior covariance matrices")
+
+        # Once we have the matrices, we can start to make the full covariance matrix, and then decompose it
+
+        fullcov = np.zeros((self.nparams,self.nparams),float)
+
+        nocn    = cov_ocn.shape[0]
+        nbio    = cov_bio.shape[0]
+
+        fullcov[0:nbio,0:nbio]                  = cov_bio
+        fullcov[nbio:nbio+nocn,nbio:nbio+nocn]  = 0.16 *np.dot(cov_ocn,cov_ocn)
+        fullcov[nocn+nbio,nocn+nbio]            = 1.e-10
+
+
+        try:
+            plt.imshow(fullcov)
+            plt.colorbar()
+            plt.savefig('fullcovariancematrix.png')
+            plt.close('all')
+            dummy   = logging.debug("Covariance matrix visualized for inspection")
+        except:
+            pass
+
+        return fullcov
+
+    def MakeNewEnsemble(self,lag, covariancematrixlist = [None]):
+        """ 
+        :param lag: an integer indicating the time step in the lag order
+        :param covariancematrix: a list of matrices specifying the covariance distribution to draw from
+        :rtype: None
+    
+        Make a new ensemble, 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]
+
+        The covariance list object to be passed holds a list of matrices with a total number of dimensions [nparams, nparams], which is
+        used to draw ensemblemembers from. Each draw is done on a matrix from the list, to make the computational burden smaller when
+        the StateVector nparams becomes very large.
+
+        """    
+
+        if not isistance(covariancematrixlist,list):
+            msg = "The covariance matrix or matrices must be passed as a list of array objects, exiting..." ; logging.error(msg)
+            raise ValueError
+
+        # Check dimensions of covariance matrix list, must add up to nparams
+
+        dims = 0
+        for matrix in covariancematrixlist: dims += matrix.shape[0]
+
+        if dims != self.nparams:
+            msg = "The total dimension of the covariance matrices passed (%d) does not add up to the prescribed nparams (%d), exiting..."%(dims,self.nparams) ; logging.error(msg)
+            raise ValueError
+
+        # Loop over list if identity matrices and create a matrix of (nparams,nmembers) with the deviations
+
+        all_devs = []
+        for matrix in covariancematrixlist:
+
+            # Make a cholesky decomposition of the covariance matrix
+
+            U,s,Vh  = np.linalg.svd(covariancematrix)
+            dof     = np.sum(s)**2/sum(s**2)
+            C       = np.linalg.cholesky(covariancematrix)
+
+            dummy   = logging.debug("Covariance matrix visualized for inspection")
+
+            msg =   'Cholesky decomposition has succeeded offline'                    ; logging.debug(msg)
+            msg =   'Appr. degrees of freedom in covariance matrix is %s'%(int(dof))  ; logging.info(msg)
+
+            # Draw nmembers instances of this distribution
+
+            npoints      = matrix.shape[0]
+            randstate    = np.random.get_state()
+            for member in range(1,self.nmembers):
+                rands        = np.random.randn(npoints)
+                deviations   = np.dot(C,rands)
+
+                all_devs.append(deviations)
+
+            sys.exit(2)
+
+            # Create mean values 
+
+            NewMean                     = np.ones(self.nparams,float) # standard value for a new time step is 1.0
+
+            # If this is not the start of the filter, average previous two optimized steps into the mix
+
+            if lag == self.nlag and self.nlag >= 3:
+                NewMean                     += self.EnsembleMembers[lag-2][0].ParameterValues + \
+                                               self.EnsembleMembers[lag-3][0].ParameterValues 
+                NewMean                      = NewMean/3.0
+
+            # Create the first ensemble member with a deviation of 0.0 and add to list
+
+            NewMember                   = self.GetNewMember(0)
+            NewMember.ParameterValues   = NewMean.flatten()  # no deviations
+            dummy                       = self.EnsembleMembers[lag-1].append(NewMember)
+
+            # Create members 1:nmembers and add to EnsembleMembers list
+
+            for member in range(1,self.nmembers):
+                randstate                   = np.random.get_state()
+                rands                       = np.random.randn(self.nparams)
+
+                NewMember                   = self.GetNewMember(member)
+                NewMember.ParameterValues   = np.dot(C,rands) + NewMean
+                dummy                       = self.EnsembleMembers[lag-1].append(NewMember)
+
+            msg =   '%d new ensemble members were added to the state vector # %d'%(self.nmembers,lag)  ; logging.debug(msg)
+
+
+    def GetNewMember(self,memberno):
+        """ Return an CTensemblemember object """
+
+        return CTEnsembleMember(memberno)
+
+
+################### End Class CtStateVector ###################
+
+
+if __name__ == "__main__":
+
+    import os
+    import sys
+    from da.tools.initexit import StartLogger 
+    from da.tools.pipeline import JobStart
+    import numpy as np
+
+    sys.path.append(os.getcwd())
+
+    opts = ['-v']
+    args = {'rc':'da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'}
+
+    StartLogger()
+
+    DaCycle = JobStart(opts,args)
+
+    DaCycle.Initialize()
+
+    StateVector = CtStateVector()
+
+    StateVector.Initialize()
+
+    for n in range(dims[0]):
+        cov   = StateVector.GetCovariance()
+        dummy = StateVector.MakeNewEnsemble(n+1,cov)
+
+    StateVector.Propagate()
+
+    savedir         = DaCycle['dir.output']
+    filename        = os.path.join(savedir,'savestate.nc')
+
+    dummy = StateVector.WriteToFile()
+
+    StateVector.ReadFromFile(filename)
+