statevector.py 44.4 KB
Newer Older
brunner's avatar
brunner committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters. 
Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
updates of the code. See also: http://www.carbontracker.eu. 

This program is free software: you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software Foundation, 
version 3. This program is distributed in the hope that it will be useful, but 
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

You should have received a copy of the GNU General Public License along with this 
program. If not, see <http://www.gnu.org/licenses/>."""
#!/usr/bin/env python
# ct_statevector_tools.py

"""
.. module:: statevector
.. moduleauthor:: Wouter Peters 

Revision History:
File created on 28 Jul 2010.

The module statevector implements the data structure and methods needed to work with state vectors (a set of unknown parameters to be optimized by a DA system) of different lengths, types, and configurations. Two baseclasses together form a generic framework:
    * :class:`~da.baseclasses.statevector.StateVector`
    * :class:`~da.baseclasses.statevector.EnsembleMember`

As usual, specific implementations of StateVector objects are done through inheritance form these baseclasses. An example of designing 
your own baseclass StateVector we refer to :ref:`tut_chapter5`.

.. autoclass:: da.baseclasses.statevector.StateVector 

.. autoclass:: da.baseclasses.statevector.EnsembleMember 

"""

import os
import logging
import numpy as np
from datetime import timedelta
import da.cosmo.io4 as io
from netCDF4 import Dataset

identifier = 'Baseclass Statevector '
version = '0.0'

################### Begin Class EnsembleMember ###################

class EnsembleMember(object):
    """ 
        An ensemble member object consists of:
           * a member number
           * parameter values
           * an observation object to hold sampled values for this member

        Ensemble members are initialized by passing only an ensemble member number, all data is added by methods 
        from the :class:`~da.baseclasses.statevector.StateVector`. Ensemble member objects have almost no functionality 
        except to write their data to file using method :meth:`~da.baseclasses.statevector.EnsembleMember.write_to_file`

        .. automethod:: da.baseclasses.statevector.EnsembleMember.__init__ 
        .. automethod:: da.baseclasses.statevector.EnsembleMember.write_to_file 
        .. automethod:: da.baseclasses.statevector.EnsembleMember.AddCustomFields 

    """

    def __init__(self, membernumber):
        """
           :param memberno: integer ensemble number
           :rtype: None

           An EnsembleMember object is initialized with only a number, and holds two attributes as containter for later
           data:
                * param_values, will hold the actual values of the parameters for this data
                * ModelSample, will hold an :class:`~da.baseclasses.obs.Observation` object and the model samples resulting from this members' data

        """
        self.membernumber = membernumber   # the member number
        self.param_values = None           # Parameter values of this member

################### End Class EnsembleMember ###################

################### Begin Class StateVector ###################


class StateVector(object):
    """ 
    The StateVector object first of all contains the data structure of a statevector, defined by 3 attributes that define the 
    dimensions of the problem in parameter space:
        * nlag
        * nparameters
        * nmembers

    The fourth important dimension `nobs` is not related to the StateVector directly but is initialized to 0, and later on 
    modified to be used in other parts of the pipeline:
        * nobs

    These values are set as soon as the :meth:`~da.baseclasses.statevector.StateVector.setup` is called from the :ref:`pipeline`. 
    Additionally, the value of attribute `isOptimized` is set to `False` indicating that the StateVector holds a-priori values 
    and has not been modified by the :ref:`optimizer`.

    StateVector objects can be filled with data in two ways
        1. By reading the data from file
        2. By creating the data through a set of method calls

    Option (1) is invoked using method :meth:`~da.baseclasses.statevector.StateVector.read_from_file`. 
    Option (2) consists of a call to method :meth:`~da.baseclasses.statevector.StateVector.make_new_ensemble`

    Once the StateVector object has been filled with data, it is used in the pipeline and a few more methods are
    invoked from there:
        * :meth:`~da.baseclasses.statevector.StateVector.propagate`, to advance the StateVector from t=t to t=t+1
        * :meth:`~da.baseclasses.statevector.StateVector.write_to_file`, to write the StateVector to a NetCDF file for later use

    The methods are described below:

    .. automethod:: da.baseclasses.statevector.StateVector.setup 
    .. automethod:: da.baseclasses.statevector.StateVector.read_from_file
    .. automethod:: da.baseclasses.statevector.StateVector.write_to_file
    .. automethod:: da.baseclasses.statevector.StateVector.make_new_ensemble
    .. automethod:: da.baseclasses.statevector.StateVector.propagate
    .. automethod:: da.baseclasses.statevector.StateVector.write_members_to_file

    Finally, the StateVector can be mapped to a gridded array, or to a vector of TransCom regions, using:

    .. automethod:: da.baseclasses.statevector.StateVector.grid2vector
    .. automethod:: da.baseclasses.statevector.StateVector.vector2grid
    .. automethod:: da.baseclasses.statevector.StateVector.vector2tc
    .. automethod:: da.baseclasses.statevector.StateVector.state2tc

    """

    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 setup(self, dacycle):
        """
        setup the object by specifying the dimensions. 
        There are two major requirements for each statvector that you want to build:
        
            (1) is that the statevector can map itself onto a regular grid
            (2) is that the statevector can map itself (mean+covariance) onto TransCom regions

        An example is given below.
        """

        self.nlag = int(dacycle['time.nlag'])
        self.nmembers = int(dacycle['da.optimizer.nmembers'])
