linearization.py 15.9 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
Martin Reinecke's avatar
Martin Reinecke committed
23
from . import utilities
Martin Reinecke's avatar
Martin Reinecke committed
24
25
26


class Linearization(object):
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
27
    """Let `A` be an operator and `x` a field. `Linearization` stores the value
28
29
30
31
32
    of the operator application (i.e. `A(x)`), the local Jacobian
    (i.e. `dA(x)/dx`) and, optionally, the local metric.

    Parameters
    ----------
Philipp Arras's avatar
Philipp Arras committed
33
34
    val : Field or MultiField
        The value of the operator application.
35
    jac : LinearOperator
Philipp Arras's avatar
Philipp Arras committed
36
37
38
39
40
41
        The Jacobian.
    metric : LinearOperator or None
        The metric. Default: None.
    want_metric : bool
        If True, the metric will be computed for other Linearizations derived
        from this one. Default: False.
42
    """
43
    def __init__(self, val, jac, metric=None, want_metric=False):
Martin Reinecke's avatar
Martin Reinecke committed
44
45
        self._val = val
        self._jac = jac
Martin Reinecke's avatar
Martin Reinecke committed
46
47
        if self._val.domain != self._jac.target:
            raise ValueError("domain mismatch")
48
        self._want_metric = want_metric
Martin Reinecke's avatar
Martin Reinecke committed
49
50
        self._metric = metric

51
    def new(self, val, jac, metric=None):
52
53
54
55
56
        """Create a new Linearization, taking the `want_metric` property from
           this one.

        Parameters
        ----------
Philipp Arras's avatar
Philipp Arras committed
57
        val : Field or MultiField
58
59
60
            the value of the operator application
        jac : LinearOperator
            the Jacobian
Philipp Arras's avatar
Philipp Arras committed
61
62
        metric : LinearOperator or None
            The metric. Default: None.
63
        """
64
65
        return Linearization(val, jac, metric, self._want_metric)

Philipp Arras's avatar
Philipp Arras committed
66
67
68
    def trivial_jac(self):
        return Linearization.make_var(self._val, self._want_metric)

Philipp Arras's avatar
Philipp Arras committed
69
70
71
72
73
74
75
    def prepend_jac(self, jac):
        metric = None
        if self._metric is not None:
            from .operators.sandwich_operator import SandwichOperator
            metric = None if self._metric is None else SandwichOperator.make(jac, self._metric)
        return self.new(self._val, self._jac @ jac, metric)

Martin Reinecke's avatar
Martin Reinecke committed
76
77
    @property
    def domain(self):
Philipp Arras's avatar
Philipp Arras committed
78
        """DomainTuple or MultiDomain : the Jacobian's domain"""
Martin Reinecke's avatar
Martin Reinecke committed
79
80
81
82
        return self._jac.domain

    @property
    def target(self):
Philipp Arras's avatar
Philipp Arras committed
83
        """DomainTuple or MultiDomain : the Jacobian's target (i.e. the value's domain)"""
Martin Reinecke's avatar
Martin Reinecke committed
84
85
86
87
        return self._jac.target

    @property
    def val(self):
Philipp Arras's avatar
Philipp Arras committed
88
        """Field or MultiField : the value"""
Martin Reinecke's avatar
Martin Reinecke committed
89
90
91
92
        return self._val

    @property
    def jac(self):
93
        """LinearOperator : the Jacobian"""
Martin Reinecke's avatar
Martin Reinecke committed
94
95
        return self._jac

Martin Reinecke's avatar
Martin Reinecke committed
96
97
    @property
    def gradient(self):
Philipp Arras's avatar
Philipp Arras committed
98
        """Field or MultiField : the gradient
99
100
101
102
103

        Notes
        -----
        Only available if target is a scalar
        """
Martin Reinecke's avatar
Martin Reinecke committed
104
        return self._jac.adjoint_times(Field.scalar(1.))
Martin Reinecke's avatar
Martin Reinecke committed
105

