dynamic_operator.py 7.52 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
# 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, either version 3 of the License, or
# (at your option) any later version.
#
# 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/>.
#
Martin Reinecke's avatar
Martin Reinecke committed
14
# Copyright(C) 2013-2019 Max-Planck-Society
15
#
Martin Reinecke's avatar
Martin Reinecke committed
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
17
18
19

import numpy as np

Philipp Frank's avatar
Philipp Frank committed
20
from ..domain_tuple import DomainTuple
Philipp Arras's avatar
Philipp Arras committed
21
from ..domains.unstructured_domain import UnstructuredDomain
Philipp Frank's avatar
Philipp Frank committed
22
from ..domains.rg_space import RGSpace
Philipp Arras's avatar
Philipp Arras committed
23
from ..field import Field
Philipp Frank's avatar
cleanup  
Philipp Frank committed
24
from ..operators.diagonal_operator import DiagonalOperator
Philipp Arras's avatar
Philipp Arras committed
25
26
27
28
from ..operators.field_zero_padder import FieldZeroPadder
from ..operators.harmonic_operators import FFTOperator
from ..operators.scaling_operator import ScalingOperator
from ..operators.simple_linear_operators import FieldAdapter, Realizer
29
from ..sugar import makeOp
Philipp Arras's avatar
Fixup  
Philipp Arras committed
30
from .light_cone_operator import LightConeOperator, _field_from_function
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
31
32


Philipp Frank's avatar
Philipp Frank committed
33
def _make_dynamic_operator(domain, harmonic_padding, sm_s0, sm_x0, keys=['f', 'c'],
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
34
35
                          causal=True, cone=True, minimum_phase=False, sigc=3.,
                          quant=5.):
Philipp Frank's avatar
Philipp Frank committed
36
37
    dom = DomainTuple.make(domain)
    if not isinstance(dom[0], RGSpace):
Philipp Frank's avatar
Philipp Frank committed
38
        raise TypeError("RGSpace required")
39
    ops = {}
Philipp Frank's avatar
Philipp Frank committed
40
41
    FFT = FFTOperator(dom)
    Real = Realizer(dom)
Philipp Frank's avatar
Philipp Frank committed
42
43
    ops['FFT'] = FFT
    ops['Real'] = Real
44
    if harmonic_padding is None:
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
45
        CentralPadd = ScalingOperator(1., FFT.target)
46
    else:
Philipp Frank's avatar
Philipp Frank committed
47
48
49
50
        if isinstance(harmonic_padding, int):
            harmonic_padding = list((harmonic_padding,)*len(FFT.target.shape))
        elif len(harmonic_padding) != len(FFT.target.shape):
            raise (ValueError, "Shape mismatch!")
51
52
53
        shp = ()
        for i in range(len(FFT.target.shape)):
            shp += (FFT.target.shape[i] + harmonic_padding[i],)
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
54
        CentralPadd = FieldZeroPadder(FFT.target, shp, central=True)
Philipp Frank's avatar
Philipp Frank committed
55
    ops['central_padding'] = CentralPadd
56
57
    sdom = CentralPadd.target[0].get_default_codomain()
    FFTB = FFTOperator(sdom)(Realizer(sdom))
Philipp Frank's avatar
Philipp Frank committed
58

59
    m = FieldAdapter(sdom, keys[0])
Philipp Frank's avatar
Philipp Frank committed
60

61
    dists = m.target[0].distances
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
62
    if isinstance(sm_x0, float):
Philipp Frank's avatar
Philipp Frank committed
63
        sm_x0 = list((sm_x0,)*len(dists))
Philipp Frank's avatar
cleanup  
Philipp Frank committed
64
    elif len(sm_x0) != len(dists):
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
65
66
        raise (ValueError, "Shape mismatch!")

Philipp Frank's avatar
cleanup  
Philipp Frank committed
67
    def smoothness_prior_func(x):
68
        res = 1.
Philipp Frank's avatar
Philipp Frank committed
69
        for i in range(len(dists)):
70
71
            res = res + (x[i]/sm_x0[i]/dists[i])**2
        return sm_s0/res
Philipp Frank's avatar
Philipp Frank committed
72

Philipp Arras's avatar
Cleanup  
Philipp Arras committed
73
    Sm = makeOp(_field_from_function(m.target, smoothness_prior_func))
Philipp Frank's avatar
cleanup  
Philipp Frank committed
74
    m = CentralPadd.adjoint(FFTB(Sm(m)))
Philipp Frank's avatar
Philipp Frank committed
75
    ops['smoothed_dynamics'] = m
Philipp Frank's avatar
Philipp Frank committed
76

77
78
79
    m = -m.log()
    if not minimum_phase:
        m = m.exp()
Philipp Frank's avatar
Philipp Frank committed
80
81
    if causal or minimum_phase:
        m = Real.adjoint(FFT.inverse(Realizer(FFT.target).adjoint(m)))
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
82
83
        kernel = makeOp(
            _field_from_function(FFT.domain, (lambda x: 1. + np.sign(x[0]))))
