optimizer.py 19.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
"""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/>."""
13
14
15
16
#!/usr/bin/env python
# optimizer.py

"""
Peters, Wouter's avatar
Peters, Wouter committed
17
18
.. module:: optimizer
.. moduleauthor:: Wouter Peters 
19
20
21
22
23
24
25

Revision History:
File created on 28 Jul 2010.

"""

import logging
karolina's avatar
karolina committed
26
27
28
import numpy as np
import numpy.linalg as la
import da.tools.io4 as io
29
30

identifier = 'Optimizer baseclass'
karolina's avatar
karolina committed
31
version = '0.0'
32
33
34
35
36
37
38
39
40
41
42
43

################### Begin Class Optimizer ###################

class Optimizer(object):
    """
        This creates an instance of an optimization object. It handles the minimum least squares optimization
        of the state vector given a set of sample objects. Two routines will be implemented: one where the optimization
        is sequential and one where it is the equivalent matrix solution. The choice can be made based on considerations of speed
        and efficiency.
    """

    def __init__(self):
karolina's avatar
karolina committed
44
45
        self.ID = identifier
        self.version = version
46

karolina's avatar
karolina committed
47
        logging.info('Optimizer object initialized: %s' % self.ID)
48

49
    def setup(self, dims):
karolina's avatar
karolina committed
50
51
52
53
        self.nlag = dims[0]
        self.nmembers = dims[1]
        self.nparams = dims[2]
        self.nobs = dims[3]
54
        self.create_matrices()
55

56
    def create_matrices(self):
57
58
59
        """ Create Matrix space needed in optimization routine """

        # mean state  [X]
karolina's avatar
karolina committed
60
        self.x = np.zeros((self.nlag * self.nparams,), float)
61
        # deviations from mean state  [X']
karolina's avatar
karolina committed
62
        self.X_prime = np.zeros((self.nlag * self.nparams, self.nmembers,), float)
63
        # mean state, transported to observation space [ H(X) ]
karolina's avatar
karolina committed
64
        self.Hx = np.zeros((self.nobs,), float)
65
        # deviations from mean state, transported to observation space [ H(X') ]
karolina's avatar
karolina committed
66
        self.HX_prime = np.zeros((self.nobs, self.nmembers), float)
67
        # observations
karolina's avatar
karolina committed
68
        self.obs = np.zeros((self.nobs,), float)
69
        # observation ids
karolina's avatar
karolina committed
70
        self.obs_ids = np.zeros((self.nobs,), float)
71
        # covariance of observations
72
73
74
75
76
        # Total covariance of fluxes and obs in units of obs [H P H^t + R]
        if self.algorithm == 'Serial':
            self.R = np.zeros((self.nobs,), float)
            self.HPHR = np.zeros((self.nobs,), float)
        else:
karolina's avatar
karolina committed
77
            self.R = np.zeros((self.nobs, self.nobs,), float)
78
            self.HPHR = np.zeros((self.nobs, self.nobs,), float)
79
        # localization of obs
karolina's avatar
karolina committed
80
        self.may_localize = np.zeros(self.nobs, bool)
81
        # rejection of obs
karolina's avatar
karolina committed
82
        self.may_reject = np.zeros(self.nobs, bool)
83
        # flags of obs
karolina's avatar
karolina committed
84
        self.flags = np.zeros(self.nobs, int)
85
        # species type
karolina's avatar
karolina committed
86
        self.species = np.zeros(self.nobs, str)
87
        # species type
karolina's avatar
karolina committed
88
        self.sitecode = np.zeros(self.nobs, str)
89
90

        # species mask
karolina's avatar
karolina committed
91
        self.speciesmask = {}
92

93
        # Kalman Gain matrix
94
95
        #self.KG = np.zeros((self.nlag * self.nparams, self.nobs,), float)
        self.KG = np.zeros((self.nlag * self.nparams,), float)
karolina's avatar
karolina committed
96

97
    def state_to_matrix(self, statevector):
