dynamic_operator.py 8.68 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 Frank's avatar
Philipp Frank committed
21
from ..domains.rg_space import RGSpace
Philipp Arras's avatar
Philipp Arras committed
22
from ..domains.unstructured_domain import UnstructuredDomain
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 Arras's avatar
Philipp Arras committed
33
34
35
36
def _float_or_listoffloat(inp):
    return [float(x) for x in inp] if isinstance(inp, list) else float(inp)


Philipp Arras's avatar
Philipp Arras committed
37
38
39
40
41
42
43
44
45
46
def _make_dynamic_operator(domain,
                           harmonic_padding,
                           sm_s0,
                           sm_x0,
                           cone,
                           keys,
                           causal,
                           minimum_phase,
                           sigc=None,
                           quant=None):
Philipp Arras's avatar
Philipp Arras committed
47
48
49
50
51
52
53
54
55
56
57
58
59
    if not isinstance(domain, RGSpace):
        raise TypeError("RGSpace required")
    if not (isinstance(harmonic_padding, int) or harmonic_padding is None
            or all(isinstance(ii, int) for ii in harmonic_padding)):
        raise TypeError
    sm_s0 = float(sm_s0)
    sm_x0 = _float_or_listoffloat(sm_x0)
    cone = bool(cone)
    if not all(isinstance(ss, str) for ss in keys):
        raise TypeError
    causal, minimum_phase = bool(causal), bool(minimum_phase)
    if sigc is not None:
        sigc = _float_or_listoffloat(sigc)
Philipp Arras's avatar
Philipp Arras committed
60
61
    if quant is not None:
        quant = float(quant)
Philipp Arras's avatar
Philipp Arras committed
62
63
64
    if cone and (sigc is None or quant is None):
        raise RuntimeError

Philipp Frank's avatar
Philipp Frank committed
65
    dom = DomainTuple.make(domain)
66
    ops = {}
Philipp Frank's avatar
Philipp Frank committed
67
68
    FFT = FFTOperator(dom)
    Real = Realizer(dom)
Philipp Frank's avatar
Philipp Frank committed
69
70
    ops['FFT'] = FFT
    ops['Real'] = Real
71
    if harmonic_padding is None:
Philipp Arras's avatar
Cleanup    
Philipp Arras committed
72
        CentralPadd = ScalingOperator(1., FFT.target)
73
    else:
Philipp Frank's avatar
Philipp Frank committed
74
75
76
77
        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!")
Philipp Arras's avatar
Philipp Arras committed
78
79
80
        shp = [
            sh + harmonic_padding[ii] for ii, sh in enumerate(FFT.target.shape)
        ]
Philipp Arras's avatar
Cleanup    
Philipp Arras committed
81
        CentralPadd = FieldZeroPadder(FFT.target, shp, central=True)
Philipp Frank's avatar
Philipp Frank committed
82
    ops['central_padding'] = CentralPadd
83
84
    sdom = CentralPadd.target[0].get_default_codomain()
    FFTB = FFTOperator(sdom)(Realizer(sdom))
Philipp Frank's avatar
Philipp Frank committed
85

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

88
    dists = m.target[0].distances
Philipp Arras's avatar
Cleanup    
Philipp Arras committed
89
    if isinstance(sm_x0, float):
Philipp Frank's avatar
Philipp Frank committed
90
        sm_x0 = list((sm_x0,)*len(dists))
Philipp Frank's avatar
cleanup    
Philipp Frank committed
91
    elif len(sm_x0) != len(dists):
Philipp Arras's avatar
Cleanup    
Philipp Arras committed
92
93
        raise (ValueError, "Shape mismatch!")

Philipp Frank's avatar
cleanup    
Philipp Frank committed
94
    def smoothness_prior_func(x):
95
        res = 1.
Philipp Frank's avatar
Philipp Frank committed
96
        for i in range(len(dists)):
97
98
            res = res + (x[i]/sm_x0[i]/dists[i])**2
        return sm_s0/res
Philipp Frank's avatar
Philipp Frank committed
99

Philipp Arras's avatar
Cleanup    
Philipp Arras committed
100
    Sm = makeOp(_field_from_function(m.target, smoothness_prior_func))
