statevectors and ensemble members can now write data to file

ct_tools.py

@@ -243,13 +243,40 @@ class CtMember():
         self.membernumber           = membernumber   # the member number
         self.randomstate            = (None,)        # the state of the random number generator needed to reproduce the random sequence
         self.ModelSamples           = None           # CtObservations object that holds model samples for this member
-        self.ParameterValues        = None           # Parameter values of this member
+        self.MeanDeviations         = None           # Parameter values of this member
         self.ParameterOutputFile    = None           # File to which the member info should be written
         self.SampleInputFile        = None           # File from which the sample info can be read
 
     def __str__(self):
         return "%03d"%self.membernumber
 
+    def WriteToFile(self, outdir):
+        """ Write the information needed by an external model to a netcdf file for future use """
+        from ct_netcdf import CT_CDF, std_savedict
+        import pycdf as cdf
+        import mysettings
+        import numpy as np
+
+        filename        = os.path.join(outdir,'parameters.%03d.nc'% self.membernumber)
+        f               = CT_CDF(filename,'create')
+        dimparams       = f.AddParamsDim(len(self.MeanDeviations))
+
+        data =  self.MeanDeviations
+
+        savedict                = std_savedict.copy()
+        savedict['name']        = "parametervalues"
+        savedict['long_name']   = "parameter_values_for_member_%d"%self.membernumber
+        savedict['units']       = "unitless"
+        savedict['dims']        = dimparams 
+        savedict['values']      = data.tolist()
+        savedict['comment']     = 'These are parameter values to use for member %d'%self.membernumber
+        dummy                   = f.AddData(savedict)
+
+        dummy           = f.close()
+
+        msg     = 'Successfully wrote data from ensemble member %d to file (%s) ' % (self.membernumber,filename,) ; logging.info(msg)
+
+
 
 ################### End Class CtMember ###################
 
@@ -266,7 +293,12 @@ class CtStateVector():
         self.nparams            = int(DaInfo.da_settings['nparameters'])
         self.isOptimized        = False
 
+        # 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 CtMember objects, whereas the MeanValues will simply be a list of array values.
+
+        self.MeanValues         = range(self.nlag)
         self.EnsembleMembers    = range(self.nlag)
+
         for n in range(self.nlag):
             self.EnsembleMembers[n] = []
 
@@ -345,10 +377,15 @@ class CtStateVector():
 
         dummy = np.random.seed()
 
-        # mean does not hold data yet
+        # Create mean values for the new time step (default = 1.0) and add to the MeanValues list
+
+        NewMean                     = np.ones(self.nparams,float)
+        self.MeanValues[nlag]       = NewMean
+
+        # Create the first ensemble member with a deviation of 0.0 and add to list
 
         NewMember                   = CtMember(0)
-        NewMember.ParameterValues   = np.ones((self.nparams,self.nmembers),float)  
+        NewMember.MeanDeviations    = np.zeros((self.nparams),float)  
         dummy                       = self.EnsembleMembers[nlag].append(NewMember)
 
         # Create members 1:nmembers and add to EnsembleMembers list
@@ -358,12 +395,80 @@ class CtStateVector():
             rands                       = np.random.randn(self.nparams)
 
             NewMember                   = CtMember(member)
-            NewMember.ParameterValues   = np.dot(C,rands)
+            NewMember.MeanDeviations    = np.dot(C,rands)
             NewMember.randomstate       = randstate
             dummy                       = self.EnsembleMembers[nlag].append(NewMember)
 
         msg =   '%d new ensemble members were added to the state vector'%(self.nmembers)  ; logging.info(msg)
 
