obsoper.py 8.69 KB
Newer Older
1
2
3
4
import pandas as pd
import datetime
import numpy as np

Antoine Berchet's avatar
Antoine Berchet committed
5
6
7
from pycif.utils import path
from pycif.utils.check import verbose
from pycif.utils.datastores.dump import dump_datastore, read_datastore
8
9
from check import check_inputs

10
import pdb
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

def obsoper(self, inputs, mode,
            run_id=0, datei=datetime.datetime(1979, 1, 1),
            datef=datetime.datetime(2100, 1, 1),
            workdir='./',
            reload_results=False,
            **kwargs):
    """The observation operator.
    This function maps information from the control space to the observation
    space and conversely depending on the running mode.
    
    Generates model inputs from the control vector
    inputs(x)
    
    Turns observations into model compatible extraction points
    i.e. y0 (fwd) or dy* (adj)
    
    Runs the model to get model extractions
    i.e. Hx (fwd) or H^T.dy* (adj)
    
    if FWD
    Turns back model-equivalents into observation space
    Hx(native) -> Hx(obs)

    Generates departures and multiply by R-1
    -> Hx - y0, R^(-1).(Hx - y0)
    
    if ADJ
    Turns native increments to control space increments
    H^T.dy* (native) -> H^T.dy* (control)

    Translates control increments to chi values
    -> B^(1/2) . H^T.dy*
    
    Args:
        inputs: can be a control vector or an observation vector, depending
                on the mode
        mode (str): the running mode; should be one of 'fwd', 'tl' or 'adj'
        run_id (int): the ID number of the current run; is used to define
                the current sub-directory
        datei, datef (datetime.datetime): start and end dates of the
                simulation window
        workdir (str): the parent directory of the simulation to run
        reload_results (bool): look for results from pre-computed simulation
                if any
    
    Returns:
        observation or control vector depending on the running mode
    """
60

61
62
63
64
65
66
67
68
69
    # Check that inputs are consistent with the mode to run
    check_inputs(inputs, mode)
    
    # Create of sub- working directory for the present run
    rundir = "{}/obsoperator/{}_{:04d}/".format(workdir, mode, run_id)
    path.init_dir(rundir)
    
    # Create save directory for chaining sub-simulations
    path.init_dir(
70
        "{}/chain/".format(rundir, run_id))
71
72
73
74
75
    
    # Return results from previous versions if existing
    if reload_results:
        if mode in ['fwd', 'tl']:
            try:
76
77
78
79
                # Saving the directory for possible later use by the adjoint
                self.model.adj_refdir = rundir
                
                # Trying reading the monitor if any
80
                obsvect = self.obsvect
81
82
                obsvect.datastore = \
                    read_datastore("{}/monitor.nc".format(rundir))
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
                return obsvect
            
            except IOError:
                verbose("There is no monitor file to be recovered. "
                        "Compute the full forward simulation")
        
        elif mode == 'adj':
            try:
                controlvect = self.controlvect
                controlvect.load("{}/controlvect.pickle".format(rundir))
                return controlvect
            
            except IOError:
                verbose("There is no controlvect file to be recovered. "
                        "Compute the full adjoint simulation.")
98
    
99
100
101
102
    # Initializing modules and variables from the setup
    model = self.model
    controlvect = self.controlvect
    obsvect = self.obsvect
103
    subsimu_dates = model.subsimu_dates
104
105
106
107

    print "obsoper/model", model
    print "obsoper/mode", mode
    pdb.set_trace()
108
109
110
    
    if mode in ['fwd', 'tl']:
        controlvect = inputs
111
112
        obsvect.datastore.loc[:, 'sim'] = 0.
        obsvect.datastore.loc[:, 'sim_tl'] = 0.
113
        
114
        # Dumps control vector in forward and tl modes
115
116
117
118
        controlvect.dump(
            "{}/controlvect.pickle".format(rundir),
            to_netcdf=getattr(controlvect, 'save_out_netcdf', False),
            dir_netcdf='{}/controlvect/'.format(rundir))
119
120
121
122
    
    elif mode == 'adj':
        obsvect = inputs
        controlvect.dx = 0 * controlvect.x
123
        subsimu_dates = subsimu_dates[::-1]
124

125
126
127
128
129
130
131
132
133
134
135
136
137
138

    #import code
    #code.interact(local=dict(locals(), **globals()))


    elif mode == 'footprint':
        print "obsoper; footprint"
        obsvect = inputs
        controlvect.dx = 0 * controlvect.x
        subsimu_dates = subsimu_dates[::-1]
        read_only = getattr(self, 'read_only', False)
        print "read_only", read_only
    

