optimizer.py 3.36 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#!/usr/bin/env python
# optimizer.py

"""
Author : peters 

Revision History:
File created on 28 Jul 2010.

"""

import os
import sys
import logging
karolina's avatar
karolina committed
15
sys.path.append(os.getcwd())
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
from da.baseclasses.optimizer import Optimizer


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

################### Begin Class CtOptimizer ###################

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

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

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

karolina's avatar
karolina committed
35
        if loctype == 'CT2007':
36
37
            self.localization = True
            self.localizetype = 'CT2007'
38
39
40
41
42
43
44
45
46
47
            #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    
48
49
50
        else:
            self.localization = False
            self.localizetype = 'None'
51
            self.tvalue = 0
52
    
karolina's avatar
karolina committed
53
        logging.info("Current localization option is set to %s" % self.localizetype)
54
55
56
        if self.tvalue == 0:
            logging.error("Critical tvalue for localization not set for %i ensemble members"%(self.nmembers))
            sys.exit(2)
57
        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))
58

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

karolina's avatar
karolina committed
63
        if not self.localization: 
64
            logging.debug('Not localized observation %i' % self.obs_ids[n])
karolina's avatar
karolina committed
65
            return 
66
        if self.localizetype == 'CT2007':
67
68
69
            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]
70
                prob = corr / np.sqrt((1.000000001 - corr ** 2) / (self.nmembers - 2))
71
72
73
74
                if abs(prob) < self.tvalue:
                    self.KG[r, n] = 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)))
75
76
77
78
79
80
81


################### End Class CtOptimizer ###################

if __name__ == "__main__":

    from da.tools.initexit import StartLogger 
karolina's avatar
karolina committed
82
    from da.tools.pipeline import start_job 
83
84
85
86

    sys.path.append(os.getcwd())

    opts = ['-v']
karolina's avatar
karolina committed
87
    args = {'rc':'da.rc', 'logfile':'da_initexit.log', 'jobrcfilename':'test.rc'}
88
89
90

    StartLogger()

karolina's avatar
karolina committed
91
    DaCycle = start_job(opts, args)
92
93
94
95
    DaCycle.Initialize()

    opt = CtOptimizer()

karolina's avatar
karolina committed
96
97
    nobs = 100
    dims = (int(DaCycle['time.nlag']),
98
                  int(DaCycle['da.optimizer.nmembers']),
99
                  int(DaCycle.DaSystem['nparameters']),
karolina's avatar
karolina committed
100
                  nobs,)
101
102

    opt.Initialize(dims)
103
    opt.set_localization(type='CT2007')