execute.py 3.56 KB
Newer Older
1
2
import numpy as np
import copy
Antoine Berchet's avatar
Antoine Berchet committed
3
from pycif.utils.check import verbose
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18


def execute(self, **kwargs):
    """Runs the model in forward mode

    Args:
        self (Plugin): the mode Plugin with all set-up arguments

    """
    
    # Working directory
    workdir = self.workdir
    
    # Control vector
    controlvect = self.controlvect
19

20
21
    # Observation operator
    obsoper = self.obsoperator
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
    # Simulation window
    datei = self.datei
    datef = self.datef
    
    # Increments in x
    increments = getattr(self, 'increments', 0.)
    incrmode = getattr(self, 'incrmode', 'cst')
    testspace = getattr(self, 'testspace', 'control')
    
    if testspace == 'control':
        if incrmode == 'cst':
            controlvect.dx = increments * controlvect.std
        
        elif incrmode == 'rand':
            controlvect.dx = np.random.normal(0, increments, controlvect.dim) \
                             * controlvect.std
    
    elif testspace == 'chi':
        if incrmode == 'cst':
            controlvect.chi[:] = increments
        
        elif incrmode == 'rand':
            controlvect.chi = \
                np.random.normal(0, increments, controlvect.chi_dim)
47
        
48
49
50
51
52
53
54
55
56
57
        controlvect.dx = \
            controlvect.sqrtbprod(controlvect.chi, **kwargs) - controlvect.xb
    
    if testspace not in ['control', 'chi'] \
            or incrmode not in ['cst', 'rand']:
        verbose("The increment mode you specified can't be parsed: {}"
                .format(incrmode))
        verbose("Please check the definition of the running mode "
                "in your Yaml file")
        raise Exception
58
    
59
60
    controlvect.x = copy.deepcopy(controlvect.xb)
    controlvect.xb += controlvect.dx
61

62
63
    # Some verbose
    verbose("Computing the test of the adjoint")
64
    
65
66
    # Get the machine accuracy
    accuracy = np.finfo(np.float64).eps
67

68
69
70
71
    # Running the tangent linear code of the model
    obsvect = obsoper.obsoper(controlvect, 'tl',
                              datei=datei, datef=datef,
                              workdir=workdir,
72
73
                              reload_results=getattr(self, 'reload_results',
                                                     False),
74
                              **kwargs)
75

76
    # Computing < H.dx, H.dx >
Espen Sollum's avatar
Espen Sollum committed
77
    scaleprod1 = (obsvect.datastore['sim_tl'] ** 2).sum()
78

79
80
81
82
    # Putting increments in the observation vector
    obsvect.datastore['obs_incr'] = obsvect.datastore['sim_tl']
    obsvect.datastore.loc[
        np.isnan(obsvect.datastore['obs_incr']), 'obs_incr'] = 0
83

84
85
86
87
    # Running the observation operator
    controlvect = obsoper.obsoper(obsvect, 'adj',
                                  datei=datei, datef=datef,
                                  workdir=workdir,
88
89
                                  reload_results=getattr(self, 'reload_results',
                                                         False),
90
                                  **kwargs)
91
    
92
93
    # Computing < dx, H*(H.dx) >
    if testspace == 'control':
Espen Sollum's avatar
Espen Sollum committed
94
95
        scaleprod2 = (controlvect.dx
                      * (controlvect.xb - controlvect.x)).sum()
96

97
98
99
100
    elif testspace == 'chi':
        scaleprod2 = \
            (controlvect.sqrtbprod_ad(controlvect.dx, **kwargs) *
             controlvect.chi).sum()
101

102
103
    else:
        scaleprod2 = np.nan
104

105
106
107
108
    # Final verbose
    verbose('Machine accuracy: {}'.format(accuracy))
    verbose('< H.dx, H.dx >   = {:.17E}'.format(scaleprod1))
    verbose('< dx, H*(H.dx) > = {:.17E}'.format(scaleprod2))
Espen Sollum's avatar
Espen Sollum committed
109
    verbose('The difference is {:.1E} time the machine accuracy'
110
            .format(np.abs(scaleprod2 / scaleprod1 - 1) / accuracy))
111