106
107
    @property
    def want_metric(self):
Martin Reinecke's avatar
Martin Reinecke committed
108
        """bool : True iff the metric was requested in the constructor"""
109
110
        return self._want_metric

Martin Reinecke's avatar
Martin Reinecke committed
111
112
    @property
    def metric(self):
113
114
115
116
117
118
        """LinearOperator : the metric

        Notes
        -----
        Only available if target is a scalar
        """
Martin Reinecke's avatar
Martin Reinecke committed
119
120
        return self._metric

Martin Reinecke's avatar
Martin Reinecke committed
121
    def __getitem__(self, name):
Martin Reinecke's avatar
Martin Reinecke committed
122
        return self.new(self._val[name], self._jac.ducktape_left(name))
Martin Reinecke's avatar
Martin Reinecke committed
123

Martin Reinecke's avatar
Martin Reinecke committed
124
    def __neg__(self):
125
126
        return self.new(-self._val, -self._jac,
                        None if self._metric is None else -self._metric)
Martin Reinecke's avatar
Martin Reinecke committed
127

Martin Reinecke's avatar
Martin Reinecke committed
128
    def conjugate(self):
129
        return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
130
131
132
133
134
            self._val.conjugate(), self._jac.conjugate(),
            None if self._metric is None else self._metric.conjugate())

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

Martin Reinecke's avatar
Martin Reinecke committed
137
    def _myadd(self, other, neg):
Martin Reinecke's avatar
Martin Reinecke committed
138
139
140
        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
141
                met = self._metric._myadd(other._metric, neg)
142
            return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
143
144
                self._val.flexible_addsub(other._val, neg),
                self._jac._myadd(other._jac, neg), met)
Martin Reinecke's avatar
Martin Reinecke committed
145
        if isinstance(other, (int, float, complex, Field, MultiField)):
Martin Reinecke's avatar
Martin Reinecke committed
146
            if neg:
147
                return self.new(self._val-other, self._jac, self._metric)
Martin Reinecke's avatar
Martin Reinecke committed
148
            else:
149
                return self.new(self._val+other, self._jac, self._metric)
Martin Reinecke's avatar
Martin Reinecke committed
150
151
152

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

    def __radd__(self, other):
Martin Reinecke's avatar
Martin Reinecke committed
155
        return self._myadd(other, False)
Martin Reinecke's avatar
Martin Reinecke committed
156
157

    def __sub__(self, other):
Martin Reinecke's avatar
Martin Reinecke committed
158
        return self._myadd(other, True)
Martin Reinecke's avatar
Martin Reinecke committed
159
160
161
162

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

163
164
    def __truediv__(self, other):
        if isinstance(other, Linearization):
Philipp Frank's avatar
Philipp Frank committed
165
            return self.__mul__(other.one_over())
166
167
168
        return self.__mul__(1./other)

    def __rtruediv__(self, other):
Philipp Frank's avatar
Philipp Frank committed
169
        return self.one_over().__mul__(other)
170

Martin Reinecke's avatar
Martin Reinecke committed
171
172
173
    def __pow__(self, power):
        if not np.isscalar(power):
            return NotImplemented
Martin Reinecke's avatar
Martin Reinecke committed
174
175
        return self.new(self._val**power,
                        makeOp(self._val**(power-1)).scale(power)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
176

Martin Reinecke's avatar
Martin Reinecke committed
177
178
179
    def __mul__(self, other):
        from .sugar import makeOp
        if isinstance(other, Linearization):
Martin Reinecke's avatar
Martin Reinecke committed
180
181
            if self.target != other.target:
                raise ValueError("domain mismatch")
182
            return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
183
                self._val*other._val,
Martin Reinecke's avatar
Martin Reinecke committed
184
185
                (makeOp(other._val)(self._jac))._myadd(
                 makeOp(self._val)(other._jac), False))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
186
187
188
189
        if np.isscalar(other):
            if other == 1:
                return self
            met = None if self._metric is None else self._metric.scale(other)
