updates from Emma's code

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 ###################