diff --git a/da/ct/optimizer.py b/da/ct/optimizer.py index 9bc1cc299db85e4a5a2525bb65d91f708cac9942..53f31350792732dc1603ce4aea2be464fbdc2f55 100755 --- a/da/ct/optimizer.py +++ b/da/ct/optimizer.py @@ -21,7 +21,7 @@ identifier = 'CarbonTracker CO2' class CtOptimizer(Optimizer): """ - This creates an instance of a CarbonTracker optimization object. The base class it derives from is the optimezer object. + This creates an instance of a CarbonTracker optimization object. The base class it derives from is the optimizer object. Additionally, this CtOptimizer implements a special localization option following the CT2007 method. All other methods are inherited from the base class Optimizer. diff --git a/da/ct/statevector.py b/da/ct/statevector.py index 4b58ee961c95425f2a1c72fa0b0dda0b2a79a88c..08eaa9ad3ab9db69cc847d87a91bca108ce776c6 100755 --- a/da/ct/statevector.py +++ b/da/ct/statevector.py @@ -11,93 +11,31 @@ File created on 28 Jul 2010. import os import sys +sys.path.append('../../') + import logging import datetime +from da.baseclasses.statevector import EnsembleMember, StateVector +import numpy as np -################### Begin Class CtMember ################### - -class CtMember(object): - """ - An ensemble member for CarbonTracker. This consists of - * parameter values - * file names to write/read the data for this member - """ - - def __init__(self, membernumber): - self.membernumber = membernumber # the member number - self.ParameterValues = None # Parameter values of this member - self.ModelSample = None # Model Sampled Parameter values of this member - - 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 da.tools.io import CT_CDF, std_savedict - import numpy as np - - filename = os.path.join(outdir,'parameters.%03d.nc'% self.membernumber) - f = CT_CDF(filename,'create') - dimparams = f.AddParamsDim(len(self.ParameterValues)) - - data = self.ParameterValues - - 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.debug(msg) - - - -################### End Class CtMember ################### - ################### Begin Class CtStateVector ################### -class CtStateVector(object): - """ an object that holds data + methods and attributes needed to manipulate state vector values """ - def __init__(self,CycleInfo): - - self.CycleInfo = CycleInfo - self.DaSystem = CycleInfo.DaSystem - self.nlag = int(self.CycleInfo.da_settings['time.nlag']) - self.nmembers = int(self.CycleInfo.da_settings['forecast.nmembers']) - self.nparams = int(self.DaSystem.da_settings['nparameters']) - self.isOptimized = False +class CtStateVector(StateVector): + """ This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """ - # 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, 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] = [] - - def __str__(self): - return "This is a CarbonTracker state vector object" - - def MakeNewEnsemble(self,lag): + def GetCovariance(self,DaSystem): """ 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 Nio - import numpy as np - import matplotlib.pyplot as plt # Get the needed matrices from the specified covariance files - file_ocn_cov = self.DaSystem.da_settings['ocn.covariance'] - file_bio_cov = self.DaSystem.da_settings['bio.covariance'] + file_ocn_cov = DaSystem.da_settings['ocn.covariance'] + file_bio_cov = DaSystem.da_settings['bio.covariance'] for file in [file_ocn_cov,file_bio_cov]: @@ -133,140 +71,15 @@ class CtStateVector(object): # Visually inspect full covariance matrix if verbose - if self.CycleInfo.da_settings['verbose']: - + try: plt.imshow(fullcov) plt.colorbar() plt.savefig('fullcovariancematrix.png') plt.close('all') + except: + pass - # Make a cholesky decomposition of the covariance matrix - - U,s,Vh = np.linalg.svd(fullcov) - dof = np.sum(s)**2/sum(s**2) - - C = np.linalg.cholesky(cov_bio) - C = np.linalg.cholesky(cov_ocn) - print fullcov.max(),fullcov.min() - C = np.linalg.cholesky(fullcov) - - msg = 'Cholesky decomposition has succeeded offline' ; logging.debug(msg) - msg = 'Appr. degrees of freedom in covariance matrix is %s'%(int(dof)) ; logging.info(msg) - - if self.CycleInfo.da_settings['verbose']: - - plt.plot(s) - plt.savefig('eigenvaluespectrum.png') - plt.close('all') - - - # Create mean values for the new time step (default = 1.0) - - NewMean = np.ones(self.nparams,float) - - # Create the first ensemble member with a deviation of 0.0 and add to list - - NewMember = CtMember(0) - NewMember.ParameterValues = np.zeros((self.nparams),float) + NewMean - 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 = CtMember(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 Propagate(self): - """ - - 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 using for instance a biosphere model - - """ - - # Remove State Vector n=1 by simply "popping" it from the list and appending a new empty list at the front. This empty list will - # hold the new ensemble for the new cycle - - dummy = self.EnsembleMembers.pop(0) - - dummy = self.EnsembleMembers.append([]) - - # And now create a new week of mean + members for n=nlag - - dummy = self.MakeNewEnsemble(self.nlag) - - msg = 'The state vector has been propagated by one cycle ' ; logging.info(msg) - - return None - - def WriteToFile(self): - """ Write the StateVector object to a netcdf file for future use """ - from da.tools.io import CT_CDF - import numpy as np - - outdir = self.CycleInfo.da_settings['dir.output'] - ifiltertime = self.CycleInfo.da_settings['time.start'].strftime('%Y%m%d') - efiltertime = self.CycleInfo.da_settings['time.end'].strftime('%Y%m%d') - filename = os.path.join(outdir,'savestate.nc' ) - f = CT_CDF(filename,'create') - dimparams = f.AddParamsDim(self.nparams) - dimmembers = f.AddMembersDim(self.nmembers) - dimlag = f.AddLagDim(self.nlag,unlimited=True) - - setattr(f,'cyclestartdate',ifiltertime) - setattr(f,'cycleenddate',efiltertime) - setattr(f,'isOptimized',str(self.isOptimized)) - - for n in range(self.nlag): - - members = self.EnsembleMembers[n] - MeanState = members[0].ParameterValues - - savedict = f.StandardVar(varname='meanstate') - savedict['dims'] = dimlag+dimparams - savedict['values'] = MeanState - savedict['count'] = n - savedict['comment'] = 'this represents the mean of the ensemble' - dummy = f.AddData(savedict) - - members = self.EnsembleMembers[n] - data = np.array([m.ParameterValues for m in members])- MeanState - - savedict = f.StandardVar(varname='ensemblestate') - savedict['dims'] = dimlag+dimmembers+dimparams - savedict['values'] = data - 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,filename): - """ Read the StateVector object from a netcdf file for future use """ - import Nio - import numpy as np - - f = Nio.open_file(filename,'r') - MeanState = f.variables['statevectormean'].get_value() - EnsembleMembers = f.variables['statevectorensemble'].get_value() - dummy = f.close() - - for n in range(self.nlag): - for m in range(self.nmembers): - NewMember = CtMember(m) - NewMember.ParameterValues = EnsembleMembers[n,m,:] + MeanState # add the mean to the deviations to hold the full parameter values - dummy = self.EnsembleMembers[n].append(NewMember) - - msg = 'Successfully read the State Vector from file (%s) ' % (filename,) ; logging.info(msg) + return fullcov ################### End Class CtStateVector ################### @@ -297,10 +110,13 @@ def PrepareState(CycleInfo): """ + dims = ( int(CycleInfo.da_settings['time.nlag']), + int(CycleInfo.da_settings['forecast.nmembers']), + int(CycleInfo.DaSystem.da_settings['nparameters']), + ) + StateVector = CtStateVector(dims) - StateVector = CtStateVector(CycleInfo) - - nlag = int(CycleInfo.da_settings['time.nlag']) + nlag = dims[0] # We now have an empty CtStateVector object that we need to populate with data. If this is a continuation from a previous cycle, we can read # the previous StateVector values from a NetCDF file in the save/ directory. If this is the first cycle, we need to populate the CtStateVector @@ -312,13 +128,14 @@ def PrepareState(CycleInfo): # Fill each week from n=1 to n=nlag with a new ensemble for n in range(0,nlag): - dummy = StateVector.MakeNewEnsemble(n+1) + cov = StateVector.GetCovariance(CycleInfo.DaSystem) + dummy = StateVector.MakeNewEnsemble(n+1,cov) else: # Read the StateVector data from file - savedir = StateVector.CycleInfo.da_settings['dir.save'] - filtertime = StateVector.CycleInfo.da_settings['time.start'].strftime('%Y%m%d') + savedir = CycleInfo.da_settings['dir.save'] + filtertime = CycleInfo.da_settings['time.start'].strftime('%Y%m%d') filename = os.path.join(savedir,'savestate.nc') StateVector.ReadFromFile(filename) @@ -332,8 +149,6 @@ def PrepareState(CycleInfo): if __name__ == "__main__": - sys.path.append('../../') - import os import sys from da.tools.general import StartLogger @@ -347,7 +162,7 @@ if __name__ == "__main__": StartLogger() DaCycle = CycleControl(opts,args) - dummy = CtMember(0) + dummy = EnsembleMember(0) print dummy DaCycle.Initialize() @@ -359,11 +174,15 @@ if __name__ == "__main__": members[0].WriteToFile(DaCycle.da_settings['dir.input']) - data = np.array([m.ParameterValues for m in members]) + StateVector.Propagate() + + savedir = DaCycle.da_settings['dir.output'] + filename = os.path.join(savedir,'savestate.nc') + + StateVector.WriteToFile(filename) - StateVector.WriteToFile() - savedir = StateVector.CycleInfo.da_settings['dir.output'] + savedir = DaCycle.da_settings['dir.output'] filename = os.path.join(savedir,'savestate.nc') StateVector.ReadFromFile(filename)