karolina's avatar
karolina committed
98
99
100
101
102
103
104
105
        allsites = []      # collect all obs for n=1,..,nlag
        allobs = []      # collect all obs for n=1,..,nlag
        allmdm = []      # collect all mdm for n=1,..,nlag
        allids = []  # collect all model samples for n=1,..,nlag
        allreject = []  # collect all model samples for n=1,..,nlag
        alllocalize = []  # collect all model samples for n=1,..,nlag
        allflags = []  # collect all model samples for n=1,..,nlag
        allspecies = []  # collect all model samples for n=1,..,nlag
106
        allsimulated = []  # collect all members model samples for n=1,..,nlag
107

108
        for n in range(self.nlag):
109
            samples = statevector.obs_to_assimilate[n]
110
111
112
            members = statevector.ensemble_members[n]
            self.x[n * self.nparams:(n + 1) * self.nparams] = members[0].param_values
            self.X_prime[n * self.nparams:(n + 1) * self.nparams, :] = np.transpose(np.array([m.param_values for m in members]))
113

karolina's avatar
karolina committed
114
115
            if samples != None:        
                self.rejection_threshold = samples.rejection_threshold
Peters, Wouter's avatar
Peters, Wouter committed
116

karolina's avatar
karolina committed
117
118
119
120
121
122
123
124
                allreject.extend(samples.getvalues('may_reject'))
                alllocalize.extend(samples.getvalues('may_localize'))
                allflags.extend(samples.getvalues('flag'))
                allspecies.extend(samples.getvalues('species'))
                allobs.extend(samples.getvalues('obs'))
                allsites.extend(samples.getvalues('code'))
                allmdm.extend(samples.getvalues('mdm'))
                allids.extend(samples.getvalues('id'))
125

karolina's avatar
karolina committed
126
                simulatedensemble = samples.getvalues('simulated')
brunner's avatar
brunner committed
127
128
                for s in range(simulatedensemble.shape[0]):
                        allsimulated.append(simulatedensemble[s])
129

karolina's avatar
karolina committed
130
131
132
133
        self.obs[:] = np.array(allobs)
        self.obs_ids[:] = np.array(allids)
        self.HX_prime[:, :] = np.array(allsimulated)
        self.Hx[:] = self.HX_prime[:, 0]
134

karolina's avatar
karolina committed
135
136
137
138
139
        self.may_reject[:] = np.array(allreject)
        self.may_localize[:] = np.array(alllocalize)
        self.flags[:] = np.array(allflags)
        self.species[:] = np.array(allspecies)
        self.sitecode = allsites
140

karolina's avatar
karolina committed
141
142
        self.X_prime = self.X_prime - self.x[:, np.newaxis] # make into a deviation matrix
        self.HX_prime = self.HX_prime - self.Hx[:, np.newaxis] # make a deviation matrix
143

144
145
146
147
148
        if self.algorithm == 'Serial':
            for i, mdm in enumerate(allmdm):
                self.R[i] = mdm ** 2
        else:
            for i, mdm in enumerate(allmdm):
karolina's avatar
karolina committed
149
                self.R[i, i] = mdm ** 2
150
                    
151
    def matrix_to_state(self, statevector):
152
        for n in range(self.nlag):
153
            members = statevector.ensemble_members[n]
karolina's avatar
karolina committed
154
            for m, mem in enumerate(members):
155
                members[m].param_values[:] = self.X_prime[n * self.nparams:(n + 1) * self.nparams, m] + self.x[n * self.nparams:(n + 1) * self.nparams]     
156

karolina's avatar
karolina committed
157
        logging.debug('Returning optimized data to the StateVector, setting "StateVector.isOptimized = True" ')
158

karolina's avatar
karolina committed
159
    def write_diagnostics(self, filename, type):
160
161
162
163
164
165
166
167
168
169
        """
            Open a NetCDF file and write diagnostic output from optimization process:

                - calculated residuals
                - model-data mismatches
                - HPH^T
                - prior ensemble of samples
                - posterior ensemble of samples
                - prior ensemble of fluxes
                - posterior ensemble of fluxes
170
171

            The type designation refers to the writing of prior or posterior data and is used in naming the variables"
172
173
        """