#        self.nparams = int(dacycle.dasystem['nparameters'])
brunner's avatar
brunner committed
153
        self.nparams = 23
brunner's avatar
brunner committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
        self.nobs = 0
        
        self.obs_to_assimilate = ()  # empty containter to hold observations to assimilate later on

        # 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 EnsembleMember objects, we define member 0 as the mean of the distribution and n=1,...,nmembers as the spread.

        self.ensemble_members = list(range(self.nlag))

        for n in range(self.nlag):
            self.ensemble_members[n] = []


        # This specifies the file to read with the gridded mask at 1x1 degrees. Each gridbox holds a number that specifies the parametermember
        #  that maps onto it. From this map, a dictionary is created that allows a reverse look-up so that we can map parameters to a grid.

        mapfile = os.path.join(dacycle.dasystem['regionsfile'])
#        ncf = io.ct_read(mapfile, 'read')
#        self.gridmap = ncf.get_variable('regions')
#        self.tcmap = ncf.get_variable('transcom_regions')
#        ncf.close()
        regfile = Dataset(mapfile,mode='r')
        self.gridmap = np.squeeze(regfile.variables['regions'])
        self.tcmap = np.squeeze(regfile.variables['transcom_regions'])
        regfile.close()

        logging.debug("A TransCom  map on 1x1 degree was read from file %s" % dacycle.dasystem['regionsfile'])
        logging.debug("A parameter map on 1x1 degree was read from file %s" % dacycle.dasystem['regionsfile'])

        # Create a dictionary for state <-> gridded map conversions

brunner's avatar
brunner committed
185
        nparams = 23
brunner's avatar
brunner committed
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
        #nparams = self.gridmap.max()
        self.griddict = {}
        for r in range(1, 11):