Philipp Frank's avatar
cleanup    
Philipp Frank committed
101
    m = CentralPadd.adjoint(FFTB(Sm(m)))
Philipp Frank's avatar
Philipp Frank committed
102
    ops['smoothed_dynamics'] = m
Philipp Frank's avatar
Philipp Frank committed
103

104
105
106
    m = -m.log()
    if not minimum_phase:
        m = m.exp()
Philipp Frank's avatar
Philipp Frank committed
107
108
    if causal or minimum_phase:
        m = Real.adjoint(FFT.inverse(Realizer(FFT.target).adjoint(m)))
Philipp Arras's avatar
Cleanup    
Philipp Arras committed
109
110
        kernel = makeOp(
            _field_from_function(FFT.domain, (lambda x: 1. + np.sign(x[0]))))
111
        m = kernel(m)
Philipp Frank's avatar
Philipp Frank committed
112
113

    if cone and len(m.target.shape) > 1:
Philipp Arras's avatar
Cleanup    
Philipp Arras committed
114
115
116
117
        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
118
119
120
        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
121
        lightspeed = ScalingOperator(-0.5, c.target)(c).exp()
Philipp Frank's avatar
hotfix    
Philipp Frank committed
122
        scaling = np.array(m.target[0].distances[1:])/m.target[0].distances[0]
Philipp Arras's avatar
Cleanup    
Philipp Arras committed
123
        scaling = DiagonalOperator(Field.from_global_data(c.target, scaling))
Philipp Frank's avatar
cleanup    
Philipp Frank committed
124
125
126
        ops['lightspeed'] = scaling(lightspeed)

        c = LightConeOperator(c.target, m.target, quant)(c.exp())
Philipp Frank's avatar
Philipp Frank committed
127
        ops['light_cone'] = c
128
129
        m = c*m

Philipp Frank's avatar
hotfix    
Philipp Frank committed
130
    if causal or minimum_phase:
Philipp Frank's avatar
Philipp Frank committed
131
        m = FFT(Real(m))
132
133
134
    if minimum_phase:
        m = m.exp()
    return m, ops
Philipp Frank's avatar
Philipp Frank committed
135

Philipp Arras's avatar
Philipp Arras committed
136
137
138
139
140
141
142
143
144
145

def dynamic_operator(domain,
                     harmonic_padding,
                     sm_s0,
                     sm_x0,
                     key,
                     causal=True,
                     minimum_phase=False):
    '''Constructs an operator encoding the Green's function of a linear
    homogeneous dynamic system.
Philipp Frank's avatar
Philipp Frank committed
146
147
148
149

    Parameters
    ----------
    domain : RGSpace
Philipp Arras's avatar
Philipp Arras committed
150
        The space under consideration.
Philipp Frank's avatar
Philipp Frank committed
151
    harmonic_padding : None, int, list of int
Philipp Arras's avatar
Philipp Arras committed
152
153
        Amount of central padding in harmonic space in pixels. If None the
        field is not padded at all.
Philipp Frank's avatar
Philipp Frank committed
154
    sm_s0 : float
Philipp Arras's avatar
Philipp Arras committed
155
        Cutoff for dynamic smoothness prior.
Philipp Arras's avatar
Philipp Arras committed
156
    sm_x0 : float, list of float
Philipp Arras's avatar
Philipp Arras committed
157
        Scaling of dynamic smoothness along each axis.
Philipp Frank's avatar
Philipp Frank committed
158
159
160
    key : String
        key for dynamics encoding parameter.
    causal : boolean
Philipp Arras's avatar
Philipp Arras committed
161
        Whether or not the reconstructed dynamics should be causal in time.
Philipp Frank's avatar
Philipp Frank committed
162
    minimum_phase: boolean
Philipp Arras's avatar
Philipp Arras committed
163
        Whether or not the reconstructed dynamics should be minimum phase.
Philipp Frank's avatar
Philipp Frank committed
164
165
166
167

    Returns
    -------
    Operator
Philipp Arras's avatar
Philipp Arras committed
168
        The Operator encoding the dynamic Green's function in harmonic space.
Philipp Frank's avatar
Philipp Frank committed
169
    Dictionary of Operator
Philipp Arras's avatar
Philipp Arras committed
170
        A collection of sub-chains of Operator which can be used for plotting
Philipp Arras's avatar
Philipp Arras committed
171
        and evaluation.
Philipp Frank's avatar
Philipp Frank committed
172
173
174

    Notes
    -----
Philipp Arras's avatar
Philipp Arras committed
175
    The first axis of the domain is interpreted the time axis.
Philipp Frank's avatar
Philipp Frank committed
176
    '''