190
            return self.new(self._val*other, self._jac.scale(other), met)
Martin Reinecke's avatar
Martin Reinecke committed
191
        if isinstance(other, (Field, MultiField)):
Martin Reinecke's avatar
Martin Reinecke committed
192
193
            if self.target != other.domain:
                raise ValueError("domain mismatch")
194
            return self.new(self._val*other, makeOp(other)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
195
196

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

199
    def outer(self, other):
Martin Reinecke's avatar
Martin Reinecke committed
200
201
202
203
204
205
206
207
208
209
210
211
        """Computes the outer product of this Linearization with a Field or
        another Linearization

        Parameters
        ----------
        other : Field or MultiField or Linearization

        Returns
        -------
        Linearization
            the outer product of self and other
        """
212
213
214
        from .operators.outer_product_operator import OuterProduct
        if isinstance(other, Linearization):
            return self.new(
Sebastian Hutschenreuter's avatar
Sebastian Hutschenreuter committed
215
216
217
                OuterProduct(self._val, other.target)(other._val),
                OuterProduct(self._jac(self._val), other.target)._myadd(
                    OuterProduct(self._val, other.target)(other._jac), False))
218
        if np.isscalar(other):
Martin Reinecke's avatar
Martin Reinecke committed
219
            return self.__mul__(other)
220
        if isinstance(other, (Field, MultiField)):
Sebastian Hutschenreuter's avatar
Sebastian Hutschenreuter committed
221
222
            return self.new(OuterProduct(self._val, other.domain)(other),
                            OuterProduct(self._jac(self._val), other.domain))
223

Martin Reinecke's avatar
Martin Reinecke committed
224
    def vdot(self, other):
Martin Reinecke's avatar
Martin Reinecke committed
225
226
227
228
229
230
231
232
233
234
235
236
        """Computes the inner product of this Linearization with a Field or
        another Linearization

        Parameters
        ----------
        other : Field or MultiField or Linearization

        Returns
        -------
        Linearization
            the inner product of self and other
        """
Martin Reinecke's avatar
Martin Reinecke committed
237
        from .operators.simple_linear_operators import VdotOperator
Martin Reinecke's avatar
Martin Reinecke committed
238
        if isinstance(other, (Field, MultiField)):
239
            return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
240
                self._val.vdot(other),
Martin Reinecke's avatar
Martin Reinecke committed
241
                VdotOperator(other)(self._jac))
