diagonal_operator.py 5.67 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 ``````from __future__ import division `````` theos committed Aug 22, 2016 20 ``````import numpy as np `````` Martin Reinecke committed Oct 07, 2017 21 22 23 24 ``````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 Nov 09, 2017 25 ``````from .. import dobj `````` Martin Reinecke committed Oct 28, 2017 26 `````` `````` 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 `````` 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") `````` Martin Reinecke committed Oct 31, 2017 90 `````` # if nspc==len(self.diagonal.domain), `````` Martin Reinecke committed Oct 28, 2017 91 92 `````` # we could do some optimization for i, j in enumerate(self._spaces): `````` Martin Reinecke committed Sep 25, 2017 93 94 `````` if diagonal.domain[i] != self._domain[j]: raise ValueError("domain mismatch") `````` Martin Reinecke committed Nov 06, 2017 95 96 97 98 99 100 101 102 103 104 `````` if self._spaces == tuple(range(len(self._domain.domains))): self._spaces = None # shortcut if self._spaces is not None: active_axes = [] for space_index in self._spaces: active_axes += self._domain.axes[space_index] self._reshaper = [shp if i in active_axes else 1 for i, shp in enumerate(self._domain.shape)] `````` Martin Reinecke committed Sep 25, 2017 105 `````` `````` Martin Reinecke committed Oct 18, 2017 106 `````` self._diagonal = diagonal.copy() `````` Theo Steininger committed Aug 26, 2017 107 108 `````` self._self_adjoint = None self._unitary = None `````` theos committed Aug 22, 2016 109 `````` `````` Martin Reinecke committed Sep 25, 2017 110 `````` def _times(self, x): `````` Martin Reinecke committed Oct 29, 2017 111 `````` return self._times_helper(x, self._diagonal) `````` theos committed Aug 22, 2016 112 `````` `````` Martin Reinecke committed Sep 25, 2017 113 `````` def _adjoint_times(self, x): `````` Martin Reinecke committed Oct 29, 2017 114 `````` return self._times_helper(x, self._diagonal.conj()) `````` theos committed Aug 22, 2016 115 `````` `````` Martin Reinecke committed Sep 25, 2017 116 `````` def _inverse_times(self, x): `````` Martin Reinecke committed Oct 29, 2017 117 `````` return self._times_helper(x, 1./self._diagonal) `````` theos committed Aug 22, 2016 118 `````` `````` Martin Reinecke committed Sep 25, 2017 119 `````` def _adjoint_inverse_times(self, x): `````` Martin Reinecke committed Oct 29, 2017 120 `````` return self._times_helper(x, 1./self._diagonal.conj()) `````` theos committed Aug 22, 2016 121 `````` `````` Martin Reinecke committed Oct 03, 2017 122 `````` def diagonal(self): `````` Theo Steininger committed Jun 21, 2017 123 124 125 126 127 128 129 `````` """ Returns the diagonal of the Operator. Returns ------- out : Field The diagonal of the Operator. """ `````` Martin Reinecke committed Oct 18, 2017 130 `````` return self._diagonal.copy() `````` Theo Steininger committed Jun 21, 2017 131 `````` `````` theos committed Sep 17, 2016 132 133 `````` @property def domain(self): `````` Martin Reinecke committed Sep 25, 2017 134 `````` return self._domain `````` theos committed Sep 17, 2016 135 `````` `````` theos committed Aug 22, 2016 136 `````` @property `````` Martin Reinecke committed Apr 28, 2017 137 138 `````` def self_adjoint(self): if self._self_adjoint is None: `````` Martin Reinecke committed Sep 16, 2017 139 140 141 142 `````` 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 143 `````` return self._self_adjoint `````` theos committed Aug 22, 2016 144 145 146 `````` @property def unitary(self): `````` Theo Steininger committed Feb 16, 2017 147 `````` if self._unitary is None: `````` Martin Reinecke committed Sep 16, 2017 148 `````` self._unitary = (abs(self._diagonal.val) == 1.).all() `````` theos committed Aug 22, 2016 149 150 `````` return self._unitary `````` Martin Reinecke committed Oct 29, 2017 151 `````` def _times_helper(self, x, diag): `````` Martin Reinecke committed Sep 25, 2017 152 `````` if self._spaces is None: `````` Martin Reinecke committed Oct 29, 2017 153 `````` return diag*x `````` theos committed Sep 13, 2016 154 `````` `````` Martin Reinecke committed Nov 09, 2017 155 156 `````` reshaped_local_diagonal = np.reshape(dobj.to_global_data(diag.val), self._reshaper) if 0 in self._spaces: `````` Martin Reinecke committed Nov 12, 2017 157 `````` reshaped_local_diagonal = dobj.local_data(dobj.from_global_data(reshaped_local_diagonal)) `````` Martin Reinecke committed Oct 29, 2017 158 `` return Field(x.domain, val=x.val*reshaped_local_diagonal)``