Commit 9ae301fa authored by Peters, Wouter's avatar Peters, Wouter
Browse files

changed back init and call to MakeNewEnsemble to include a covariance matrix

parent 80a11873
......@@ -63,6 +63,7 @@ class EnsembleMember(object):
################### Begin Class StateVector ###################
import numpy as np
class StateVector(object):
""" an object that holds data + methods and attributes needed to manipulate state vector values """
def __init__(self,dims):
......@@ -79,18 +80,34 @@ class StateVector(object):
for n in range(self.nlag):
self.EnsembleMembers[n] = []
dummy = self.MakeNewEnsemble(n+1)
def __str__(self):
return "This is a base class derived state vector object"
def MakeNewEnsemble(self,lag):
def MakeNewEnsemble(self,lag, covariancematrix = 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 optional DaCycle object to be passed holds info on the data assimilation system settings needed
for instance to read covariance data from file. This is not part of the base class object here, hence, this
mehtod will almost always be overwritten in the derived class for the specific application.
The code below is just an example of what could be used
"""
import numpy as np
if covariancematrix == None:
covariancematrix=np.identity(self.nparams)
# 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)
msg = 'Cholesky decomposition has succeeded offline' ; logging.debug(msg)
msg = 'Appr. degrees of freedom in covariance matrix is %s'%(int(dof)) ; logging.info(msg)
# Create mean values
......@@ -109,7 +126,7 @@ class StateVector(object):
rands = np.random.randn(self.nparams)
NewMember = EnsembleMember(member)
NewMember.ParameterValues = rands + NewMean # Note that the mean is added to the param values for each 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)
......@@ -224,6 +241,9 @@ if __name__ == "__main__":
StateVector = StateVector(dims)
StateVector.MakeNewEnsemble(lag=1)
StateVector.MakeNewEnsemble(lag=2)
members = StateVector.EnsembleMembers[1]
members[0].WriteToFile(DaCycle.da_settings['dir.input'])
......
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