wiener_process_integrated_amplitude.py 3.7 KB
Newer Older
1
2
3
4
import nifty5 as ift
import numpy as np


Philipp Arras's avatar
Work    
Philipp Arras committed
5
class _TwoLogIntegrations(ift.LinearOperator):
6
7
8
9
10
    def __init__(self, target):
        self._target = ift.makeDomain(target)
        self._domain = ift.makeDomain(
            ift.UnstructuredDomain(self.target.shape[0] - 2))
        self._capability = self.TIMES | self.ADJOINT_TIMES
Philipp Arras's avatar
Work    
Philipp Arras committed
11
12
13
14
        if not isinstance(self._target[0], ift.PowerSpace):
            raise TypeError
        logk_lengths = np.log(self._target[0].k_lengths[1:])
        self._logvol = logk_lengths[1:] - logk_lengths[:-1]
15
16
17
18
19
20
21
22

    def apply(self, x, mode):
        self._check_input(x, mode)
        if mode == self.TIMES:
            x = x.to_global_data()
            res = np.empty(self._target.shape)
            res[0] = 0
            res[1] = 0
Philipp Arras's avatar
Work    
Philipp Arras committed
23
24
            res[2:] = np.cumsum(x*self._logvol)
            res[2:] = np.cumsum(res[2:]*self._logvol)
25
26
27
28
            return ift.from_global_data(self._target, res)
        else:
            x = x.to_global_data()
            res = np.empty(self._target.shape)
Philipp Arras's avatar
Work    
Philipp Arras committed
29
30
            res[2:] = np.cumsum(x[2:][::-1])[::-1]*self._logvol
            res[2:] = np.cumsum(res[2:][::-1])[::-1]*self._logvol
31
32
33
            return ift.from_global_data(self._domain, res[2:])


Philipp Arras's avatar
Work    
Philipp Arras committed
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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
class _Rest(ift.LinearOperator):
    def __init__(self, target):
        self._target = ift.makeDomain(target)
        self._domain = ift.makeDomain(ift.UnstructuredDomain(3))
        self._logk_lengths = np.log(self._target[0].k_lengths[1:])
        self._logk_lengths -= self._logk_lengths[0]
        self._capability = self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
        self._check_input(x, mode)
        x = x.to_global_data()
        res = np.empty(self._tgt(mode).shape)
        if mode == self.TIMES:
            res[0] = x[0]
            res[1:] = x[1]*self._logk_lengths + x[2]
        else:
            res[0] = x[0]
            res[1] = np.vdot(self._logk_lengths, x[1:])
            res[2] = np.sum(x[1:])
        return ift.from_global_data(self._tgt(mode), res)


def LogIntegratedWienerProcess(target, means, stddevs, keys):
    # means and stddevs: zm, slope, yintercept, wienersigma
    # keys: rest smooth wienersigma
    if not (len(means) == 4 and len(stddevs) == 4 and len(keys) == 3):
        raise ValueError
    means = np.array(means)
    stddevs = np.array(stddevs)
    # FIXME More checks
    rest = _Rest(target)
    restmeans = ift.from_global_data(rest.domain, means[:-1])
    reststddevs = ift.from_global_data(rest.domain, stddevs[:-1])
    rest = rest @ ift.Adder(restmeans) @ ift.makeOp(reststddevs)

    expander = ift.VdotOperator(ift.full(target, 1.)).adjoint
    sigma = ift.Adder(ift.full(expander.domain, means[3])) @ (
        stddevs[3]*ift.ducktape(expander.domain, None, keys[2]))
    sigma = expander @ sigma.exp()
    smooth = _TwoLogIntegrations(target).ducktape(keys[1])*sigma
    return rest.ducktape(keys[0]) + smooth


77
78
if __name__ == '__main__':
    np.random.seed(42)
Philipp Arras's avatar
Philipp Arras committed
79
    ndim = 2
Philipp Arras's avatar
Work    
Philipp Arras committed
80
    sspace = ift.RGSpace((512, 512))
81
    hspace = sspace.get_default_codomain()
Philipp Arras's avatar
Work    
Philipp Arras committed
82
83
84
85
86
87
88
89
90
91
    target = ift.PowerSpace(hspace,
                            ift.PowerSpace.useful_binbounds(hspace, True))
    test0 = _Rest(target)
    test1 = _TwoLogIntegrations(target)
    ift.extra.consistency_check(test0)
    ift.extra.consistency_check(test1)
    op = LogIntegratedWienerProcess(target,
                                    [0, -4, 1, 0],
                                    [1, 1, 1, 0.5],
                                    ['rest', 'smooth', 'wienersigma']).exp()
92
    fld = ift.from_random('normal', op.domain)
Philipp Arras's avatar
Work    
Philipp Arras committed
93
94
95
96
97
98
    ift.extra.check_jacobian_consistency(op, fld)
    plts = []
    for _ in range(50):
        fld = ift.from_random('normal', op.domain)
        plts.append(op(fld))
    ift.single_plot(plts, name='debug.png')