#        for r in range(1, int(nparams) + 1):
            sel = np.nonzero(self.gridmap.flat == r)
            if len(sel[0]) > 0: 
                self.griddict[r] = sel
                self.griddict[r+11] = sel   # pavle - expand dictionary for nparam values because of RESP
        # pavle: sel sorts out regions by PFT

        logging.debug("A dictionary to map grids to states and vice versa was created")

        # Create a matrix for state <-> TransCom conversions

        self.tcmatrix = np.zeros((self.nparams, 23), 'float') 

        for r in range(1, self.nparams + 1):
            sel = np.nonzero(self.gridmap.flat == r)
            if len(sel[0]) < 1: 
                continue
            else:
                n_tc = set(self.tcmap.flatten().take(sel[0]))
                if len(n_tc) > 1: 
                    logging.error("Parameter %d seems to map to multiple TransCom regions (%s), I do not know how to handle this" % (r, n_tc))
                    raise ValueError
                self.tcmatrix[r - 1, int(n_tc.pop()) - 1] = 1.0

        logging.debug("A matrix to map states to TransCom regions and vice versa was created")

        # Create a mask for species/unknowns

        self.make_species_mask()

    def make_species_mask(self):

        """

        This method creates a dictionary with as key the name of a tracer, and as values an array of 0.0/1.0 values 
        specifying which StateVector elements are constrained by this tracer. This mask can be used in 
        the optimization to ensure that certain types of osbervations only update certain unknowns.

        An example would be that the tracer '14CO2' can be allowed to only map onto fossil fuel emissions in the state

        The form of the mask is:

        {'co2': np.ones(self.nparams), 'co2c14', np.zeros(self.nparams)  }

        so that 'co2' maps onto all parameters, and 'co2c14' on none at all. These arrays are used in the Class 
        optimizer when state updates are actually performed

        """
        self.speciesdict = {'co2': np.ones(self.nparams)}
        logging.debug("A species mask was created, only the following species are recognized in this system:")
        for k in self.speciesdict.keys(): 
            logging.debug("   ->    %s" % k)


    def make_new_ensemble(self, lag, 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 covariancematrix == None: 
 #           covariancematrix = np.identity(self.nparams)

        # Make a cholesky decomposition of the covariance matrix


        try:
            _, s, _ = np.linalg.svd(covariancematrix)
        except:
            s = np.linalg.svd(covariancematrix, full_matrices=1, compute_uv=0) #Cartesius fix
        dof = np.sum(s) ** 2 / sum(s ** 2)
brunner's avatar
brunner committed
270
#        covariancematrix = np.identity(self.nparams)   # pavle - for testing, final matrix will be 23x23
brunner's avatar
brunner committed
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
        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 mean values 

        newmean = np.ones(self.nparams, float) # standard value for a new time step is 1.0

        # If this is not the start of the filter, average previous two optimized steps into the mix
        if lag == self.nlag - 1 and self.nlag >= 3:
            newmean += self.ensemble_members[lag - 1][0].param_values + \
                                           self.ensemble_members[lag - 2][0].param_values 
            newmean = newmean / 3.0

        # 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_members[lag].append(newmember)

        # Create members 1:nmembers and add to ensemble_members list
        for member in range(1, self.nmembers):
brunner's avatar
brunner committed
295
296
297
            rands = np.random.uniform(low=-1., high=1., size=self.nparams-1)
            rands_bg = np.random.uniform(low=-0.05, high=0.05, size=1)

brunner's avatar
brunner committed
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
#            if member == 1 and lag == 0:
 #               rands=np.array([-0.2236048,0.1492591,-0.4134172,0.2030501,-0.3056341,0.6603037,0.3247675,0.6550926,0.8024981,0.1593491,-1,0.08171216,-0.4006799,-0.3123662, 0.1022556, 0.1977177, -0.4678199, 0.8751256, 0.9531838,0.9820679,0.3594339,-1,0.03129737])
  #          if member == 2 and lag == 0:
   #             rands=np.array([-0.1477133,-0.6238207,-0.8987545,-1,-1,0.1647741,0.6627979,0.4307449,-0.4548403,-0.1154654,-1,-0.05427124,-0.9482519,0.7116647,0.03384373,-0.09207293,0.4631868,-0.5646629,0.7750968,-0.8704156,0.1906742,-1,0.004112379])
    #        if member == 3 and lag == 0:
     #           rands=np.array([-0.5120361,-0.6635633,0.4805151,0.3135487,0.02116868,0.0115789,-0.2708648, 0.6396611,-0.4579386,-0.8980371,-1,0.4571423,0.2513106,0.3329439,0.5895487,-0.2265133,-0.3762723,0.9784772,0.6362864,-0.0876907,-0.1079749,-1,0.00497188])
      #      if member == 4 and lag == 0:
       #         rands=np.array([0.5773118,0.8732759,1.126802,-0.3782958,0.662484,1.199195,-0.5346411,-0.0550764,-0.6245533,0.6759356,-1,-0.8961901,0.3352093,-0.5532659,-0.7597447,0.5648068,-0.3713143,0.124998,-0.6321908,0.2655768,-0.02420997,-1,-0.01510416])
        #    if member == 5 and lag == 0:
         #       rands=np.array([0.1289302,-0.08508042,-0.3619103,-0.5414423,-0.07692706,-0.9641049,-0.268249,-0.7405945,-0.2146955,-0.1976647,-1,-0.1085423,-0.946457,-0.008138934,-0.530806,0.315332,-0.2215554,0.7034806,0.527541,-0.4504287,-0.3679401,-1,0.02846291])
          #  if member == 6 and lag == 0:
           #     rands=np.array([-0.2574656,-0.3259587,0.8530074,0.04455876,0.5508205,0.553733,0.5282572,-0.7782328,0.7970823,-0.4580753,-1,-0.1998393,0.04793931,0.01853116,0.6679623,-0.4434216,-0.1505937,-0.1531717,-0.1313033,0.3774659,-0.8304669,-1,0.0197916,])
            #if member == 7 and lag == 0:
             #   rands=np.array([0.6881989,0.3826058,-0.5495766,0.1132997,-0.57166,0.3214818,-0.9898017,0.6837509,-0.3501504,0.7928728,-1,-0.2640067,0.5886974,0.262185,0.2781147,0.7389106,0.6068769,0.1230468,-0.4471386,0.09954312,-0.4604041,-1,-0.04947798])
#            if member == 8 and lag == 0:
 #               rands=np.array([-0.8718347,-0.03896444,-0.384823,-0.5786896,0.2542168,-0.5035606,0.7728764,-0.6899376,-0.5994669,0.1059537,-1,-0.859535,0.2382655,-0.3679568,0.450976,0.6301627,-0.6471004,-0.8724778,0.7950776,-0.1000319,-0.1237619,-1,-0.04825296])
  #          if member == 9 and lag == 0:
   #             rands=np.array([0.06249616,0.78708,-0.5142601,-0.1251998,-0.5136412,-0.9106983,0.7841995,-0.8327041,-0.03150041,-0.7065354,-1,-0.1541739,0.7130719,0.6997306,0.5323737,-0.6959949,0.1497954,0.2615451,-0.451424,-0.8742322,-0.4250237,-1,0.03196495])
    #        if member == 10 and lag == 0:
     #           rands=np.array([-0.1827562,0.3512708,0.3367097,-0.6161442,0.2948794,0.3662989,0.7885535,-0.6752307,0.4505339,0.9947404,-1,-0.1427709,-0.4374289,-0.4191162,-0.9822996,-0.2934023,-0.2808521,-0.8965437, 0.7602426,-0.1778743,-0.000989437,-1,-0.01274748])
      #      if member == 11 and lag == 0:
       #         rands=np.array([0.3784615,-0.7083653,-0.821948,0.06502473,-0.6035533,-0.7092476,0.9460906,0.2481164,0.772471,0.6445355,-1,-0.1532128,-0.5447619,-0.5863768,0.2508805,-0.2562239,-0.9050175,-0.6196863,-1,-0.3799991,0.5817245,-1,-0.02425112])
        #    if member == 12 and lag == 0:
         #       rands=np.array([-0.3772987,-0.09633347,-0.9725795,-0.9913203,-0.1097328,0.3683293,-0.3654553,0.3282615,0.5734389,-0.03295924,-1,0.6370476,-0.4982804,0.174771,-0.740661,-0.02542321,-0.7887223,0.7312356,-0.3222232,-0.1027924,0.02633788,-1,-0.040399])
          #  if member == 13 and lag == 0:
           #     rands=np.array([0.7854867,0.3317518,0.6469278,0.5409229,-0.7125678,-0.5082978,0.4885957,-0.490857,-0.1638525,-0.3768145,-1,-0.8093954,0.5902002,0.03597473,-0.7890476,0.5711794,-0.3381689,-0.67668,-1,0.7070413,-0.1368402,-1,0.01318651])
            #if member == 14 and lag == 0:
             #   rands=np.array([0.7588015,0.4496917,-0.1509738,0.5066012,0.7388234,0.310389,-0.2899449,-0.6862295,-0.5775338,0.02781044,-1,-0.3429568,-0.7947903,-0.6043629,-0.8869596,-0.5914347,-0.2657907,-0.6107077,-0.2889067,-0.283401,0.8258069,-1,-0.03961588])
#            if member == 15 and lag == 0:
 #               rands=np.array([-0.1068678,-0.6150778,0.2749358,-0.006487241,0.3801498,-0.2595585,-0.1739224,0.4673905,0.2893905,-0.7730854,-1,-0.6906102,0.3519762,0.6020219,-0.1512664,0.4392967,0.4572283,-0.109156,0.1201844,-0.8120148,-1,-1,0.008621466])
  #          if member == 16 and lag == 0:
   #             rands=np.array([-0.4203033,0.6120592,-0.07144696,0.07225665,-0.7220178,-0.5390107,-0.3260586,0.7489389,0.4817865,0.04893194,-1,0.3428221,0.3335567,1.011872,0.3968997,1.008683,0.6700837,0.886711,0.5548246,0.4422186,-0.3110386,-1,0.04391268])
    #        if member == 17 and lag == 0:
     #           rands=np.array([0.6939667,-0.02326318,-0.6011316,0.5709327,-0.6587268,-0.2600601,0.09742739,-0.6605019,0.8611858,-0.3578327,-1,0.8069043,1.003004,0.5246789,-0.2310518,-0.3621024,-0.1325169,-0.2894324,0.5083101,-0.04920054,0.2199543,-1,0.009945367])
      #      if member == 18 and lag == 0:
       #         rands=np.array([0.3914793,-0.06666035,0.7248741,-0.4003048,-0.6606267,-0.807146,-0.5127164,-0.3519507,0.3180182,0.1996123,-1,-0.9029648,0.1214889,0.02123361,0.6533288,-0.3655816,0.7919229,-0.02215478,0.246939,0.4821936,-0.5654917,-1,-0.003499507])
        #    if member == 19 and lag == 0:
         #       rands=np.array([-0.8513823,-1,-0.9796664,0.3713398,-0.4331474,0.3044456,-0.761591,-0.5368128,0.2107361,0.3060099,-1,-0.2859453,-0.5175856,-0.4946607,0.3281591,-0.2430111,0.05990436,0.01103085,-0.6508834,-0.002251296,0.3083046,-1,0.005067795])
          #  if member == 20 and lag == 0:
           #     rands=np.array([0.1491016,-0.5305333,-0.4203315,-0.08089796,0.4151005,-0.7401157,0.3687969,-0.5454835,0.714487,1.023472,-1,0.3852504,0.2557997,0.1142526,-0.3280363,0.6893721,-0.4994218,0.7723718,1.014927,0.1936809,0.4526471,-1,0.01181981])
            #if member == 1 and lag == 1:
             #   rands=np.array([-0.8534008,-1,-0.01288013,-0.9141065,0.4783214,-0.1723187,-0.8843083,-0.1162029,-0.009133766,-0.6961404,-1,-0.992416,-0.2270698,-0.795909,-0.9905862,-1,-1,0.5997077,-0.08789548,0.9465016,-0.3178105,-1,0.02080303])
#            if member == 2 and lag == 1:
 #               rands=np.array([0.2119585,-0.06496529,-0.113645,-0.3219482,0.3354161,-0.2561546,-0.8578026,0.488537,0.4687303,0.9537283,-1,0.2903501,-0.2037528,0.2844993,-0.2640329,0.6626308,1.016451,-0.2284834,-0.9198995,-0.772432,-0.9456592,-1,-0.00193539])
  #          if member == 3 and lag == 1:
   #             rands=np.array([-0.8951237,-0.03746896,0.6338864,-0.7385425,-0.5215907,-0.2300481,-1,-0.5272307,-0.3447395,-0.161121,-1,-0.3376194,-0.5966941,0.5816291,0.4801027,0.8620167,0.02183707,-0.5763896,0.3380404,-0.4263453,0.5804913,-1,0.04326432])
    #        if member == 4 and lag == 1:
     #           rands=np.array([-0.3111334,0.1525611,-0.5140871,-0.9283414,0.08702877,-0.8027726,0.7558329,-0.3327861,-0.3074336,0.4950007,-1,-0.3270844,0.6640525,-0.2391291,-0.08126307,-0.01146206,0.01725924,0.7988746,1.104426,0.1875366,0.5716885,-1,-0.01126271])
      #      if member == 5 and lag == 1:
       #         rands=np.array([0.7825249,-0.4134479,-0.4143417,-0.9416251,0.497882,-0.7672758,0.4174937,0.8023338,0.9017133,0.4341536,-1,0.1975201,0.7627142,0.2342599,-0.5487284,-0.8285513,-0.1792318,0.7551883,0.7257428,0.6304003,0.4330811,-1,-0.03357022])
        #    if member == 6 and lag == 1:
         #       rands=np.array([-0.3209096,-0.5156574,0.4877877,-0.235484,-0.4256767,0.6646661,0.2886034,0.9879024,-0.06839389,0.6481704,-1,-0.5804082,-0.6514803,0.7308455,0.4101343,-0.8554199,-1,-0.254794,-0.1424363,-0.1296123,0.1170656,-1,0.007876628])
          #  if member == 7 and lag == 1:
           #     rands=np.array([-0.4625655,0.6606569,-0.1525902,-0.627054,-0.5481565,-0.7688485,-0.6856943,-0.3053931,-0.03701457,-0.3725091,-1,0.155129,0.136141,-0.7843549,-0.0782493,-0.5983688,0.4498052,0.7143801,-0.1611456,1.007675,0.9826942,-1,0.04915765])
            #if member == 8 and lag == 1:
             #   rands=np.array([-0.3430965,-0.4227155,-0.07735894,0.1337408,0.3767924,-0.3956357,-0.2084762,0.2861364,0.9596422,0.718588,-1,0.7803586,-0.2091809,0.02144417,-0.269459,0.5019733,0.4792759,0.4641314,1.016747,1.19465,0.3354242,-1,0.02544297])
#            if member == 9 and lag == 1:
 #               rands=np.array([-0.7412553,-0.5466392,0.6302599,0.6477938,0.3541165,-0.3394596,0.346593,0.7322681,-0.7179851,-0.731463,-1,0.5184832,-0.1218895,0.7145109,1.139154,0.1859014,-0.6142078,-0.6477512,0.799903,-0.2647903,0.264973,-1,0.01599488])
  #          if member == 10 and lag == 1:
   #             rands=np.array([-0.03720611,-0.4162885,0.3607477,0.223355,0.1092008,-0.3264467,-0.721734,-0.6520889,-0.6977112,-0.9723021,-1,-0.693766,-0.4895594,-0.6218721,-0.9569306,-0.009422544,-0.09222607,0.6750132,-0.7697481,-0.6925188,0.6922073,-1,0.04024281])
    #        if member == 11 and lag == 1:
     #           rands=np.array([-0.5116026,0.343556,-0.9518977,0.4320086,0.7536819,-0.6906458,-0.8697743,-0.5209014,-0.0006806303,-0.7983823,-1,-0.6035937,-0.6045027,0.2035089,-0.01446889,-0.4622703,0.04925987,-0.1290628,-0.8433856,-0.5047902,-0.4363325,-1,0.04609365])
      #      if member == 12 and lag == 1:
       #         rands=np.array([-0.03203867,0.1992727,0.01177486,-0.8682783,-0.6371214,-0.8578926,-0.5281639,0.4448915,0.05798108,0.4604711,-1,-0.6847421,-1,-0.2765961,-0.4123928,-0.1756036,0.459426,0.303864,-0.6993556,0.1463735,0.06817089,-1,-0.02779672])
        #    if member == 13 and lag == 1:
         #       rands=np.array([-0.1434858,0.6530933,0.4572614,1.004458,-0.1688031,0.2899839,-0.4890676,-0.08017457,0.2399375,0.8467381,-1,-0.9818414,-0.7904143,0.6110988,-0.2375934,0.008658743,-0.6596891,0.804915,0.8108714,1.026871,-0.632525,-1,0.008359899])
          #  if member == 14 and lag == 1:
           #     rands=np.array([0.1951734,0.1356803,-0.7984297,0.6222939,-0.3783839,-0.2109271,-0.5047454,0.7663568,-0.2467663,0.7969002,-1,-0.1039866,-0.0799861,-0.9506351,-0.2622426,-0.5555511,-0.6443623,-0.3867293,0.8174118,-0.2623511,-0.6504129,-1,0.02077365])
            #if member == 15 and lag == 1:
             #   rands=np.array([0.5399131,0.9709167,0.1353998,0.9513429,0.1061087,-0.1906423,0.2546594,0.6244627,0.4253108,-0.5599692,-1,-0.6188933,-0.009267714,0.5136899,0.7686136,-0.7976044,0.6507949,0.2095091,-0.8000021,0.7896537,-0.6757631,-1,-0.04675981])
#            if member == 16 and lag == 1:
 #               rands=np.array([0.2039032,0.4437117,-0.8730499,-0.8436514,-0.3197937,0.1300565,0.2079234,0.07507716,-0.1368318,-0.8760222,-1,0.7775679,0.1503317,-0.6096121,-0.4738224,-0.4769458,0.01603237,-0.5856133,0.1603515,0.639724,0.8044817,-1,-0.03076992])
  #          if member == 17 and lag == 1:
   #             rands=np.array([0.6912341,0.03829452,0.0645803,-0.1688701,-0.06769498,0.5935946,0.794823,-0.1764786,-0.5644639,0.1006729,-1,0.3818233,0.8611205,0.8617951,1.227797,0.06026031,-0.3915816,-0.2597351,0.03994489,0.04319654,0.4627483,-1,0.02341671])
    #        if member == 18 and lag == 1:
     #           rands=np.array([0.1631074,-0.7144772,0.5909197,0.8623705,-0.157497,0.5982791,-0.2239446,0.5731422,0.6889938,-0.04528522,-1,0.07149924,-0.6744667,0.4025916,0.9270787,0.1585069,0.8808218,0.5500813,1.050181,0.6608209,-0.4115532,-1,0.0001819044])
      #      if member == 19 and lag == 1:
       #         rands=np.array([-0.4331279,-0.2030248,-0.2923387,-0.1422065,0.04285587,-0.6101277,-0.5285338,0.820045,0.4597438,-0.7786282,-1,0.3376869,0.1703605,-0.5970429,-0.2520662,-0.6760406,0.3587383,0.1665516,-0.351487,-0.3638605,0.4790523,-1,-0.01560173])
        #    if member == 20 and lag == 1:
         #       rands=np.array([0.4845566,0.6727741,0.9107583,-0.2242058,0.8334284,-0.0900023,-0.1698219,-0.1709382,-0.1608661,-0.9621966,-1,0.9589482,-0.5061898,0.02963674,-0.5164318,-0.05424574,-0.8521295,0.7594227,-0.5421007,0.8706058,0.7356571,-1,-0.03166761])
            #newmember.param_values = rands + newmean
brunner's avatar
brunner committed
379
            newmember = EnsembleMember(member)
brunner's avatar
brunner committed
380
            newmember.param_values = (np.hstack((np.dot(C, rands[0:11]),np.dot(C, rands[11:22]), rands_bg)) + newmean).ravel()
brunner's avatar
brunner committed
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
            newmember.param_values[10] = 0.
            newmember.param_values[21] = 0.
            newmember.param_values[newmember.param_values<0.] = 0.
            self.ensemble_members[lag].append(newmember)

        logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, (lag + 1)))


    def propagate(self, dacycle):
        """ 
        :rtype: None

        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.

        """
        
        # Remove State Vector n=1 by simply "popping" it from the list and appending a new empty list at the front. This empty list will
        # hold the new ensemble for the new cycle 

        self.ensemble_members.pop(0)
        self.ensemble_members.append([])

        # And now create a new time step of mean + members for n=nlag
        date = dacycle['time.start'] + timedelta(days=(self.nlag - 0.5) * int(dacycle['time.cycle']))
        cov = self.get_covariance(date, dacycle)
        self.make_new_ensemble(self.nlag - 1, cov)

        logging.info('The state vector has been propagated by one cycle')


    def write_to_file(self, filename, qual):
        """
        :param filename: the full filename for the output NetCDF file
        :rtype: None

        Write the StateVector information to a NetCDF file for later use. 
        In principle the output file will have only one two datasets inside 
        called:
            * `meanstate`, dimensions [nlag, nparamaters]
            * `ensemblestate`, dimensions [nlag,nmembers, nparameters]

        This NetCDF information can be read back into a StateVector object using 
        :meth:`~da.baseclasses.statevector.StateVector.read_from_file`

        """
        #import da.tools.io4 as io
        #import da.tools.io as io

        if qual == 'prior':
            f = io.CT_CDF(filename, method='create')
            logging.debug('Creating new StateVector output file (%s)' % filename)
            #qual = 'prior'
        else:
            f = io.CT_CDF(filename, method='write')
            logging.debug('Opening existing StateVector output file (%s)' % filename)
            #qual = 'opt'

        dimparams = f.add_params_dim(self.nparams)
        dimmembers = f.add_members_dim(self.nmembers)
        dimlag = f.add_lag_dim(self.nlag, unlimited=True)

        for n in range(self.nlag):
            members = self.ensemble_members[n]
            mean_state = members[0].param_values

            savedict = f.standard_var(varname='meanstate_%s' % qual)
            savedict['dims'] = dimlag + dimparams 
            savedict['values'] = mean_state
            savedict['count'] = n
            savedict['comment'] = 'this represents the mean of the ensemble'
            f.add_data(savedict)

            members = self.ensemble_members[n]
            devs = np.asarray([m.param_values.flatten() for m in members])
            data = devs - np.asarray(mean_state)

            savedict = f.standard_var(varname='ensemblestate_%s' % qual)
            savedict['dims'] = dimlag + dimmembers + dimparams 
            savedict['values'] = data
            savedict['count'] = n
            savedict['comment'] = 'this represents deviations from the mean of the ensemble'
            f.add_data(savedict)
        f.close()

        logging.info('Successfully wrote the State Vector to file (%s) ' % filename)

    def read_from_file(self, filename, qual='opt'):
        """ 
        :param filename: the full filename for the input NetCDF file
        :param qual: a string indicating whether to read the 'prior' or 'opt'(imized) StateVector from file
        :rtype: None

        Read the StateVector information from a NetCDF file and put in a StateVector object
        In principle the input file will have only one four datasets inside 
        called:
            * `meanstate_prior`, dimensions [nlag, nparamaters]
            * `ensemblestate_prior`, dimensions [nlag,nmembers, nparameters]
            * `meanstate_opt`, dimensions [nlag, nparamaters]
            * `ensemblestate_opt`, dimensions [nlag,nmembers, nparameters]

        This NetCDF information can be written to file using 
        :meth:`~da.baseclasses.statevector.StateVector.write_to_file`

        """

        #import da.tools.io as io
        f = io.ct_read(filename, 'read')
        meanstate = f.get_variable('statevectormean_' + qual)
        ensmembers = f.get_variable('statevectorensemble_' + qual)
        f.close()

        for n in range(self.nlag):
            if not self.ensemble_members[n] == []:
                self.ensemble_members[n] = []
                logging.warning('Existing ensemble for lag=%d was removed to make place for newly read data' % (n + 1))

            for m in range(self.nmembers):
                newmember = EnsembleMember(m)
                newmember.param_values = ensmembers[n, m, :].flatten() + meanstate[n]  # add the mean to the deviations to hold the full parameter values
                self.ensemble_members[n].append(newmember)

        logging.info('Successfully read the State Vector from file (%s) ' % filename)

    def write_members_to_file(self, lag, outdir,endswith='.nc'):
        """ 
           :param: lag: Which lag step of the filter to write, must lie in range [1,...,nlag]
           :param: outdir: Directory where to write files
           :param: endswith: Optional label to add to the filename, default is simply .nc
           :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%s' % (mem.membernumber, endswith))
            filename = os.path.join(outdir, 'parameters_lag'+str(lag)+'.%03d%s' % (mem.membernumber, endswith))
            ncf = io.CT_CDF(filename, method='create')
            dimparams = ncf.add_params_dim(self.nparams)
            dimgrid = ncf.add_latlon_dim()

            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' % 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' % mem.membernumber
            ncf.add_data(savedict)

            ncf.close()

            logging.debug('Successfully wrote data from ensemble member %d to file (%s) ' % (mem.membernumber, filename))

    def write_members_for_cosmo(self, lag, outdir,endswith='.nc'):
        members = self.ensemble_members[lag]

        self.nparams=int(self.nparams)
        for mem in members:

# GPP
            filename_gpp = os.path.join(outdir, 'parameters_gpp_lag'+str(lag)+'.%03d%s' % (mem.membernumber, endswith))
            ncf = io.CT_CDF(filename_gpp, method='create')
brunner's avatar
brunner committed
576
577
            dimparams = ncf.add_params_dim(11)
            #dimparams = ncf.add_params_dim((self.nparams-1)//2)
brunner's avatar
brunner committed
578
579
            dimgrid = ncf.add_latlon_dim()

brunner's avatar
brunner committed
580
            data = mem.param_values[0:11]
brunner's avatar
brunner committed
581
            #data = 2*np.ones(11)-mem.param_values[0:11]       # unblock this maybe
brunner's avatar
brunner committed
582
            #data = mem.param_values[0:(self.nparams-1)//2]
brunner's avatar
brunner committed
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607

            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' % 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' % mem.membernumber
            ncf.add_data(savedict)

            ncf.close()
# RESP
            filename_resp = os.path.join(outdir, 'parameters_resp_lag'+str(lag)+'.%03d%s' % (mem.membernumber, endswith))
            ncf = io.CT_CDF(filename_resp, method='create')
brunner's avatar
brunner committed
608
609
            dimparams = ncf.add_params_dim(11)
            #dimparams = ncf.add_params_dim((self.nparams-1)//2)
brunner's avatar
brunner committed
610
611
            dimgrid = ncf.add_latlon_dim()

brunner's avatar
brunner committed
612
613
            data = mem.param_values[11:22]
            #data = mem.param_values[(self.nparams-1)//2:self.nparams-1]
brunner's avatar
brunner committed
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636

            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' % 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' % mem.membernumber
            ncf.add_data(savedict)

            ncf.close()

brunner's avatar
brunner committed
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
# BACKGROUND CO2
#            filename_bg = os.path.join(outdir, 'parameters_bg_lag'+str(lag)+'.%03d%s' % (mem.membernumber, endswith))
 #           ncf = io.CT_CDF(filename_bg, method='create')
  #          dimparams = ncf.add_params_dim(11)
   #         dimgrid = ncf.add_latlon_dim()
#
 #           data = mem.param_values[-1]
  #          data = np.hstack((data, data, data, data, data, data, data, data, data, data, data)).ravel()
#
 #           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' % 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' % mem.membernumber
        #    ncf.add_data(savedict)
#
 #           ncf.close()
#
#
brunner's avatar
brunner committed
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
    def grid2vector(self, griddata=None, method='avg'):
    # not used --pavle
        """ 
            Map gridded data onto a vector of length (nparams,)

           :param griddata: a gridded dataset to use. This dataset is mapped onto a vector of length `nparams`
           :param method: a string that specifies the method to combine grid boxes in case reverse=True. Must be either ['avg','sum','minval']
           :rtype: ndarray: size (nparameters,)

           This method makes use of a dictionary that links every parameter number [1,...,nparams] to a series of gridindices. These 
           indices specify a location on a 360x180 array, stretched into a vector using `array.flat`. There are multiple ways of calling 
           this method::

               values       = self.grid2vector(griddata=mygriddeddata,method='minval') # 
                                                                    using the minimum value of all datapoints covered by that parameter index

               values       = self.grid2vector(griddata=mygriddeddata,method='avg') # 
                                                                    using the average value of all datapoints covered by that parameter index

               values       = self.grid2vector(griddata=mygriddeddata,method='sum') # 
                                                                    using the sum of values of all datapoints covered by that parameter index

           .. note:: This method uses a DaSystem object that must be initialized with a proper parameter map. See :class:`~da.baseclasses.dasystem` for details

        """

        methods = ['avg', 'sum', 'minval']
        if method not in methods:
            logging.error("To put data from a map into the statevector, please specify the method to use (%s)" % methods)
            raise ValueError

        result = np.zeros((self.nparams,), float)
        for k, v in self.griddict.items():
#            print(k,k-1,result.shape, v)
            if method == "avg": 
                result[k - 1] = griddata.take(v).mean()
            elif method == "sum" : 
                result[k - 1] = griddata.take(v).sum()
            elif method == "minval" : 
                result[k - 1] = griddata.take(v).min()
        return result # Note that the result is returned, but not yet placed in the member.param_values attrtibute!


    def vector2grid(self, vectordata=None):
        """ 
            Map vector elements to a map or vice cersa

           :param vectordata: a vector dataset to use in case `reverse = False`. This dataset is mapped onto a 1x1 grid and must be of length `nparams`
           :rtype: ndarray: an array of size (360,180,) 

           This method makes use of a dictionary that links every parameter number [1,...,nparams] to a series of gridindices. These 
           indices specify a location on a 360x180 array, stretched into a vector using `array.flat`. There are multiple ways of calling 
           this method::

               griddedarray = self.vector2grid(vectordata=param_values) # simply puts the param_values onto a (180,360,) array

           .. note:: This method uses a DaSystem object that must be initialzied with a proper parameter map. See :class:`~da.baseclasses.dasystem` for details

        """
        result = np.zeros(self.gridmap.shape, float)
        for k, v in self.griddict.items():
#            print(k,v)
            if k<=11:
                result.put(v, vectordata[k - 1])
        return result         

    def vector2tc(self, vectordata, cov=False):
        """ 
            project Vector onto TransCom regions 

           :param vectordata: a vector dataset to use, must be of length `nparams`
           :param cov: a Boolean to specify whether the input dataset is a vector (mean), or a matrix (covariance)
           :rtype: ndarray: an array of size (23,) (cov:F) or of size (23,23,) (cov:T)
        """

        M = self.tcmatrix
        if cov:
            return np.dot(np.transpose(M), np.dot(vectordata, M))
        else:
            return np.dot(vectordata.squeeze(), M)

    def state_to_grid(self, fluxvector=None, lag=1):
        """ 
            Transforms the StateVector information (mean + covariance) to a 1x1 degree grid.

            :param: fluxvector: a vector of length (nparams,) that holds the fluxes associated with each parameter in the StateVector
            :param: lag: the lag at which to evaluate the StateVector
            :rtype: a tuple of two arrays (gridmean,gridvariance) with dimensions (180,360,)

            If the attribute `fluxvector` is not passed, the function will return the mean parameter value and its variance on a 1x1 map.
            
            ..note:: Although we can return the variance information for each gridbox, the covariance information contained in the original ensemble is lost when mapping to 1x1 degree!

        """

        if fluxvector == None:
            fluxvector = np.ones(self.nparams)

        ensemble = self.ensemble_members[lag - 1]
        ensemblemean = ensemble[0].param_values

        # First transform the mean
        gridmean = self.vector2grid(vectordata=ensemblemean * fluxvector)

        # And now the covariance, first create covariance matrix (!), and then multiply
        deviations = np.array([mem.param_values * fluxvector - ensemblemean for mem in ensemble])
        ensemble = []
        for mem in deviations:
            ensemble.append(self.vector2grid(mem))

        return (gridmean, np.array(ensemble))

    def state2tc(self, fluxvector=None, lag=1):
        """ 
            Transforms the StateVector information (mean + covariance) to the TransCom regions.

            :param: fluxvector: a vector of length (nparams,) that holds the fluxes associated with each parameter in the StateVector
            :param: lag: the lag at which to evaluate the StateVector
            :rtype: a tuple of two arrays (mean,covariance) with dimensions ((23,), (23,23,) )

        """
        ensemble = self.ensemble_members[lag - 1]
        ensemblemean = ensemble[0].param_values

        # First transform the mean

        mean = self.vector2tc(vectordata=ensemble[0].param_values * fluxvector)

        # And now the covariance, first create covariance matrix (!), and then multiply

        deviations = np.array([mem.param_values * fluxvector - ensemblemean for mem in ensemble])
        covariance = np.dot(np.transpose(deviations), deviations) / (self.nmembers - 1)
        cov = self.vector2tc(covariance, cov=True)

        return (mean, cov)

    def get_covariance(self, date, cycleparams):
        pass
    
################### End Class StateVector ###################

if __name__ == "__main__":
    pass