Philipp Arras's avatar
Philipp Arras committed
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
    dct = {
        'domain': domain,
        'harmonic_padding': harmonic_padding,
        'sm_s0': sm_s0,
        'sm_x0': sm_x0,
        'keys': [key],
        'causal': causal,
        'cone': False,
        'minimum_phase': minimum_phase,
    }
    return _make_dynamic_operator(**dct)


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 Green's function of a linear
Philipp Arras's avatar
Philipp Arras committed
201
    homogeneous dynamic system. The Green's function is constrained to be
Philipp Arras's avatar
Philipp Arras committed
202
    within a light cone.
Philipp Frank's avatar
Philipp Frank committed
203
204
205
206
207
208

    Parameters
    ----------
    domain : RGSpace
        The space under consideration. Must have dim > 1.
    harmonic_padding : None, int, list of int
Philipp Arras's avatar
Philipp Arras committed
209
210
        Amount of central padding in harmonic space in pixels. If None the
        field is not padded at all.
Philipp Frank's avatar
Philipp Frank committed
211
    sm_s0 : float
Philipp Arras's avatar
Philipp Arras committed
212
        Cutoff for dynamic smoothness prior.
Philipp Arras's avatar
Philipp Arras committed
213
    sm_x0 : float, list of float
Philipp Arras's avatar
Philipp Arras committed
214
        Scaling of dynamic smoothness along each axis.
Philipp Frank's avatar
Philipp Frank committed
215
    key : String
Philipp Arras's avatar
Philipp Arras committed
216
        Key for dynamics encoding parameter.
Philipp Frank's avatar
Philipp Frank committed
217
    lightcone_key: String
Philipp Arras's avatar
Philipp Arras committed
218
        Key for lightspeed paramteter.
Philipp Arras's avatar
Philipp Arras committed
219
    sigc : float, list of float
Philipp Arras's avatar
Philipp Arras committed
220
        Variance of lightspeed parameter.
Philipp Frank's avatar
Philipp Frank committed
221
222
223
    quant : float
        Quantization of the light cone in pixels.
    causal : boolean
Philipp Arras's avatar
Philipp Arras committed
224
        Whether or not the reconstructed dynamics should be causal in time.
Philipp Frank's avatar
Philipp Frank committed
225
    minimum_phase: boolean
Philipp Arras's avatar
Philipp Arras committed
226
        Whether or not the reconstructed dynamics should be minimum phase.
Philipp Frank's avatar
Philipp Frank committed
227
228
229
230

    Returns
    -------
    Operator
Philipp Arras's avatar
Philipp Arras committed
231
        The Operator encoding the dynamic Green's function in harmonic space
Philipp Arras's avatar
Philipp Arras committed
232
233
        when evaluated.
    dict
Philipp Arras's avatar
Philipp Arras committed
234
        A collection of sub-chains of Operator which can be used
Philipp Arras's avatar
Philipp Arras committed
235
        for plotting and evaluation.
Philipp Frank's avatar
Philipp Frank committed
236
237
238

    Notes
    -----
Philipp Arras's avatar
Philipp Arras committed
239
    The first axis of the domain is interpreted the time axis.
Philipp Frank's avatar
Philipp Frank committed
240
    '''
Philipp Arras's avatar
Philipp Arras committed
241

Philipp Frank's avatar
Philipp Frank committed
242
243
    if len(domain.shape) < 2:
        raise ValueError("Space must be at least 2 dimensional!")
Philipp Arras's avatar
Philipp Arras committed
244
245
246
247
248
249
250
251
252
253
254
255
256
    dct = {
        'domain': domain,
        'harmonic_padding': harmonic_padding,
        'sm_s0': sm_s0,
        'sm_x0': sm_x0,
        'keys': [key, lightcone_key],
        'causal': causal,
        'cone': True,
        'minimum_phase': minimum_phase,
        'sigc': sigc,
        'quant': quant
    }
    return _make_dynamic_operator(**dct)