diff --git a/da/sf6/statevector.py b/da/sf6/statevector.py
index 3819a51db6c02ac3c8a061f6284e25e3c7561799..8f83f33eecf827230e46ba9849730711c5a7b88f 100755
--- a/da/sf6/statevector.py
+++ b/da/sf6/statevector.py
@@ -41,7 +41,8 @@ class SF6StateVector(StateVector):
fullcov = np.zeros((self.nparams, self.nparams), float)
- fullcov[:,:] = 0.2**2
+ for n in range(self.nparams):
+ fullcov[n,n] = 0.00001**2
return fullcov
@@ -141,16 +142,134 @@ class SF6StateVector(StateVector):
"""
- ensemble_deviations = np.array([p.param_values for p in self.ensemble_members[0] ])
- if ensemble_deviations.std(axis=0) < 0.05 :
- for m in self.ensemble_members[0]:
- m = 3.0 * m # inflate deviations by a factor of 3.0
-
- logging.info('The state vector covariance was inflated to 1-sigma of %5.2f ' % (ensemble_deviations.std(axis=0)*3.0))
+ self.ensemble_members.append([])
+ cov = self.get_covariance(None)
+ newmean = self.ensemble_members[self.nlag-1][0].param_values
+# newmean = np.ones(self.nparams)*0.1 #VDV
+ self.make_new_ensemble(self.nlag+1, newmean = newmean, covariancematrix=cov)
+ self.ensemble_members.pop(0)
logging.info('The state vector remains the same in the SF6 run')
logging.info('The state vector has been propagated by one cycle')
+ def make_new_ensemble(self, lag, newmean=None, covariancematrix=None):
+ """
+ :param lag: an integer indicating the time step in the lag order
+ :param covariancematrix: a matrix to draw random values from
+ :rtype: 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 covariance object to be passed holds a matrix of dimensions [nparams, nparams] which is
+ used to draw ensemblemembers from. If this argument is not passed it will ne substituted with an
+ identity matrix of the same dimensions.
+
+ """
+
+ if newmean == None:
+ newmean = np.ones(self.nparams)*0.1 # emma 05-02
+
+ if covariancematrix == None:
+ covariancematrix = np.identity(self.nparams)
+
+ # Make a cholesky decomposition of the covariance matrix
+
+
+ _, s, _ = np.linalg.svd(covariancematrix)
+ dof = np.sum(s) ** 2 / sum(s ** 2)
+ C = np.linalg.cholesky(covariancematrix)
+
+ logging.debug('Cholesky decomposition has succeeded ')
+ logging.info('Appr. degrees of freedom in covariance matrix is %s' % (int(dof)))
+
+
+ # Create the first ensemble member with a deviation of 0.0 and add to list
+
+ newmember = EnsembleMember(0)
+ newmember.param_values = newmean.flatten() # no deviations
+ self.ensemble_member[lag-1].append(newmember)
+
+ # Create members 1:nmembers and add to EnsembleMembers list
+
+ for member in range(1, self.nmembers):
+ rands = np.random.randn(self.nparams)
+
+ newmember = EnsembleMember(member)
+ # routine to avoids that members < 0.0
+ dev = np.dot(C, rands) #VDV
+ dummy = np.zeros(len(dev)) #VDV
+ for i in range(len(dev)): #VDV
+ if dev[i] < 0.0: #VDV
+ dummy[i] = newmean[i]*np.exp(dev[i]) #VDV
+ else: #VDV
+ dummy[i] = newmean[i]*(1+dev[i]) #VDV
+ newmember.param_values = dummy #VDV
+ #newmember.ParameterValues = np.dot(C, rands) + newmean
+ self.ensemble_members[lag-1].append(newmember)
+
+ logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, lag))
+
+ def write_members_to_file(self, lag, outdir):
+ """
+ :param: lag: Which lag step of the filter to write, must lie in range [1,...,nlag]
+ :rtype: None
+
+ Write ensemble member information to a NetCDF file for later use. The standard output filename is
+ *parameters.DDD.nc* where *DDD* is the number of the ensemble member. Standard output file location
+ is the `dir.input` of the dacycle object. In principle the output file will have only two datasets inside
+ called `parametervalues` which is of dimensions `nparameters` and `parametermap` which is of dimensions (180,360).
+ This dataset can be read and used by a :class:`~da.baseclasses.observationoperator.ObservationOperator` object.
+
+ .. note:: if more, or other information is needed to complete the sampling of the ObservationOperator you
+ can simply inherit from the StateVector baseclass and overwrite this write_members_to_file function.
+
+ """
+
+ # These import statements caused a crash in netCDF4 on MacOSX. No problems on Jet though. Solution was
+ # to do the import already at the start of the module, not just in this method.
+
+ #import da.tools.io as io
+ #import da.tools.io4 as io
+
+ members = self.ensemble_members[lag]
+
+ for mem in members:
+ filename = os.path.join(outdir, 'parameters.%03d.nc' % mem.membernumber)
+ ncf = io.CT_CDF(filename, method='create')
+ dimparams = ncf.add_params_dim(self.nparams)
+ dimgrid = ncf.add_latlon_dim()
+ logging.warning('Parameters for TM5 are NOT taken with exponential')
+
+ data = mem.param_values
+
+ savedict = io.std_savedict.copy()
+ savedict['name'] = "parametervalues"
+ savedict['long_name'] = "parameter_values_for_member_%d" % mem.membernumber
+ savedict['units'] = "unitless"
+ savedict['dims'] = dimparams
+ savedict['values'] = data
+ savedict['comment'] = 'These are parameter values to use for member %d, note: they result from an exponential function' % mem.membernumber
+ ncf.add_data(savedict)
+
+ griddata = self.vector2grid(vectordata=data)
+
+ savedict = io.std_savedict.copy()
+ savedict['name'] = "parametermap"
+ savedict['long_name'] = "parametermap_for_member_%d" % mem.membernumber
+ savedict['units'] = "unitless"
+ savedict['dims'] = dimgrid
+ savedict['values'] = griddata.tolist()
+ savedict['comment'] = 'These are gridded parameter values to use for member %d, note: they result from an exponential function' % mem.membernumber
+ ncf.add_data(savedict)
+
+ ncf.close()
+
+ logging.debug('Successfully wrote data from ensemble member %d to file (%s) ' % (mem.membernumber, filename))
+
+
+
################### End Class SF6StateVector ###################