diagonal_operator.py 5.46 KB
 Theo Steininger committed Apr 13, 2017 1 2 3 4 5 6 7 8 9 10 11 12 ``````# 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 . `````` Theo Steininger committed May 24, 2017 13 14 15 16 17 ``````# # Copyright(C) 2013-2017 Max-Planck-Society # # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik # and financially supported by the Studienstiftung des deutschen Volkes. `````` theos committed Aug 22, 2016 18 `````` `````` Martin Reinecke committed Jun 30, 2017 19 20 ``````from __future__ import division from builtins import range `````` theos committed Aug 22, 2016 21 ``````import numpy as np `````` Martin Reinecke committed Oct 07, 2017 22 23 24 25 ``````from ..field import Field from ..domain_tuple import DomainTuple from .endomorphic_operator import EndomorphicOperator from ..nifty_utilities import cast_iseq_to_tuple `````` Martin Reinecke committed Oct 26, 2017 26 ``````from .. import dobj `````` theos committed Aug 22, 2016 27 28 `````` class DiagonalOperator(EndomorphicOperator): `````` Theo Steininger committed May 13, 2017 29 30 31 32 33 `````` """ NIFTY class for diagonal operators. The NIFTY DiagonalOperator class is a subclass derived from the EndomorphicOperator. It multiplies an input field pixel-wise with its diagonal. `````` theos committed Aug 22, 2016 34 `````` `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 35 36 `````` Parameters ---------- `````` Martin Reinecke committed Sep 24, 2017 37 `````` diagonal : Field `````` Martin Reinecke committed Oct 18, 2017 38 39 `````` The diagonal entries of the operator (already containing volume factors). `````` Martin Reinecke committed Oct 03, 2017 40 41 42 43 44 45 `````` domain : tuple of DomainObjects, i.e. Spaces and FieldTypes The domain on which the Operator's input Field lives. If None, use the domain of "diagonal". spaces : tuple of int The elements of "domain" on which the operator acts. If None, it acts on all elements. `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 46 47 48 `````` Attributes ---------- `````` Pumpe, Daniel (dpumpe) committed May 14, 2017 49 50 51 52 53 54 55 56 `````` domain : tuple of DomainObjects, i.e. Spaces and FieldTypes The domain on which the Operator's input Field lives. target : tuple of DomainObjects, i.e. Spaces and FieldTypes The domain in which the outcome of the operator lives. As the Operator is endomorphic this is the same as its domain. unitary : boolean Indicates whether the Operator is unitary or not. self_adjoint : boolean `````` Martin Reinecke committed Oct 18, 2017 57 `````` Indicates whether the operator is self-adjoint or not. `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 58 `````` `````` Martin Reinecke committed Oct 12, 2017 59 `````` NOTE: the fields given to __init__ and returned from .diagonal() are `````` Martin Reinecke committed Oct 18, 2017 60 61 `````` considered to be non-bare, i.e. during operator application, no additional volume factors are applied! `````` Martin Reinecke committed Oct 12, 2017 62 `````` `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 63 64 65 66 67 `````` See Also -------- EndomorphicOperator """ `````` Martin Reinecke committed Oct 03, 2017 68 `````` def __init__(self, diagonal, domain=None, spaces=None): `````` Martin Reinecke committed Sep 25, 2017 69 `````` super(DiagonalOperator, self).__init__() `````` Theo Steininger committed May 10, 2017 70 `````` `````` Martin Reinecke committed Sep 24, 2017 71 72 `````` if not isinstance(diagonal, Field): raise TypeError("Field object required") `````` Martin Reinecke committed Sep 25, 2017 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 `````` if domain is None: self._domain = diagonal.domain else: self._domain = DomainTuple.make(domain) if spaces is None: self._spaces = None if diagonal.domain != self._domain: raise ValueError("domain mismatch") else: self._spaces = cast_iseq_to_tuple(spaces) nspc = len(self._spaces) if nspc != len(diagonal.domain.domains): raise ValueError("spaces and domain must have the same length") if nspc > len(self._domain.domains): raise ValueError("too many spaces") if nspc > len(set(self._spaces)): raise ValueError("non-unique space indices") # if nspc==len(self.diagonal.domain.domains, we could do some optimization for i, j in enumerate(self._spaces): if diagonal.domain[i] != self._domain[j]: raise ValueError("domain mismatch") `````` Martin Reinecke committed Oct 18, 2017 95 `````` self._diagonal = diagonal.copy() `````` Theo Steininger committed Aug 26, 2017 96 97 `````` self._self_adjoint = None self._unitary = None `````` theos committed Aug 22, 2016 98 `````` `````` Martin Reinecke committed Sep 25, 2017 99 100 `````` def _times(self, x): return self._times_helper(x, lambda z: z.__mul__) `````` theos committed Aug 22, 2016 101 `````` `````` Martin Reinecke committed Sep 25, 2017 102 103 `````` def _adjoint_times(self, x): return self._times_helper(x, lambda z: z.conjugate().__mul__) `````` theos committed Aug 22, 2016 104 `````` `````` Martin Reinecke committed Sep 25, 2017 105 106 `````` def _inverse_times(self, x): return self._times_helper(x, lambda z: z.__rtruediv__) `````` theos committed Aug 22, 2016 107 `````` `````` Martin Reinecke committed Sep 25, 2017 108 109 `````` def _adjoint_inverse_times(self, x): return self._times_helper(x, lambda z: z.conjugate().__rtruediv__) `````` theos committed Aug 22, 2016 110 `````` `````` Martin Reinecke committed Oct 03, 2017 111 `````` def diagonal(self): `````` Theo Steininger committed Jun 21, 2017 112 113 114 115 116 117 118 `````` """ Returns the diagonal of the Operator. Returns ------- out : Field The diagonal of the Operator. """ `````` Martin Reinecke committed Oct 18, 2017 119 `````` return self._diagonal.copy() `````` Theo Steininger committed Jun 21, 2017 120 `````` `````` theos committed Sep 17, 2016 121 122 `````` @property def domain(self): `````` Martin Reinecke committed Sep 25, 2017 123 `````` return self._domain `````` theos committed Sep 17, 2016 124 `````` `````` theos committed Aug 22, 2016 125 `````` @property `````` Martin Reinecke committed Apr 28, 2017 126 127 `````` def self_adjoint(self): if self._self_adjoint is None: `````` Martin Reinecke committed Sep 16, 2017 128 129 130 131 `````` if not issubclass(self._diagonal.dtype.type, np.complexfloating): self._self_adjoint = True else: self._self_adjoint = (self._diagonal.val.imag == 0).all() `````` Martin Reinecke committed Apr 28, 2017 132 `````` return self._self_adjoint `````` theos committed Aug 22, 2016 133 134 135 `````` @property def unitary(self): `````` Theo Steininger committed Feb 16, 2017 136 `````` if self._unitary is None: `````` Martin Reinecke committed Sep 16, 2017 137 `````` self._unitary = (abs(self._diagonal.val) == 1.).all() `````` theos committed Aug 22, 2016 138 139 `````` return self._unitary `````` Martin Reinecke committed Sep 25, 2017 140 141 `````` def _times_helper(self, x, operation): if self._spaces is None: `````` Martin Reinecke committed Sep 07, 2017 142 `````` return operation(self._diagonal)(x) `````` theos committed Sep 13, 2016 143 `````` `````` Martin Reinecke committed Sep 25, 2017 144 145 146 `````` active_axes = [] for space_index in self._spaces: active_axes += x.domain.axes[space_index] `````` theos committed Sep 13, 2016 147 `````` `````` Martin Reinecke committed Oct 26, 2017 148 149 `````` reshaper = [shp if i in active_axes else 1 for i, shp in enumerate(x.shape)] `````` Martin Reinecke committed Oct 26, 2017 150 `````` reshaped_local_diagonal = self._diagonal.val.reshape(reshaper) `````` theos committed Sep 13, 2016 151 152 `````` # here the actual multiplication takes place `````` Martin Reinecke committed Sep 16, 2017 153 `` return Field(x.domain, val=operation(reshaped_local_diagonal)(x.val))``