optimizer.py 3.99 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
"""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
# optimizer.py

"""
Author : peters 

Revision History:
File created on 28 Jul 2010.

"""

import os
import sys
import logging
sys.path.append(os.getcwd())
from da.cosmo.base_optimizer import Optimizer


identifier = 'Ensemble Square Root Filter'
version = '0.0'

################### Begin Class CO2Optimizer ###################

class CO2Optimizer(Optimizer):
    """
        This creates an instance of a CarbonTracker optimization object. The base class it derives from is the optimizer object.
        Additionally, this CO2Optimizer implements a special localization option following the CT2007 method.

        All other methods are inherited from the base class Optimizer.
    """

    def set_localization(self, loctype='None'):
        """ determine which localization to use """

        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
51
52
53
            if self.nmembers == 21:
                self.tvalue = 2.08
            elif self.nmembers == 50:
brunner's avatar
brunner committed
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
                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'
    
        logging.info("Current localization option is set to %s" % self.localizetype)
        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))

    def localize(self, n):
        """ localize the Kalman Gain matrix """
        import numpy as np

        if not self.localization: 
            logging.debug('Not localized observation %i' % self.obs_ids[n])
            return 
        if self.localizetype == 'CT2007':
            count_localized = 0
            for r in range(self.nlag * self.nparams):
                corr = np.corrcoef(self.HX_prime[n, :], self.X_prime[r, :].squeeze())[0, 1]
                prob = corr / np.sqrt((1.000000001 - corr ** 2) / (self.nmembers - 2))
                if abs(prob) < self.tvalue:
                    self.KG[r] = 0.0
                    count_localized = count_localized + 1
            logging.debug('Localized observation %i, %i%% of values set to 0' % (self.obs_ids[n],count_localized*100/(self.nlag * self.nparams)))

    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)

################### End Class CO2Optimizer ###################

if __name__ == "__main__":
    pass