174
175
176
        # Open or create file

        if type == 'prior':
karolina's avatar
karolina committed
177
178
            f = io.CT_CDF(filename, method='create')
            logging.debug('Creating new diagnostics file for optimizer (%s)' % filename)
179
        elif type == 'optimized':
karolina's avatar
karolina committed
180
181
            f = io.CT_CDF(filename, method='write')
            logging.debug('Opening existing diagnostics file for optimizer (%s)' % filename)
182
183

        # Add dimensions 
184

185
186
187
188
189
190
        dimparams = f.add_params_dim(self.nparams)
        dimmembers = f.add_members_dim(self.nmembers)
        dimlag = f.add_lag_dim(self.nlag, unlimited=False)
        dimobs = f.add_obs_dim(self.nobs)
        dimstate = f.add_dim('nstate', self.nparams * self.nlag)
        dim200char = f.add_dim('string_of200chars', 200)
191

192
193
        # Add data, first the ones that are written both before and after the optimization

karolina's avatar
karolina committed
194
195
196
197
198
        savedict = io.std_savedict.copy() 
        savedict['name'] = "statevectormean_%s" % type
        savedict['long_name'] = "full_statevector_mean_%s" % type
        savedict['units'] = "unitless"
        savedict['dims'] = dimstate
karolina's avatar
karolina committed
199
        savedict['values'] = self.x.tolist()
karolina's avatar
karolina committed
200
        savedict['comment'] = 'Full %s state vector mean ' % type
201
        f.add_data(savedict)
karolina's avatar
karolina committed
202
203
204
205
206
207

        savedict = io.std_savedict.copy()
        savedict['name'] = "statevectordeviations_%s" % type
        savedict['long_name'] = "full_statevector_deviations_%s" % type
        savedict['units'] = "unitless"
        savedict['dims'] = dimstate + dimmembers
karolina's avatar
karolina committed
208
        savedict['values'] = self.X_prime.tolist()
karolina's avatar
karolina committed
209
        savedict['comment'] = 'Full state vector %s deviations as resulting from the optimizer' % type
210
        f.add_data(savedict)
karolina's avatar
karolina committed
211
212
213
214
215
216

        savedict = io.std_savedict.copy()
        savedict['name'] = "modelsamplesmean_%s" % type
        savedict['long_name'] = "modelsamplesforecastmean_%s" % type
        savedict['units'] = "mol mol-1"
        savedict['dims'] = dimobs
karolina's avatar
karolina committed
217
        savedict['values'] = self.Hx.tolist()
218
        savedict['comment'] = '%s mean mole fractions based on %s state vector' % (type, type)
219
        f.add_data(savedict)
karolina's avatar
karolina committed
220
221
222
223
224
225

        savedict = io.std_savedict.copy()
        savedict['name'] = "modelsamplesdeviations_%s" % type
        savedict['long_name'] = "modelsamplesforecastdeviations_%s" % type
        savedict['units'] = "mol mol-1"
        savedict['dims'] = dimobs + dimmembers
karolina's avatar
karolina committed
226
        savedict['values'] = self.HX_prime.tolist()
227
        savedict['comment'] = '%s mole fraction deviations based on %s state vector' % (type, type)
228
        f.add_data(savedict)
229

230
231
232
233
        # Continue with prior only data

        if type == 'prior':

karolina's avatar
karolina committed
234
235
236
237
238
            savedict = io.std_savedict.copy()
            savedict['name'] = "sitecode"
            savedict['long_name'] = "site code propagated from observation file"
            savedict['dtype'] = "char"
            savedict['dims'] = dimobs + dim200char
karolina's avatar
karolina committed
239
            savedict['values'] = self.sitecode
240
            savedict['missing_value'] = '!'
241
            f.add_data(savedict)