139
140
    # Loop through model periods and runs the model
    self.missingperiod = False
141
    for di, df in zip(subsimu_dates[:-1], subsimu_dates[1:]):
142
143
144
        # Create a sub-path for each period
        runsubdir = rundir + min(di, df).strftime("%Y-%m-%d_%H-%M")
        _, created = path.init_dir(runsubdir)
145
        # ESO: copy data to above directory
146
147
148
149
150
151
        
        # If the sub-directory was already created,
        # the observation operator considers that the sub-simulation
        # was already properly run, thus passing to next sub-periods
        # Compute the sub-simulation anyway if some previous periods
        # were missing (as stored in self.missingperiod)
Espen Sollum's avatar
Espen Sollum committed
152
        do_simu = created \
153
                  or not getattr(self, 'autorestart', False) \
Espen Sollum's avatar
Espen Sollum committed
154
                  or self.missingperiod
155
156
157
158
159
160
161
        self.missingperiod = do_simu
        
        # Some verbose
        verbose("Running {} for the period {} to {}"
                .format(model.plugin.name, di, df))
        verbose("Running mode: {}".format(mode))
        verbose("Sub-directory: {}".format(runsubdir))
162
163
        print "do_simu", do_simu
        print "read_only", read_only
164
165
        
        # Prepare observations for the model
166
167
        if not read_only:
            model.dataobs = obsvect.obsvect2native(di, df, mode, runsubdir, workdir)
168
169
170
171
172
        
        # If the simulation was already carried out, pass to next steps
        # If a sub-period was missing, following ones will be recomputed even
        # if available
        if do_simu:
173
174
175
176
            
            # Writing observations for on-line model extraction if any
            model.native2inputs(model.dataobs, 'obs', di, df, mode,
                                runsubdir, workdir)
177
            
178
179
180
            # Prepare inputs for the model
            for mod_input in model.required_inputs:
                verbose("Preparing {}".format(mod_input))
181
                
182
183
184
185
186
                # Turns state vector variables at the model resolution
                data = \
                    controlvect.control2native(
                        mod_input, di, df, mode,
                        runsubdir, workdir)
187
                
188
189
190
191
192
193
194
195
196
                # Create input files from variables at native resolution
                model.native2inputs(data, mod_input, di, df, mode,
                                    runsubdir, workdir)
        
        # If only initializing inputs, continue to next sub-period
        if getattr(self, 'onlyinit', False):
            continue
        
        # Run the model
197
198
199
        if not read_only:
            model.run(runsubdir, mode, workdir,
                      do_simu=do_simu, **kwargs)
200
        
201
202
203
204
205
        # Read outputs
        model.outputs2native(runsubdir, mode)
        
        # Update observation vector if necessary
        if mode in ['fwd', 'tl'] and obsvect.datastore.size > 0:
206
207
            obsvect.native2obsvect(model.dataobs, di, df, workdir, runsubdir)
        
208
209
210
211
212
213
214
        # Update state vector if necessary
        elif mode == 'adj':
            controlvect.native2control(
                model.datasensit, di, df, workdir, runsubdir)
        
        # Keep in memory the fact that it is (or not) a chained simulation
        model.chain = min(di, df)
215
    
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
    # If only initializing inputs, exit
    if getattr(self, 'onlyinit', False):
        verbose("The run was correctly initialized")
        return
    
    # Re-initalizing the chain argument
    if hasattr(model, 'chain'):
        del model.chain
    
    # Dump observation vector for later use in fwd and tl modes
    # Otherwise dumps the control vector
    if mode in ['fwd', 'tl']:
        dump_type = obsvect.dump_type
        dump_datastore(obsvect.datastore,
                       fic_monit='{}/monitor.{}'.format(rundir, dump_type),
                       mode='w', dump_type=dump_type)
232
    
233
234
235
236
237
238
239
240
241
242
243
244
245
246
    elif mode == 'adj':
        controlvect.dump("{}/controlvect.pickle".format(rundir))
    
    # Cleaning unnecessary files
    if getattr(model, 'autoflush', False):
        verbose("Flushing unnecessary files in {}".format(rundir))
        model.flushrun(rundir, mode)
    
    # Returning the output object depending on the running mode
    if mode in ['fwd', 'tl']:
        return obsvect
    
    if mode == 'adj':
        return controlvect
247
248
249

    if mode == 'footprint':
        return controlvect