linearization.py 10.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 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/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Martin Reinecke's avatar
Martin Reinecke committed
17
18
19
20

import numpy as np

from .field import Field
Martin Reinecke's avatar
Martin Reinecke committed
21
from .multi_field import MultiField
Martin Reinecke's avatar
Martin Reinecke committed
22
from .sugar import makeOp
Jakob Knollmueller's avatar
Jakob Knollmueller committed
23
from .operators.scaling_operator import ScalingOperator
Martin Reinecke's avatar
Martin Reinecke committed
24
25
26


class Linearization(object):
27
    def __init__(self, val, jac, metric=None, want_metric=False):
Martin Reinecke's avatar
Martin Reinecke committed
28
29
        self._val = val
        self._jac = jac
Martin Reinecke's avatar
Martin Reinecke committed
30
31
        if self._val.domain != self._jac.target:
            raise ValueError("domain mismatch")
32
        self._want_metric = want_metric
Martin Reinecke's avatar
Martin Reinecke committed
33
34
        self._metric = metric

35
36
37
    def new(self, val, jac, metric=None):
        return Linearization(val, jac, metric, self._want_metric)

Martin Reinecke's avatar
Martin Reinecke committed
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
    @property
    def domain(self):
        return self._jac.domain

    @property
    def target(self):
        return self._jac.target

    @property
    def val(self):
        return self._val

    @property
    def jac(self):
        return self._jac

Martin Reinecke's avatar
Martin Reinecke committed
54
55
56
    @property
    def gradient(self):
        """Only available if target is a scalar"""
Martin Reinecke's avatar
Martin Reinecke committed
57
        return self._jac.adjoint_times(Field.scalar(1.))
Martin Reinecke's avatar
Martin Reinecke committed
58

59
60
61
62
    @property
    def want_metric(self):
        return self._want_metric

Martin Reinecke's avatar
Martin Reinecke committed
63
64
    @property
    def metric(self):
Martin Reinecke's avatar
Martin Reinecke committed
65
        """Only available if target is a scalar"""
Martin Reinecke's avatar
Martin Reinecke committed
66
67
        return self._metric

Martin Reinecke's avatar
Martin Reinecke committed
68
    def __getitem__(self, name):
Martin Reinecke's avatar
Martin Reinecke committed
69
70
        from .operators.simple_linear_operators import ducktape
        return self.new(self._val[name], ducktape(None, self.domain, name))
Martin Reinecke's avatar
Martin Reinecke committed
71

Martin Reinecke's avatar
Martin Reinecke committed
72
    def __neg__(self):
73
74
        return self.new(-self._val, -self._jac,
                        None if self._metric is None else -self._metric)
Martin Reinecke's avatar
Martin Reinecke committed
75

Martin Reinecke's avatar
Martin Reinecke committed
76
    def conjugate(self):
77
        return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
78
79
80
81
82
            self._val.conjugate(), self._jac.conjugate(),
            None if self._metric is None else self._metric.conjugate())

    @property
    def real(self):
83
        return self.new(self._val.real, self._jac.real)
Martin Reinecke's avatar
Martin Reinecke committed
84

Martin Reinecke's avatar
Martin Reinecke committed
85
    def _myadd(self, other, neg):
Martin Reinecke's avatar
Martin Reinecke committed
86
87
88
        if isinstance(other, Linearization):
            met = None
            if self._metric is not None and other._metric is not None:
Martin Reinecke's avatar
Martin Reinecke committed
89
                met = self._metric._myadd(other._metric, neg)
90
            return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
91
92
                self._val.flexible_addsub(other._val, neg),
                self._jac._myadd(other._jac, neg), met)
Martin Reinecke's avatar
Martin Reinecke committed
93
        if isinstance(other, (int, float, complex, Field, MultiField)):
Martin Reinecke's avatar
Martin Reinecke committed
94
            if neg:
95
                return self.new(self._val-other, self._jac, self._metric)
Martin Reinecke's avatar
Martin Reinecke committed
96
            else:
97
                return self.new(self._val+other, self._jac, self._metric)
Martin Reinecke's avatar
Martin Reinecke committed
98
99
100

    def __add__(self, other):
        return self._myadd(other, False)
Martin Reinecke's avatar
Martin Reinecke committed
101
102

    def __radd__(self, other):
Martin Reinecke's avatar
Martin Reinecke committed
103
        return self._myadd(other, False)
Martin Reinecke's avatar
Martin Reinecke committed
104
105

    def __sub__(self, other):
Martin Reinecke's avatar
Martin Reinecke committed
106
        return self._myadd(other, True)
Martin Reinecke's avatar
Martin Reinecke committed
107
108
109
110

    def __rsub__(self, other):
        return (-self).__add__(other)

111
112
113
114
115
116
117
118
    def __truediv__(self, other):
        if isinstance(other, Linearization):
            return self.__mul__(other.inverse())
        return self.__mul__(1./other)

    def __rtruediv__(self, other):
        return self.inverse().__mul__(other)