karolina's avatar
karolina committed
242
243
244
245
246
247

            savedict = io.std_savedict.copy()
            savedict['name'] = "observed"
            savedict['long_name'] = "observedvalues"
            savedict['units'] = "mol mol-1"
            savedict['dims'] = dimobs
karolina's avatar
karolina committed
248
            savedict['values'] = self.obs.tolist()
karolina's avatar
karolina committed
249
            savedict['comment'] = 'Observations used in optimization'
250
            f.add_data(savedict)
karolina's avatar
karolina committed
251
252
253
254
255
256
257

            savedict = io.std_savedict.copy()
            savedict['name'] = "obspack_num"
            savedict['dtype'] = "int64"
            savedict['long_name'] = "Unique_ObsPack_observation_number"
            savedict['units'] = ""
            savedict['dims'] = dimobs
karolina's avatar
karolina committed
258
            savedict['values'] = self.obs_ids.tolist()
karolina's avatar
karolina committed
259
            savedict['comment'] = 'Unique observation number across the entire ObsPack distribution'
260
            f.add_data(savedict)
karolina's avatar
karolina committed
261
262
263
264
265

            savedict = io.std_savedict.copy()
            savedict['name'] = "modeldatamismatchvariance"
            savedict['long_name'] = "modeldatamismatch variance"
            savedict['units'] = "[mol mol-1]^2"
266
267
268
            if self.algorithm == 'Serial':
                savedict['dims'] = dimobs
            else: savedict['dims'] = dimobs + dimobs
karolina's avatar
karolina committed
269
            savedict['values'] = self.R.tolist()
karolina's avatar
karolina committed
270
            savedict['comment'] = 'Variance of mole fractions resulting from model-data mismatch'
271
            f.add_data(savedict)
272
273
274
275

        # Continue with posterior only data

        elif type == 'optimized':
karolina's avatar
karolina committed
276
            
karolina's avatar
karolina committed
277
278
279
280
            savedict = io.std_savedict.copy()
            savedict['name'] = "totalmolefractionvariance"
            savedict['long_name'] = "totalmolefractionvariance"
            savedict['units'] = "[mol mol-1]^2"
281
282
283
            if self.algorithm == 'Serial':
                savedict['dims'] = dimobs
            else: savedict['dims'] = dimobs + dimobs
karolina's avatar
karolina committed
284
            savedict['values'] = self.HPHR.tolist()
karolina's avatar
karolina committed
285
            savedict['comment'] = 'Variance of mole fractions resulting from prior state and model-data mismatch'
286
            f.add_data(savedict)
karolina's avatar
karolina committed
287
288
289
290
291
292

            savedict = io.std_savedict.copy()
            savedict['name'] = "flag"
            savedict['long_name'] = "flag_for_obs_model"
            savedict['units'] = "None"
            savedict['dims'] = dimobs
karolina's avatar
karolina committed
293
            savedict['values'] = self.flags.tolist()
karolina's avatar
karolina committed
294
            savedict['comment'] = 'Flag (0/1/2/99) for observation value, 0 means okay, 1 means QC error, 2 means rejected, 99 means not sampled'
295
            f.add_data(savedict)
karolina's avatar
karolina committed
296

297
298
299
300
301
302
303
            #savedict = io.std_savedict.copy()
            #savedict['name'] = "kalmangainmatrix"
            #savedict['long_name'] = "kalmangainmatrix"
            #savedict['units'] = "unitless molefraction-1"
            #savedict['dims'] = dimstate + dimobs
            #savedict['values'] = self.KG.tolist()
            #savedict['comment'] = 'Kalman gain matrix of all obs and state vector elements'
304
            #dummy                   = f.add_data(savedict)
305

karolina's avatar
karolina committed
306
307
        f.close()
        logging.debug('Diagnostics file closed')
308

309

310
    def serial_minimum_least_squares(self):
311
312
313
        """ Make minimum least squares solution by looping over obs"""
        for n in range(self.nobs):