242
        return self.new(
Martin Reinecke's avatar
Martin Reinecke committed
243
            self._val.vdot(other._val),
Martin Reinecke's avatar
Martin Reinecke committed
244
245
            VdotOperator(self._val)(other._jac) +
            VdotOperator(other._val)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
246

247
    def sum(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
248
249
250
251
252
253
254
255
256
257
258
259
260
        """Computes the (partial) sum over self

        Parameters
        ----------
        spaces : None, int or list of int
            - if None, sum over the entire domain
            - else sum over the specified subspaces

        Returns
        -------
        Linearization
            the (partial) sum
        """
261
        from .operators.contraction_operator import ContractionOperator
Martin Reinecke's avatar
Martin Reinecke committed
262
263
264
        return self.new(
            self._val.sum(spaces),
            ContractionOperator(self._jac.target, spaces)(self._jac))
265
266

    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
267
268
269
270
271
272
273
274
275
276
277
278
279
        """Computes the (partial) integral over self

        Parameters
        ----------
        spaces : None, int or list of int
            - if None, integrate over the entire domain
            - else integrate over the specified subspaces

        Returns
        -------
        Linearization
            the (partial) integral
        """
280
        from .operators.contraction_operator import ContractionOperator
Martin Reinecke's avatar
Martin Reinecke committed
281
282
283
        return self.new(
            self._val.integrate(spaces),
            ContractionOperator(self._jac.target, spaces, 1)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
284
285
286

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

Martin Reinecke's avatar
Martin Reinecke committed
289
290
    def clip(self, min=None, max=None):
        tmp = self._val.clip(min, max)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
291
        if (min is None) and (max is None):
Martin Reinecke's avatar
Martin Reinecke committed
292
            return self
Jakob Knollmueller's avatar
Jakob Knollmueller committed
293
294
295
296
297
298
        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))
299
300
        return self.new(tmp, tmp2(self._jac))

Philipp Arras's avatar
Philipp Arras committed
301
    def sqrt(self):
Philipp Arras's avatar
Philipp Arras committed
302
303
        tmp = self._val.sqrt()
        return self.new(tmp, makeOp(0.5/tmp)(self._jac))
Philipp Arras's avatar
Philipp Arras committed
304

305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
    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()
Philipp Arras's avatar
Philipp Arras committed
322
        tmp2 = ((np.pi*self._val).cos()-tmp)/self._val
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
323
        ind = self._val.val == 0
Martin Reinecke's avatar
Martin Reinecke committed
324
        loc = tmp2.val_rw()
Philipp Arras's avatar
Philipp Arras committed
325
        loc[ind] = 0
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
326
        tmp2 = Field(tmp.domain, loc)
327
328
        return self.new(tmp, makeOp(tmp2)(self._jac))

Martin Reinecke's avatar
Martin Reinecke committed
329
330
    def log(self):
        tmp = self._val.log()
331
        return self.new(tmp, makeOp(1./self._val)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
332

333
334
335
336
337
    def log10(self):
        tmp = self._val.log10()
        tmp2 = 1. / (self._val * np.log(10))
        return self.new(tmp, makeOp(tmp2)(self._jac))

338
    def log1p(self):
339
340
341
        tmp = self._val.log1p()
        tmp2 = 1. / (1. + self._val)
        return self.new(tmp, makeOp(tmp2)(self.jac))
342

343
344
345
346
347
    def expm1(self):
        tmp = self._val.expm1()
        tmp2 = self._val.exp()
        return self.new(tmp, makeOp(tmp2)(self.jac))

348
349
350
351
352
353
354
355
356
357
    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
358
359
    def tanh(self):
        tmp = self._val.tanh()
360
        return self.new(tmp, makeOp(1.-tmp**2)(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
361

362
    def sigmoid(self):
Martin Reinecke's avatar
Martin Reinecke committed
363
364
        tmp = self._val.tanh()
        tmp2 = 0.5*(1.+tmp)
365
        return self.new(tmp2, makeOp(0.5*(1.-tmp**2))(self._jac))
Martin Reinecke's avatar
Martin Reinecke committed
366

367
    def absolute(self):
Martin Reinecke's avatar
Martin Reinecke committed
368
369
        if utilities.iscomplextype(self._val.dtype):
            raise TypeError("Argument must not be complex")
370
371
        tmp = self._val.absolute()
        tmp2 = self._val.sign()
Philipp Arras's avatar
Philipp Arras committed
372

Martin Reinecke's avatar
stage2    
Martin Reinecke committed
373
        ind = self._val.val == 0
Martin Reinecke's avatar
Martin Reinecke committed
374
        loc = tmp2.val_rw().astype(float)
Philipp Arras's avatar
Philipp Arras committed
375
        loc[ind] = np.nan
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
376
        tmp2 = Field(tmp.domain, loc)
Philipp Arras's avatar
Philipp Arras committed
377

378
379
380
381
382
383
384
        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
385
    def add_metric(self, metric):
386
        return self.new(self._val, self._jac, metric)
Martin Reinecke's avatar
Martin Reinecke committed
387

Martin Reinecke's avatar
Martin Reinecke committed
388
389
390
    def with_want_metric(self):
        return Linearization(self._val, self._jac, self._metric, True)

Martin Reinecke's avatar
Martin Reinecke committed
391
    @staticmethod
392
    def make_var(field, want_metric=False):
Martin Reinecke's avatar
Martin Reinecke committed
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
        """Converts a Field to a Linearization, with a unity Jacobian

        Parameters
        ----------
        field : Field or Multifield
            the field to be converted
        want_metric : bool
            If True, the metric will be computed for other Linearizations
            derived from this one. Default: False.

        Returns
        -------
        Linearization
            the requested Linearization
        """
Martin Reinecke's avatar
Martin Reinecke committed
408
        from .operators.scaling_operator import ScalingOperator
409
        return Linearization(field, ScalingOperator(field.domain, 1.),
410
                             want_metric=want_metric)
Martin Reinecke's avatar
Martin Reinecke committed
411
412

    @staticmethod
413
    def make_const(field, want_metric=False):
Martin Reinecke's avatar
Martin Reinecke committed
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
        """Converts a Field to a Linearization, with a zero Jacobian

        Parameters
        ----------
        field : Field or Multifield
            the field to be converted
        want_metric : bool
            If True, the metric will be computed for other Linearizations
            derived from this one. Default: False.

        Returns
        -------
        Linearization
            the requested Linearization

        Notes
        -----
        The Jacobian is square and contains only zeroes.
        """
Martin Reinecke's avatar
Martin Reinecke committed
433
        from .operators.simple_linear_operators import NullOperator
434
435
        return Linearization(field, NullOperator(field.domain, field.domain),
                             want_metric=want_metric)
Martin Reinecke's avatar
Martin Reinecke committed
436

Martin Reinecke's avatar
Martin Reinecke committed
437
438
    @staticmethod
    def make_const_empty_input(field, want_metric=False):
Martin Reinecke's avatar
Martin Reinecke committed
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
        """Converts a Field to a Linearization, with a zero Jacobian

        Parameters
        ----------
        field : Field or Multifield
            the field to be converted
        want_metric : bool
            If True, the metric will be computed for other Linearizations
            derived from this one. Default: False.

        Returns
        -------
        Linearization
            the requested Linearization

        Notes
        -----
        The Jacobian has an empty input domain, i.e. its matrix representation
        has 0 columns.
        """
Martin Reinecke's avatar
Martin Reinecke committed
459
460
        from .operators.simple_linear_operators import NullOperator
        from .multi_domain import MultiDomain
Martin Reinecke's avatar
Martin Reinecke committed
461
462
463
        return Linearization(
            field, NullOperator(MultiDomain.make({}), field.domain),
            want_metric=want_metric)
Martin Reinecke's avatar
Martin Reinecke committed
464

Martin Reinecke's avatar
Martin Reinecke committed
465
466
    @staticmethod
    def make_partial_var(field, constants, want_metric=False):
Martin Reinecke's avatar
Martin Reinecke committed
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
        """Converts a MultiField to a Linearization, with a Jacobian that is
        unity for some MultiField components and a zero matrix for others.

        Parameters
        ----------
        field : Multifield
            the field to be converted
        constants : list of string
            the MultiField components for which the Jacobian should be
            a zero matrix.
        want_metric : bool
            If True, the metric will be computed for other Linearizations
            derived from this one. Default: False.

        Returns
        -------
        Linearization
            the requested Linearization

        Notes
        -----
        The Jacobian is square.
        """
Martin Reinecke's avatar
Martin Reinecke committed
490
        from .operators.scaling_operator import ScalingOperator
Philipp Arras's avatar
Typo    
Philipp Arras committed
491
        from .operators.block_diagonal_operator import BlockDiagonalOperator
Martin Reinecke's avatar
Martin Reinecke committed
492
493
494
        if len(constants) == 0:
            return Linearization.make_var(field, want_metric)
        else:
495
            ops = {key: ScalingOperator(dom, 0. if key in constants else 1.)
496
497
                   for key, dom in field.domain.items()}
            bdop = BlockDiagonalOperator(field.domain, ops)
Martin Reinecke's avatar
Martin Reinecke committed
498
            return Linearization(field, bdop, want_metric=want_metric)