Martin Reinecke's avatar
Martin Reinecke committed
119
120
121
    def __pow__(self, power):
        if not np.isscalar(power):
            return NotImplemented
Martin Reinecke's avatar
Martin Reinecke committed
122
123
        return self.new(self._val**power,
                        makeOp(self._val**(power-1)).scale(power)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
124

125
126
127
    def inverse(self):
        return self.new(1./self._val, makeOp(-1./(self._val**2))(self._jac))

Martin Reinecke's avatar
Martin Reinecke committed
128
129
130
    def __mul__(self, other):
        from .sugar import makeOp
        if isinstance(other, Linearization):
Martin Reinecke's avatar
Martin Reinecke committed
131
132
            if self.target != other.target:
                raise ValueError("domain mismatch")
133
            return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
134
                self._val*other._val,
Martin Reinecke's avatar
Martin Reinecke committed
135
136
                (makeOp(other._val)(self._jac))._myadd(
                 makeOp(self._val)(other._jac), False))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
137
138
139
140
        if np.isscalar(other):
            if other == 1:
                return self
            met = None if self._metric is None else self._metric.scale(other)
141
            return self.new(self._val*other, self._jac.scale(other), met)
Martin Reinecke's avatar
Martin Reinecke committed
142
        if isinstance(other, (Field, MultiField)):
Martin Reinecke's avatar
Martin Reinecke committed
143
144
            if self.target != other.domain:
                raise ValueError("domain mismatch")
145
            return self.new(self._val*other, makeOp(other)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
146
147

    def __rmul__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
148
        return self.__mul__(other)
Martin Reinecke's avatar
Martin Reinecke committed
149

150
151
152
153
    def outer(self, other):
        from .operators.outer_product_operator import OuterProduct
        if isinstance(other, Linearization):
            return self.new(
Sebastian Hutschenreuter's avatar
Sebastian Hutschenreuter committed
154
155
156
                OuterProduct(self._val, other.target)(other._val),
                OuterProduct(self._jac(self._val), other.target)._myadd(
                    OuterProduct(self._val, other.target)(other._jac), False))
157
        if np.isscalar(other):
Martin Reinecke's avatar
Martin Reinecke committed
158
            return self.__mul__(other)
159
        if isinstance(other, (Field, MultiField)):
Sebastian Hutschenreuter's avatar
Sebastian Hutschenreuter committed
160
161
            return self.new(OuterProduct(self._val, other.domain)(other),
                            OuterProduct(self._jac(self._val), other.domain))
162

Martin Reinecke's avatar
Martin Reinecke committed
163
    def vdot(self, other):
Martin Reinecke's avatar
Martin Reinecke committed
164
        from .operators.simple_linear_operators import VdotOperator
Martin Reinecke's avatar
Martin Reinecke committed
165
        if isinstance(other, (Field, MultiField)):
166
            return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
167
                Field.scalar(self._val.vdot(other)),
Martin Reinecke's avatar
Martin Reinecke committed
168
                VdotOperator(other)(self._jac))
169
        return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
170
            Field.scalar(self._val.vdot(other._val)),
Martin Reinecke's avatar
Martin Reinecke committed
171
172
            VdotOperator(self._val)(other._jac) +
            VdotOperator(other._val)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
173

174
    def sum(self, spaces=None):
175
        from .operators.contraction_operator import ContractionOperator
176
177
178
        if spaces is None:
            return self.new(
                Field.scalar(self._val.sum()),
179
                ContractionOperator(self._jac.target, None)(self._jac))
180
181
182
        else:
            return self.new(
                self._val.sum(spaces),
183
                ContractionOperator(self._jac.target, spaces)(self._jac))
184
185

    def integrate(self, spaces=None):
186
        from .operators.contraction_operator import ContractionOperator
187
188
189
        if spaces is None:
            return self.new(
                Field.scalar(self._val.integrate()),
190
                ContractionOperator(self._jac.target, None, 1)(self._jac))