+    def Propagate(self,nlag):
+
+        return self.MakeNewEnsemble(nlag)
+
+    def WriteToFile(self):
+        """ Write the StateVector object to a netcdf file for future use """
+        from ct_netcdf import CT_CDF
+        import pycdf as cdf
+        import mysettings
+        import numpy as np
+
+        outdir          = self.CycleInfo.da_settings['dir.output']
+        filtertime      = self.CycleInfo.da_settings['startdate'].strftime('%Y%m%d')
+        filename        = os.path.join(outdir,'savestate.%s.nc'% filtertime)
+        f               = CT_CDF(filename,'create')
+        dimparams       = f.AddParamsDim(self.nparams)
+        dimmembers      = f.AddMembersDim(self.nmembers)
+        dimlag          = f.AddLagDim(self.nlag,unlimited=True)
+
+        setattr(f,'date',filtertime)
+        setattr(f,'isOptimized',str(self.isOptimized))
+
+        for n in range(self.nlag):
+            data =  StateVector.MeanValues[n]
+
+            savedict                = f.StandardVar(varname='meanstate')
+            savedict['dims']        = dimlag+dimparams 
+            savedict['values']      = data.tolist()
+            savedict['count']       = n
+            savedict['comment']     = 'this represents the mean of the ensemble'
+            dummy                   = f.AddData(savedict)
+
+
+            members =  StateVector.EnsembleMembers[n]
+            data = np.array([m.MeanDeviations for m in members])
+
+            savedict                = f.StandardVar(varname='ensemblestate')
+            savedict['dims']        = dimlag+dimmembers+dimparams 
+            savedict['values']      = data.tolist()
+            savedict['count']       = n
+            savedict['comment']     = 'this represents deviations from the mean of the ensemble'
+            dummy                   = f.AddData(savedict)
+
+        dummy           = f.close()
+
+        msg     = 'Successfully wrote the State Vector to file (%s) ' % (filename,) ; logging.info(msg)
+
+    def ReadFromFile(self):
+        """ Read the StateVector object from a netcdf file for future use """
+        import pyhdf.SD as SD
+        import numpy as np
+
+        outdir          = self.CycleInfo.da_settings['dir.output']
+        filtertime      = self.CycleInfo.da_settings['startdate'].strftime('%Y%m%d')
+        filename        = os.path.join(outdir,'savestate.%s.nc'% filtertime)
+
+        f               = SD.SD(filename)
+        MeanState       = f.select('statevectormean').get()
+        EnsembleMembers = f.select('statevectorensemble').get()
+        dummy           = f.end()
+
+        for n in range(self.nlag):
+            self.MeanValues[n] = MeanState[n,:]
+            members            = self.EnsembleMembers[n]
+            for i,m in enumerate(members):
+                m.MeanDeviations = EnsembleMembers[n,i,:] 
+
+        msg     = 'Successfully read the State Vector from file (%s) ' % (filename,) ; logging.info(msg)
 
 ################### End Class CtStateVector ###################
 
@@ -391,11 +496,17 @@ class CtOptimizer():
         """ Create Matrix space needed in optimization routine """
         import numpy as np
 
-        self.X_prior         = np.array( (self.nlag*self.nparams,), 'float')
-        self.X_prior_prime   = np.array( (self.nmembers,self.nlag*self.nparams,), 'float')
+        self.X_prior         = np.zeros( (self.nlag*self.nparams,), float)
+        self.X_prior_prime   = np.zeros( (self.nmembers,self.nlag*self.nparams,), float)
 
     def FillMatrices(self,StateVector):
+        import numpy as np
 
+        for n in range(self.nlag):
+
+            self.X_prior[n*self.nparams:(n+1)*self.nparams]            = StateVector.MeanValues[n]
+            members                                                    = StateVector.EnsembleMembers[n]
+            self.X_prior_prime[:,n*self.nparams:(n+1)*self.nparams]    = np.array([m.MeanDeviations for m in members])
         
         return None
 
@@ -480,7 +591,7 @@ def PrepareState(CycleInfo):
 
             # Route B: Propagate existing ensemble
 
-            Xprime,Qprior = [0,0]
+            dummy = StateVector.Propagate(n-1)
 
     return StateVector
 
@@ -491,6 +602,7 @@ if __name__ == "__main__":
     import sys
     from tools_da import StartLogger 
     from da_initexit import CycleControl 
+    import numpy as np
 
     opts = ['-v']
     args = {'rc':'da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'}
@@ -508,7 +620,7 @@ if __name__ == "__main__":
     samples = CtObservations(Da_Info)
 
     dummy = samples.AddObs()
-    dummy = samples.AddSimulations()
+    dummy = samples.AddSimulations('/Users/peters/tmp/test_da/output/samples.000.nc')
 
     print samples.MrData.getvalues('code')
     print samples.MrData.getvalues('obs')
@@ -517,9 +629,22 @@ if __name__ == "__main__":
 
     dummy = PrepareObs(DaCycle,'forecast')
 
-    dummy = PrepareState(DaCycle)
+    StateVector = PrepareState(DaCycle)
+
 
+    members =  StateVector.EnsembleMembers[1]
 
+    members[0].WriteToFile(DaCycle.da_settings['dir.input'])
+
+    data = np.array([m.MeanDeviations for m in members])
+
+    StateVector.WriteToFile()
+
+    StateVector.ReadFromFile()
+
+    opt = CtOptimizer(DaCycle,Da_Info)
+    opt.FillMatrices(StateVector)
+    
 
 
 

da_initexit.py

@@ -16,6 +16,7 @@ needed_da_items=[
     'time.cycle',
     'dir.da_run',
     'forecast.model',
+    'forecast.nmembers',
     'forecast.model.rc',
     'da.system']
 
@@ -56,7 +57,6 @@ class CycleControl():
         self.da_settings['da.crash.recover'] = '-r' in opts
         self.da_settings['verbose']          = '-v' in opts
 
-
     def __str__(self):
         """
         String representation of a CycleControl object