314
315
316
            # Screen for flagged observations (for instance site not found, or no sample written from model)

            if self.flags[n] != 0:
317
                logging.debug('Skipping observation (%s,%i) because of flag value %d' % (self.sitecode[n], self.obs_ids[n], self.flags[n]))
318
319
                continue

320
321
            # Screen for outliers greather than 3x model-data mismatch, only apply if obs may be rejected

karolina's avatar
karolina committed
322
            res = self.obs[n] - self.Hx[n]
323
324

            if self.may_reject[n]:
325
                threshold = self.rejection_threshold * np.sqrt(self.R[n])
326
                if np.abs(res) > threshold:
327
                    logging.debug('Rejecting observation (%s,%i) because residual (%f) exceeds threshold (%f)' % (self.sitecode[n], self.obs_ids[n], res, threshold))
328
                    self.flags[n] = 2
329
                    continue
330

331
332
            logging.debug('Proceeding to assimilate observation %s, %i' % (self.sitecode[n], self.obs_ids[n]))

karolina's avatar
karolina committed
333
            PHt = 1. / (self.nmembers - 1) * np.dot(self.X_prime, self.HX_prime[n, :])
334
            self.HPHR[n] = 1. / (self.nmembers - 1) * (self.HX_prime[n, :] * self.HX_prime[n, :]).sum() + self.R[n]
335

336
            self.KG[:] = PHt / self.HPHR[n]
337

338
            if self.may_localize[n]:
karolina's avatar
karolina committed
339
                logging.debug('Trying to localize observation %s, %i' % (self.sitecode[n], self.obs_ids[n]))
340
                self.localize(n)
341
            else:
karolina's avatar
karolina committed
342
                logging.debug('Not allowed to localize observation %s, %i' % (self.sitecode[n], self.obs_ids[n]))
343

344
            alpha = np.double(1.0) / (np.double(1.0) + np.sqrt((self.R[n]) / self.HPHR[n]))
345

346
            self.x[:] = self.x + self.KG[:] * res
347
348

            for r in range(self.nmembers):
349
                self.X_prime[:, r] = self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r])
350
351
352
353

#WP !!!! Very important to first do all obervations from n=1 through the end, and only then update 1,...,n. The current observation
#WP      should always be updated last because it features in the loop of the adjustments !!!!

karolina's avatar
karolina committed
354
355
            for m in range(n + 1, self.nobs):
                res = self.obs[n] - self.Hx[n]
356
                fac = 1.0 / (self.nmembers - 1) * (self.HX_prime[n, :] * self.HX_prime[m, :]).sum() / self.HPHR[n]
karolina's avatar
karolina committed
357
358
                self.Hx[m] = self.Hx[m] + fac * res
                self.HX_prime[m, :] = self.HX_prime[m, :] - alpha * fac * self.HX_prime[n, :]
359

karolina's avatar
karolina committed
360
361
            for m in range(1, n + 1):
                res = self.obs[n] - self.Hx[n]
362
                fac = 1.0 / (self.nmembers - 1) * (self.HX_prime[n, :] * self.HX_prime[m, :]).sum() / self.HPHR[n]
karolina's avatar
karolina committed
363
364
                self.Hx[m] = self.Hx[m] + fac * res
                self.HX_prime[m, :] = self.HX_prime[m, :] - alpha * fac * self.HX_prime[n, :]
365
366

            
367
    def bulk_minimum_least_squares(self):
368
        """ Make minimum least squares solution by solving matrix equations"""
karolina's avatar
karolina committed
369
        
370
371
372

        # Create full solution, first calculate the mean of the posterior analysis

karolina's avatar
karolina committed
373
374
375
376
        HPH = np.dot(self.HX_prime, np.transpose(self.HX_prime)) / (self.nmembers - 1)   # HPH = 1/N * HX' * (HX')^T
        self.HPHR[:, :] = HPH + self.R                                                            # HPHR = HPH + R
        HPb = np.dot(self.X_prime, np.transpose(self.HX_prime)) / (self.nmembers - 1)    # HP = 1/N X' * (HX')^T
        self.KG[:, :] = np.dot(HPb, la.inv(self.HPHR))                                         # K = HP/(HPH+R)