191
192
193
        else:
            return self.new(
                self._val.integrate(spaces),
194
                ContractionOperator(self._jac.target, spaces, 1)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
195
196
197

    def exp(self):
        tmp = self._val.exp()
198
        return self.new(tmp, makeOp(tmp)(self._jac))
Philipp Arras's avatar
Philipp Arras committed
199

Martin Reinecke's avatar
Martin Reinecke committed
200
201
    def clip(self, min=None, max=None):
        tmp = self._val.clip(min, max)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
202
        if (min is None) and (max is None):
Martin Reinecke's avatar
Martin Reinecke committed
203
            return self
Jakob Knollmueller's avatar
Jakob Knollmueller committed
204
205
206
207
208
209
        elif max is None:
            tmp2 = makeOp(1. - (tmp == min))
        elif min is None:
            tmp2 = makeOp(1. - (tmp == max))
        else:
            tmp2 = makeOp(1. - (tmp == min) - (tmp == max))
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
        return self.new(tmp, tmp2(self._jac))

    def sin(self):
        tmp = self._val.sin()
        tmp2 = self._val.cos()
        return self.new(tmp, makeOp(tmp2)(self._jac))

    def cos(self):
        tmp = self._val.cos()
        tmp2 = - self._val.sin()
        return self.new(tmp, makeOp(tmp2)(self._jac))

    def tan(self):
        tmp = self._val.tan()
        tmp2 = 1./(self._val.cos()**2)
        return self.new(tmp, makeOp(tmp2)(self._jac))

    def sinc(self):
        tmp = self._val.sinc()
        tmp2 = (self._val.cos()-tmp)/self._val
        return self.new(tmp, makeOp(tmp2)(self._jac))

Martin Reinecke's avatar
Martin Reinecke committed
232
233
    def log(self):
        tmp = self._val.log()
234
        return self.new(tmp, makeOp(1./self._val)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
235

236
237
238
239
240
241
242
243
244
245
    def sinh(self):
        tmp = self._val.sinh()
        tmp2 = self._val.cosh()
        return self.new(tmp, makeOp(tmp2)(self._jac))

    def cosh(self):
        tmp = self._val.cosh()
        tmp2 = self._val.sinh()
        return self.new(tmp, makeOp(tmp2)(self._jac))

Martin Reinecke's avatar
Martin Reinecke committed
246
247
    def tanh(self):
        tmp = self._val.tanh()
248
        return self.new(tmp, makeOp(1.-tmp**2)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
249

250
    def sigmoid(self):
Martin Reinecke's avatar
Martin Reinecke committed
251
252
        tmp = self._val.tanh()
        tmp2 = 0.5*(1.+tmp)
253
        return self.new(tmp2, makeOp(0.5*(1.-tmp**2))(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
254

255
256
257
258
259
260
261
262
263
264
    def absolute(self):
        tmp = self._val.absolute()
        tmp2 = self._val.sign()
        return self.new(tmp, makeOp(tmp2)(self._jac))

    def one_over(self):
        tmp = 1./self._val
        tmp2 = - tmp/self._val
        return self.new(tmp, makeOp(tmp2)(self._jac))

Martin Reinecke's avatar
Martin Reinecke committed
265
    def add_metric(self, metric):
266
        return self.new(self._val, self._jac, metric)
Martin Reinecke's avatar
Martin Reinecke committed
267

Martin Reinecke's avatar
Martin Reinecke committed
268
269
270
    def with_want_metric(self):
        return Linearization(self._val, self._jac, self._metric, True)

Martin Reinecke's avatar
Martin Reinecke committed
271
    @staticmethod
272
    def make_var(field, want_metric=False):
Martin Reinecke's avatar
Martin Reinecke committed
273
        from .operators.scaling_operator import ScalingOperator
274
275
        return Linearization(field, ScalingOperator(1., field.domain),
                             want_metric=want_metric)
Martin Reinecke's avatar
Martin Reinecke committed
276
277

    @staticmethod
278
    def make_const(field, want_metric=False):
Martin Reinecke's avatar
Martin Reinecke committed
279
        from .operators.simple_linear_operators import NullOperator
280
281
        return Linearization(field, NullOperator(field.domain, field.domain),
                             want_metric=want_metric)
Martin Reinecke's avatar
Martin Reinecke committed
282

Martin Reinecke's avatar
Martin Reinecke committed
283
284
285
286
    @staticmethod
    def make_const_empty_input(field, want_metric=False):
        from .operators.simple_linear_operators import NullOperator
        from .multi_domain import MultiDomain
Martin Reinecke's avatar
Martin Reinecke committed
287
288
289
        return Linearization(
            field, NullOperator(MultiDomain.make({}), field.domain),
            want_metric=want_metric)
Martin Reinecke's avatar
Martin Reinecke committed
290

Martin Reinecke's avatar
Martin Reinecke committed
291
292
293
    @staticmethod
    def make_partial_var(field, constants, want_metric=False):
        from .operators.scaling_operator import ScalingOperator
Philipp Arras's avatar
Typo    
Philipp Arras committed
294
        from .operators.block_diagonal_operator import BlockDiagonalOperator
Martin Reinecke's avatar
Martin Reinecke committed
295
296
297
298
299
        if len(constants) == 0:
            return Linearization.make_var(field, want_metric)
        else:
            ops = [ScalingOperator(0. if key in constants else 1., dom)
                   for key, dom in field.domain.items()]
Philipp Arras's avatar
Typo    
Philipp Arras committed
300
            bdop = BlockDiagonalOperator(field.domain, tuple(ops))
Martin Reinecke's avatar
Martin Reinecke committed
301
            return Linearization(field, bdop, want_metric=want_metric)