test_optimizer.py 7.28 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
17
18
19
20
21
22
23
#!/usr/bin/env python
# test_optimizer.py

"""
Author : peters 

Revision History:
File created on 04 Aug 2010.

"""

24
25
26
print("Ingrid is testing git here")


27
def serial_py_against_serial_fortran():
28
29
30
31
32
33
34
35
36
37
38
39
    """ Test the solution of the serial algorithm against the CT cy2 fortran generated one """

    # get data from the savestate.hdf file from the first cycle of CarbonTracker 2009 release

    print "WARNING: The optimization algorithm has changed from the CT2009 release because of a bug"
    print "WARNING: in the fortran code. Hence, the two solutions calculated are no longer the same."
    print "WARNING: To change the python algorithm so that it corresponds to the fortran, change the"
    print "WARNING: loop from m=n+1,nfcast to m=1,nfcast"

    savefile = '/data/CO2/peters/carbontracker/raw/ct09rc0i/20000108/savestate.hdf'
    print savefile

40
    f = Nio.open_file(savefile, 'r')
41
42
43
    obs = f.variables['co2_obs_fcast'].get_value()
    sel_obs = obs.shape[0]

44
45
46
    dims = (int(dacycle.da_settings['time.nlag']),
                  int(dacycle.da_settings['forecast.nmembers']),
                  int(dacycle.dasystem.da_settings['nparameters']),
47
                  sel_obs,)
48

49
    nlag, nmembers, nparams, nobs = dims
50
51
52
53

    optserial = CtOptimizer(dims)
    opt = optserial

54
    opt.set_localization('CT2007')
55
56
57
58
59
60
61
62
63
64

    obs = f.variables['co2_obs_fcast'].get_value()[0:nobs]
    opt.obs = obs
    sim = f.variables['co2_sim_fcast'].get_value()[0:nobs]
    opt.Hx = sim
    error = f.variables['error_sim_fcast'].get_value()[0:nobs]
    flags = f.variables['flag_sim_fcast'].get_value()[0:nobs]
    opt.flags = flags
    simana = f.variables['co2_sim_ana'].get_value()[0:nobs]

65
    for n in range(nobs): opt.R[n, n] = np.double(error[n] ** 2)
66

67
68
    xac = []
    adX = []
69
    for lag in range(nlag):
70
71
72
73
74
75
        xpc = f.variables['xpc_%02d' % (lag + 1)].get_value()
        opt.x[lag * nparams:(lag + 1) * nparams] = xpc
        X = f.variables['pdX_%02d' % (lag + 1)].get_value()
        opt.X_prime[lag * nparams:(lag + 1) * nparams, :] = np.transpose(X)
        HX = f.variables['dF'][:, 0:sel_obs]
        opt.HX_prime[:, :] = np.transpose(HX)
76
77
78

        # Also create arrays of the analysis of the fortran code for later comparison

79
80
        xac.extend (f.variables['xac_%02d' % (lag + 1)].get_value())
        adX.append (f.variables['adX_%02d' % (lag + 1)].get_value())
81

82
83
    xac = np.array(xac)
    X_prime = np.array(adX).swapaxes(1, 2).reshape((opt.nparams * opt.nlag, opt.nmembers))
84

85
    opt.serial_minimum_least_squares()
86
87

    print "Maximum differences and correlation of 2 state vectors:"
88
    print np.abs(xac - opt.x).max(), np.corrcoef(xac, opt.x)[0, 1]
89
90
       
    plt.figure(1)
91
92
    plt.plot(opt.x, label='SerialPy')
    plt.plot(xac, label='SerialFortran')
93
94
95
96
97
    plt.grid(True)
    plt.legend(loc=0)
    plt.title('Analysis of state vector')

    print "Maximum differences of 2 state vector deviations:"
98
    print np.abs(X_prime - opt.X_prime).max()
99
100

    plt.figure(2)
101
102
    plt.plot(opt.X_prime.flatten(), label='SerialPy')
    plt.plot(X_prime.flatten(), label='SerialFortran')
103
104
105
106
107
    plt.grid(True)
    plt.legend(loc=0)
    plt.title('Analysis of state vector deviations')

    print "Maximum differences and correlation of 2 simulated obs vectors:"
108
    print np.abs(simana - opt.Hx).max(), np.corrcoef(simana, opt.Hx)[0, 1]
109
110

    plt.figure(3)
