Commit 046daa97 authored by Peters, Wouter's avatar Peters, Wouter
Browse files

base class updated to initiate a proper ensemble at initialization time

parent 2fc93b50
......@@ -15,13 +15,14 @@ import logging
import datetime
################### Begin Class CtMember ###################
################### Begin Class EnsembleMember ###################
class CtMember(object):
class EnsembleMember(object):
"""
An ensemble member for CarbonTracker. This consists of
An ensemble member. This consists of
* a member number
* parameter values
* file names to write/read the data for this member
* a ModelSample object to hold sampled values for this member
"""
def __init__(self, membernumber):
......@@ -54,119 +55,50 @@ class CtMember(object):
dummy = f.close()
msg = 'Successfully wrote data from ensemble member %d to file (%s) ' % (self.membernumber,filename,) ; logging.debug(msg)
msg = 'Successfully wrote data from ensemble member %d to file (%s) ' % (self.membernumber,filename,) ; logging.info(msg)
################### End Class CtMember ###################
################### End Class EnsembleMember ###################
################### Begin Class CtStateVector ###################
################### Begin Class StateVector ###################
class CtStateVector(object):
class StateVector(object):
""" an object that holds data + methods and attributes needed to manipulate state vector values """
def __init__(self,CycleInfo):
def __init__(self,dims):
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.nlag = dims[0]
self.nmembers = dims[1]
self.nparams = dims[2]
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, we define member 0 as the mean of the distribution and n=1,...,nmembers as the spread.
# 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] = []
dummy = self.MakeNewEnsemble(n+1)
def __str__(self):
return "This is a CarbonTracker state vector object"
return "This is a base class derived state vector object"
def MakeNewEnsemble(self,lag):
""" Make a new ensemble from specified matrices, the attribute lag refers to the position in the state vector.
""" 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]
"""
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']
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 = Nio.open_file(file_ocn_cov,'r')
f_bio = Nio.open_file(file_bio_cov,'r')
cov_ocn = f_ocn.variables['CORMAT'].get_value()
cov_bio = f_bio.variables['qprior'].get_value()
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
# Visually inspect full covariance matrix if verbose
if self.CycleInfo.da_settings['verbose']:
plt.imshow(fullcov)
plt.colorbar()
plt.savefig('fullcovariancematrix.png')
plt.close('all')
# 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)
# Create mean values
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 = EnsembleMember(0)
NewMember.ParameterValues = np.zeros((self.nparams),float) + NewMean
dummy = self.EnsembleMembers[lag-1].append(NewMember)
......@@ -176,8 +108,8 @@ class CtStateVector(object):
randstate = np.random.get_state()
rands = np.random.randn(self.nparams)
NewMember = CtMember(member)
NewMember.ParameterValues = np.dot(C,rands) + NewMean
NewMember = EnsembleMember(member)
NewMember.ParameterValues = rands + NewMean # Note that the mean is added to the param values for each member
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)
......@@ -187,7 +119,7 @@ class CtStateVector(object):
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
In the future, this routine can incorporate a formal propagation of the statevector.
"""
......@@ -206,24 +138,16 @@ class CtStateVector(object):
return None
def WriteToFile(self):
def WriteToFile(self,filename):
""" 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]
......@@ -262,73 +186,13 @@ class CtStateVector(object):
for n in range(self.nlag):
for m in range(self.nmembers):
NewMember = CtMember(m)
NewMember = EnsembleMember(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)
################### End Class CtStateVector ###################
################# Stand alone methods that manipulate the objects above ##############
def PrepareState(CycleInfo):
"""
Prepare the ensemble of parameters needed by the forecast model. There are two possible pathways:
(A) Construct a brand new ensemble from a covariance matrix
(B) Get an existing ensemble and propagate it
The steps for procedure A are:
(A1) Construct the covariance matrix
(A2) Make a Cholesky decomposition
(A3) Prepare a set of [nmembers] input parameters for each forecast member
The steps for procedure B are:
(B1) Get existing parameter values for each member
(B2) Run them through a forecast model
Both procedures finish with writing the parameter values to a NetCDF file
"""
StateVector = CtStateVector(CycleInfo)
nlag = int(CycleInfo.da_settings['time.nlag'])
# 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
# with new values for each week. After we have constructed the StateVector, it will be propagated by one cycle length so it is ready to be used
# in the current cycle
if not CycleInfo.da_settings['time.restart']:
# Fill each week from n=1 to n=nlag with a new ensemble
for n in range(0,nlag):
dummy = StateVector.MakeNewEnsemble(n+1)
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')
filename = os.path.join(savedir,'savestate.nc')
StateVector.ReadFromFile(filename)
# Now propagate the ensemble by one cycle to prepare for the current cycle
dummy = StateVector.Propagate()
return StateVector
################### End Class StateVector ###################
if __name__ == "__main__":
......@@ -347,23 +211,32 @@ if __name__ == "__main__":
StartLogger()
DaCycle = CycleControl(opts,args)
dummy = CtMember(0)
dummy = EnsembleMember(0)
print dummy
DaCycle.Initialize()
print DaCycle
StateVector = PrepareState(DaCycle)
dims = ( int(DaCycle.da_settings['time.nlag']),
int(DaCycle.da_settings['forecast.nmembers']),
int(DaCycle.DaSystem.da_settings['nparameters']),
)
StateVector = StateVector(dims)
members = StateVector.EnsembleMembers[1]
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)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment