Commit 80cad17e authored by Peters, Wouter's avatar Peters, Wouter
Browse files

parameter numbers used to fill statevector from cov files

parent 46728a88
......@@ -27,6 +27,16 @@ version = '0.0'
class CO2GriddedStateVector(StateVector):
""" This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """
def __init__(self):
self.ID = identifier
self.version = version
# The following code allows the object to be initialized with a dacycle object already present. Otherwise, it can
# be added at a later moment.
logging.info('Statevector object initialized: %s' % self.ID)
def get_covariance(self, date, dacycle):
""" 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.
......@@ -67,16 +77,19 @@ class CO2GriddedStateVector(StateVector):
if 'pco2' in file or 'cov_ocean' in file:
cov_ocn = f.get_variable('CORMAT')
parnr = range(9805,9835)
cov = cov_ocn
elif 'tropic' in file:
cov = f.get_variable('covariance')
parnr = f.get_variable('parameternumber')
else:
cov = f.get_variable('covariance')
parnr = f.get_variable('parameternumber')
cov_sf = 90.0 / np.sqrt(cov.diagonal().sum()) # this scaling factor makes the total variance close to the value of a single ecoregion
cov = cov * cov_sf
f.close()
covariancematrixlist.append(cov)
covariancematrixlist.append([cov,parnr,file])
logging.debug("Succesfully closed files after retrieving prior covariance matrices")
......@@ -112,7 +125,7 @@ class CO2GriddedStateVector(StateVector):
dims = 1 # start from 1.0 to account for the last parameter that scales Ice+Non-optimized, we have no covariance matrix for this though
for matrix in covariancematrixlist:
for matrix,parnr,file in covariancematrixlist:
dims += matrix.shape[0]
if dims != self.nparams:
......@@ -121,13 +134,11 @@ class CO2GriddedStateVector(StateVector):
# Loop over list if identity matrices and create a matrix of (nparams,nmembers) with the deviations
istart = 0
istop = 0
dof = 0.0
dev_matrix = np.zeros((self.nparams, self.nmembers,), 'float')
randstate = np.random.get_state()
for matrix in covariancematrixlist:
for matrix,parnr,file in covariancematrixlist:
# Make a cholesky decomposition of the covariance matrix
_, s, _ = np.linalg.svd(matrix)
......@@ -145,18 +156,21 @@ class CO2GriddedStateVector(StateVector):
npoints = matrix.shape[0]
istop = istop + npoints
istart = parnr[0]-1
istop = parnr[-1]
for member in range(1, self.nmembers):
rands = np.random.randn(npoints)
deviations = np.dot(C, rands)
dev_matrix[istart:istop, member - 1] = deviations
dev_matrix[istop, member - 1] = 1.e-10 * np.random.randn()
dev_matrix[-1, member - 1] = 1.e-10 * np.random.randn()
logging.debug("Inserted prior covariance matrix from file %s of %d elements into statevector location %d through %d"%(file,len(parnr),istart,istop) )
#cov2 = np.dot(dev_matrix[istart:istop,:],np.transpose(dev_matrix[istart:istop,:])) / (self.nmembers-1)
#print matrix.sum(),cov2.sum(),abs(matrix.diagonal()-cov2.diagonal()).max(), matrix.shape,cov2.shape
istart = istart + npoints
logging.debug('Successfully constructed a deviation matrix from covariance structure')
logging.info('Appr. degrees of freedom in full covariance matrix is %s' % (int(dof)))
......
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