84
        m = kernel(m)
Philipp Frank's avatar
Philipp Frank committed
85
86

    if cone and len(m.target.shape) > 1:
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
87
88
89
90
        if isinstance(sigc, float):
            sigc = list((sigc,)*(len(m.target.shape) - 1))
        elif len(sigc) != len(m.target.shape) - 1:
            raise (ValueError, "Shape mismatch!")
Philipp Frank's avatar
cleanup  
Philipp Frank committed
91
92
93
        c = FieldAdapter(UnstructuredDomain(len(sigc)), keys[1])
        c = makeOp(Field.from_global_data(c.target, np.array(sigc)))(c)

Philipp Arras's avatar
Cleanup  
Philipp Arras committed
94
        lightspeed = ScalingOperator(-0.5, c.target)(c).exp()
Philipp Frank's avatar
hotfix  
Philipp Frank committed
95
        scaling = np.array(m.target[0].distances[1:])/m.target[0].distances[0]
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
96
        scaling = DiagonalOperator(Field.from_global_data(c.target, scaling))
Philipp Frank's avatar
cleanup  
Philipp Frank committed
97
98
99
        ops['lightspeed'] = scaling(lightspeed)

        c = LightConeOperator(c.target, m.target, quant)(c.exp())
Philipp Frank's avatar
Philipp Frank committed
100
        ops['light_cone'] = c
101
102
        m = c*m

Philipp Frank's avatar
hotfix  
Philipp Frank committed
103
    if causal or minimum_phase:
Philipp Frank's avatar
Philipp Frank committed
104
        m = FFT(Real(m))
105
106
107
    if minimum_phase:
        m = m.exp()
    return m, ops
Philipp Frank's avatar
Philipp Frank committed
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196

def dynamic_operator(domain, harmonic_padding, sm_s0, sm_x0, key,
                     causal=True, minimum_phase=False):
    '''
    Constructs an operator encoding the Greens function of a linear homogeneous dynamic system.

    Parameters
    ----------
    domain : RGSpace
        The space under consideration
    harmonic_padding : None, int, list of int
        Amount of central padding in harmonic space in pixels. If None the field is not padded at all.
    sm_s0 : float
        Cutoff for dynamic smoothness prior
    sm_x0 : float, List of float
        Scaling of dynamic smoothness along each axis
    key : String
        key for dynamics encoding parameter.
    causal : boolean
        Whether or not the reconstructed dynamics should be causal in time
    minimum_phase: boolean
        Whether or not the reconstructed dynamics should be minimum phase

    Returns
    -------
    Operator
        The Operator encoding the dynamic Greens function in harmonic space.
    Dictionary of Operator
        A collection of sub-chains of Operators which can be used for plotting and evaluation.

    Notes
    -----
    Currently only supports RGSpaces.
    Note that the first axis of the space is interpreted as the time axis.
    '''
    return _make_dynamic_operator(domain, harmonic_padding, sm_s0, sm_x0,
                                  keys=[key],
                                  causal=causal, cone=False,
                                  minimum_phase=minimum_phase)

def dynamic_lightcone_operator(domain, harmonic_padding, sm_s0, sm_x0, key,
                               lightcone_key, sigc, quant,
                               causal=True, minimum_phase=False):
    '''
    Constructs an operator encoding the Greens function of a linear homogeneous dynamic system.
    The Greens function is constrained to be within a light cone.

    Parameters
    ----------
    domain : RGSpace
        The space under consideration. Must have dim > 1.
    harmonic_padding : None, int, list of int
        Amount of central padding in harmonic space in pixels. If None the field is not padded at all.
    sm_s0 : float
        Cutoff for dynamic smoothness prior
    sm_x0 : float, List of float
        Scaling of dynamic smoothness along each axis
    key : String
        key for dynamics encoding parameter.
    lightcone_key: String
        key for lightspeed paramteter.
    sigc : float, List of float
        variance of lightspeed parameter.
    quant : float
        Quantization of the light cone in pixels.
    causal : boolean
        Whether or not the reconstructed dynamics should be causal in time
    minimum_phase: boolean
        Whether or not the reconstructed dynamics should be minimum phase

    Returns
    -------
    Operator
        The Operator encoding the dynamic Greens function in harmonic space when evaluated.
    Dictionary of Operator
        A collection of sub-chains of Operators which can be used for plotting and evaluation.

    Notes
    -----
    Currently only supports RGSpaces.
    Note that the first axis of the space is interpreted as the time axis.
    '''
    if len(domain.shape) < 2:
        raise ValueError("Space must be at least 2 dimensional!")
    return _make_dynamic_operator(domain,harmonic_padding,sm_s0,sm_x0,
                                  keys=[key,lightcone_key],
                                  causal=causal, cone=True,
                                  minimum_phase = minimum_phase,
                                  sigc = sigc, quant = quant)