111
112
    plt.plot(opt.Hx, label='SerialPy')
    plt.plot(simana, label='SerialFortran')
113
114
    plt.grid(True)
    plt.legend(loc=0)
115
    plt.title('Analysis of CO2 mole fractions')
116
117
118
119
    plt.show()

    f.close()

120
def serial_vs_bulk():
121
122
123
124
125
126
127
    """ A test of the two algorithms currently implemented: serial vs bulk solution """    

    # get data from the savestate.hdf file from the first cycle of CarbonTracker 2009 release

    savefile = '/data/CO2/peters/carbontracker/raw/ct09rc0i/20000108/savestate.hdf'
    print savefile

128
    f = Nio.open_file(savefile, 'r')
129
130
131
132
    obs = f.variables['co2_obs_fcast'].get_value()

    nobs = 77

133
134
135
    dims = (int(dacycle.da_settings['time.nlag']),
                  int(dacycle.da_settings['forecast.nmembers']),
                  int(dacycle.dasystem.da_settings['nparameters']),
136
                  nobs,)
137

138
    nlag, nmembers, nparams, nobs = dims
139

140
    optbulk = CtOptimizer(dims)
141
142
    optserial = CtOptimizer(dims)

143
    for o, opt in enumerate([optbulk, optserial]):
144

145
        opt.set_localization('CT2007')
146
147
148
149
150
151
152
153
154
155

        obs = f.variables['co2_obs_fcast'].get_value()[0:nobs]
        opt.obs = obs
        sim = f.variables['co2_sim_fcast'].get_value()[0:nobs]
        opt.Hx = sim
        error = f.variables['error_sim_fcast'].get_value()[0:nobs]
        flags = f.variables['flag_sim_fcast'].get_value()[0:nobs]
        opt.flags = flags

        for n in range(nobs): 
156
            opt.R[n, n] = np.double(error[n] ** 2)
157

158
        xac = []
159
        for lag in range(nlag):
160
161
162
163
164
165
            xpc = f.variables['xpc_%02d' % (lag + 1)].get_value()
            opt.x[lag * nparams:(lag + 1) * nparams] = xpc
            X = f.variables['pdX_%02d' % (lag + 1)].get_value()
            opt.X_prime[lag * nparams:(lag + 1) * nparams, :] = np.transpose(X)
            HX = f.variables['dF'][:, 0:nobs]
            opt.HX_prime[:, :] = np.transpose(HX)
166
167

        if o == 0:
168
            opt.bulk_minimum_least_squares()
169
170
171
172
173
174
            x1 = opt.x
            xp1 = opt.X_prime
            hx1 = opt.Hx
            hxp1 = opt.HX_prime
            hphr1 = opt.HPHR
            k1 = opt.KG
175
        if o == 1:
176
            opt.serial_minimum_least_squares()
177
178
179
180
181
182
            x2 = opt.x
            xp2 = opt.X_prime
            hx2 = opt.Hx
            hxp2 = opt.HX_prime
            hphr2 = opt.HPHR
            k2 = opt.KG
183
184
185
186
           
    plt.figure()

    print "Maximum differences and correlation of 2 state vectors:"
187
    print np.abs(x2 - x1).max(), np.corrcoef(x2, x1)[0, 1]
188
189
       
    plt.figure(1)
190
191
    plt.plot(x1, label='Serial')
    plt.plot(x2, label='Bulk')
192
193
194
195
196
    plt.grid(True)
    plt.legend(loc=0)
    plt.title('Analysis of state vector')

    print "Maximum differences of 2 state vector deviations:"
197
    print np.abs(xp2 - xp1).max()
198
199

    plt.figure(2)
200
201
    plt.plot(xp1.flatten(), label='Serial')
    plt.plot(xp2.flatten(), label='Bulk')
202
203
204
205
206
    plt.grid(True)
    plt.legend(loc=0)
    plt.title('Analysis of state vector deviations')

    print "Maximum differences and correlation of 2 simulated obs vectors:"
207
    print np.abs(hx2 - hx1).max(), np.corrcoef(hx2, hx1)[0, 1]
208
209

    plt.figure(3)
210
211
    plt.plot(hx1, label='Serial')
    plt.plot(hx2, label='Bulk')
212
    plt.title('Analysis of CO2 mole fractions')
213
214
215
216
217
218
219
220
221
222
    plt.grid(True)
    plt.legend(loc=0)

    plt.show()

    f.close()



if __name__ == "__main__":
223
    pass