377
378

        for n in range(self.nobs):
379
            self.localize(n)
380

karolina's avatar
karolina committed
381
        self.x[:] = self.x + np.dot(self.KG, self.obs - self.Hx)                             # xa = xp + K (y-Hx)
382
383
384
385
386

        # And next make the updated ensemble deviations. Note that we calculate P by using the full equation (10) at once, and 
        # not in a serial update fashion as described in Whitaker and Hamill. 
        # For the current problem with limited N_obs this is easier, or at least more straightforward to do.

karolina's avatar
karolina committed
387
388
389
390
391
392
        I = np.identity(self.nlag * self.nparams)
        sHPHR = la.cholesky(self.HPHR)                                  # square root of HPH+R
        part1 = np.dot(HPb, np.transpose(la.inv(sHPHR)))                 # HP(sqrt(HPH+R))^-1
        part2 = la.inv(sHPHR + np.sqrt(self.R))                           # (sqrt(HPH+R)+sqrt(R))^-1
        Kw = np.dot(part1, part2)                                     # K~
        self.X_prime[:, :] = np.dot(I, self.X_prime) - np.dot(Kw, self.HX_prime)         # HX' = I - K~ * HX'
393
394


395
        # Now do the adjustments of the modeled mole fractions using the linearized ensemble. These are not strictly needed but can be used
396
397
        # for diagnosis.

karolina's avatar
karolina committed
398
399
400
401
        part3 = np.dot(HPH, np.transpose(la.inv(sHPHR)))                           # HPH(sqrt(HPH+R))^-1
        Kw = np.dot(part3, part2)                                               # K~
        self.Hx[:] = self.Hx + np.dot(np.dot(HPH, la.inv(self.HPHR)), self.obs - self.Hx)  # Hx  = Hx+ HPH/HPH+R (y-Hx)
        self.HX_prime[:, :] = self.HX_prime - np.dot(Kw, self.HX_prime)                            # HX' = HX'- K~ * HX'
402

karolina's avatar
karolina committed
403
        logging.info('Minimum Least Squares solution was calculated, returning')
404

405
406

    def set_localization(self, loctype='None'):
407
        """ determine which localization to use """
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425

        if loctype == 'CT2007':
            self.localization = True
            self.localizetype = 'CT2007'
            #T-test values for two-tailed student's T-test using 95% confidence interval for some options of nmembers
            if self.nmembers == 50:
                self.tvalue = 2.0086
            elif self.nmembers == 100:
                self.tvalue = 1.9840
            elif self.nmembers == 150:
                self.tvalue = 1.97591
            elif self.nmembers == 200:
                self.tvalue = 1.9719    
            else: self.tvalue = 0    
        else:
            self.localization = False
            self.localizetype = 'None'
    
karolina's avatar
karolina committed
426
        logging.info("Current localization option is set to %s" % self.localizetype)
427
428
429
430
431
        if self.localization == True:
            if self.tvalue == 0:
                logging.error("Critical tvalue for localization not set for %i ensemble members"%(self.nmembers))
                sys.exit(2)
            else: logging.info("Used critical tvalue %0.05f is based on 95%% probability and %i ensemble members in a two-tailed student's T-test"%(self.tvalue,self.nmembers))
432

433
    def localize(self, n):
434
435
        """ localize the Kalman Gain matrix """
    
436
437
438
439
440
441
442
443
444
    def set_algorithm(self, algorithm='Serial'):
        """ determine which minimum least squares algorithm to use """

        if algorithm == 'Serial':
            self.algorithm = 'Serial'
        else:
            self.algorithm = 'Bulk'
    
        logging.info("Current minimum least squares algorithm is set to %s" % self.algorithm)
445
446
447
448
449
450

################### End Class Optimizer ###################



if __name__ == "__main__":
451
    pass