diagonal_operator.py 10.8 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 19 20 21 22 23 `````` import numpy as np from d2o import distributed_data_object,\ STRATEGIES as DISTRIBUTION_STRATEGIES `````` Theo Steininger committed May 11, 2017 24 ``````from nifty.basic_arithmetics import log as nifty_log `````` theos committed Oct 06, 2016 25 ``````from nifty.config import nifty_configuration as gc `````` theos committed Aug 22, 2016 26 27 28 29 30 ``````from nifty.field import Field from nifty.operators.endomorphic_operator import EndomorphicOperator class DiagonalOperator(EndomorphicOperator): `````` Theo Steininger committed May 13, 2017 31 32 33 34 35 `````` """ 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 36 `````` `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 37 38 39 `````` Parameters ---------- `````` Theo Steininger committed May 13, 2017 40 41 42 `````` domain : tuple of DomainObjects, i.e. Spaces and FieldTypes The domain on which the Operator's input Field lives. diagonal : {scalar, list, array, Field, d2o-object} `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 43 44 `````` The diagonal entries of the operator. bare : boolean `````` Theo Steininger committed May 13, 2017 45 46 `````` Indicates whether the input for the diagonal is bare or not (default: False). `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 47 48 49 50 51 `````` copy : boolean Internal copy of the diagonal (default: True) distribution_strategy : string setting the prober distribution_strategy of the diagonal (default : None). In case diagonal is d2o-object or Field, `````` Theo Steininger committed May 13, 2017 52 `````` their distribution_strategy is used as a fallback. `````` Pumpe, Daniel (dpumpe) committed May 14, 2017 53 54 55 `````` default_spaces : tuple of ints *optional* Defines on which space(s) of a given field the Operator acts by default (default: None) `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 56 57 58 `````` Attributes ---------- `````` Pumpe, Daniel (dpumpe) committed May 14, 2017 59 60 61 62 63 64 65 66 67 `````` 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 Indicates whether the operator is self_adjoint or not. `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 68 `````` distribution_strategy : string `````` Theo Steininger committed May 13, 2017 69 70 `````` Defines the distribution_strategy of the distributed_data_object in which the diagonal entries are stored in. `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 `````` Raises ------ Notes ----- The ambiguity of bare or non-bare diagonal entries is based on the choice of a matrix representation of the operator in question. The naive choice of absorbing the volume weights into the matrix leads to a matrix-vector calculus with the non-bare entries which seems intuitive, though. The choice of keeping matrix entries and volume weights separate deals with the bare entries that allow for correct interpretation of the matrix entries; e.g., as variance in case of an covariance operator. Examples -------- >>> x_space = RGSpace(5) `````` Theo Steininger committed May 13, 2017 88 89 `````` >>> D = DiagonalOperator(x_space, diagonal=[1., 3., 2., 4., 6.]) >>> f = Field(x_space, val=2.) `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 90 91 92 `````` >>> res = D.times(f) >>> res.val `````` Theo Steininger committed May 13, 2017 93 `````` array([ 2., 6., 4., 8., 12.]) `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 94 95 96 97 98 99 100 `````` See Also -------- EndomorphicOperator """ `````` theos committed Aug 22, 2016 101 102 `````` # ---Overwritten properties and methods--- `````` Theo Steininger committed May 10, 2017 103 104 105 106 `````` def __init__(self, domain=(), diagonal=None, bare=False, copy=True, distribution_strategy=None, default_spaces=None): super(DiagonalOperator, self).__init__(default_spaces) `````` theos committed Sep 17, 2016 107 `````` self._domain = self._parse_domain(domain) `````` theos committed Aug 22, 2016 108 `````` `````` theos committed Sep 02, 2016 109 `````` if distribution_strategy is None: `````` theos committed Aug 22, 2016 110 `````` if isinstance(diagonal, distributed_data_object): `````` theos committed Sep 02, 2016 111 `````` distribution_strategy = diagonal.distribution_strategy `````` theos committed Aug 22, 2016 112 `````` elif isinstance(diagonal, Field): `````` theos committed Sep 02, 2016 113 `````` distribution_strategy = diagonal.distribution_strategy `````` theos committed Aug 22, 2016 114 `````` `````` theos committed Sep 13, 2016 115 `````` self._distribution_strategy = self._parse_distribution_strategy( `````` theos committed Sep 02, 2016 116 117 `````` distribution_strategy=distribution_strategy, val=diagonal) `````` theos committed Aug 22, 2016 118 119 120 `````` self.set_diagonal(diagonal=diagonal, bare=bare, copy=copy) `````` Theo Steininger committed Feb 07, 2017 121 122 `````` def _times(self, x, spaces): return self._times_helper(x, spaces, operation=lambda z: z.__mul__) `````` theos committed Aug 22, 2016 123 `````` `````` Theo Steininger committed Feb 07, 2017 124 125 `````` def _adjoint_times(self, x, spaces): return self._times_helper(x, spaces, `````` theos committed Sep 13, 2016 126 `````` operation=lambda z: z.adjoint().__mul__) `````` theos committed Aug 22, 2016 127 `````` `````` Theo Steininger committed Feb 07, 2017 128 129 `````` def _inverse_times(self, x, spaces): return self._times_helper(x, spaces, operation=lambda z: z.__rdiv__) `````` theos committed Aug 22, 2016 130 `````` `````` Theo Steininger committed Feb 07, 2017 131 132 `````` def _adjoint_inverse_times(self, x, spaces): return self._times_helper(x, spaces, `````` theos committed Sep 13, 2016 133 `````` operation=lambda z: z.adjoint().__rdiv__) `````` theos committed Aug 22, 2016 134 `````` `````` Theo Steininger committed Jun 21, 2017 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 `````` def diagonal(self, bare=False, copy=True): """ Returns the diagonal of the Operator. Parameters ---------- bare : boolean Whether the returned Field values should be bare or not. copy : boolean Whether the returned Field should be copied or not. Returns ------- out : Field The diagonal of the Operator. """ if bare: diagonal = self._diagonal.weight(power=-1) elif copy: diagonal = self._diagonal.copy() else: diagonal = self._diagonal return diagonal def inverse_diagonal(self, bare=False): """ Returns the inverse-diagonal of the operator. Parameters ---------- bare : boolean Whether the returned Field values should be bare or not. Returns ------- out : Field The inverse of the diagonal of the Operator. """ return 1./self.diagonal(bare=bare, copy=False) `````` theos committed Aug 22, 2016 175 176 `````` # ---Mandatory properties and methods--- `````` theos committed Sep 17, 2016 177 178 179 180 `````` @property def domain(self): return self._domain `````` theos committed Aug 22, 2016 181 `````` @property `````` Martin Reinecke committed Apr 28, 2017 182 183 184 185 `````` def self_adjoint(self): if self._self_adjoint is None: self._self_adjoint = (self._diagonal.val.imag == 0).all() return self._self_adjoint `````` theos committed Aug 22, 2016 186 187 188 `````` @property def unitary(self): `````` Theo Steininger committed Feb 16, 2017 189 190 191 `````` if self._unitary is None: self._unitary = (self._diagonal.val * self._diagonal.val.conjugate() == 1).all() `````` theos committed Aug 22, 2016 192 193 194 195 196 `````` return self._unitary # ---Added properties and methods--- @property `````` theos committed Sep 02, 2016 197 `````` def distribution_strategy(self): `````` Pumpe, Daniel (dpumpe) committed May 09, 2017 198 199 200 `````` """ distribution_strategy : string Defines the way how the diagonal operator is distributed `````` Theo Steininger committed May 13, 2017 201 202 `````` among the nodes. Available distribution_strategies are: 'fftw', 'equal' and 'not'. `````` Pumpe, Daniel (dpumpe) committed May 09, 2017 203 204 205 `````` Notes : https://arxiv.org/abs/1606.05385 `````` Theo Steininger committed May 13, 2017 206 `````` `````` Pumpe, Daniel (dpumpe) committed May 09, 2017 207 `````` """ `````` Theo Steininger committed May 13, 2017 208 `````` `````` theos committed Sep 02, 2016 209 `````` return self._distribution_strategy `````` theos committed Aug 22, 2016 210 `````` `````` theos committed Sep 02, 2016 211 212 `````` def _parse_distribution_strategy(self, distribution_strategy, val): if distribution_strategy is None: `````` theos committed Aug 22, 2016 213 `````` if isinstance(val, distributed_data_object): `````` theos committed Sep 02, 2016 214 `````` distribution_strategy = val.distribution_strategy `````` theos committed Aug 22, 2016 215 `````` elif isinstance(val, Field): `````` theos committed Sep 02, 2016 216 `````` distribution_strategy = val.distribution_strategy `````` theos committed Aug 22, 2016 217 `````` else: `````` theos committed Oct 18, 2016 218 `````` self.logger.info("Datamodel set to default!") `````` theos committed Sep 02, 2016 219 220 `````` distribution_strategy = gc['default_distribution_strategy'] elif distribution_strategy not in DISTRIBUTION_STRATEGIES['all']: `````` theos committed Oct 06, 2016 221 222 `````` raise ValueError( "Invalid distribution_strategy!") `````` theos committed Sep 02, 2016 223 `````` return distribution_strategy `````` theos committed Aug 22, 2016 224 225 `````` def set_diagonal(self, diagonal, bare=False, copy=True): `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 226 227 228 229 `````` """ Sets the diagonal of the Operator. Parameters ---------- `````` Theo Steininger committed May 13, 2017 230 `````` diagonal : {scalar, list, array, Field, d2o-object} `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 231 232 `````` The diagonal entries of the operator. bare : boolean `````` Theo Steininger committed May 13, 2017 233 234 `````` Indicates whether the input for the diagonal is bare or not (default: False). `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 235 `````` copy : boolean `````` Theo Steininger committed May 13, 2017 236 `````` Specifies if a copy of the input shall be made (default: True). `````` Pumpe, Daniel (dpumpe) committed May 08, 2017 237 238 239 `````` """ `````` theos committed Aug 22, 2016 240 241 242 `````` # use the casting functionality from Field to process `diagonal` f = Field(domain=self.domain, val=diagonal, `````` theos committed Sep 02, 2016 243 `````` distribution_strategy=self.distribution_strategy, `````` theos committed Aug 22, 2016 244 245 `````` copy=copy) `````` Pumpe, Daniel (dpumpe) committed May 09, 2017 246 `````` # weight if the given values were `bare` is True `````` theos committed Aug 24, 2016 247 `````` # do inverse weightening if the other way around `````` Pumpe, Daniel (dpumpe) committed May 09, 2017 248 `````` if bare: `````` theos committed Aug 24, 2016 249 250 251 `````` # If `copy` is True, we won't change external data by weightening # Otherwise, inplace weightening would change the external field f.weight(inplace=copy) `````` theos committed Aug 22, 2016 252 `````` `````` Martin Reinecke committed Apr 28, 2017 253 254 `````` # Reset the self_adjoint property: self._self_adjoint = None `````` theos committed Aug 22, 2016 255 `````` `````` Theo Steininger committed Feb 16, 2017 256 257 `````` # Reset the unitarity property self._unitary = None `````` theos committed Aug 22, 2016 258 259 260 `````` # store the diagonal-field self._diagonal = f `````` theos committed Sep 13, 2016 261 `````` `````` Theo Steininger committed Feb 07, 2017 262 263 `````` def _times_helper(self, x, spaces, operation): # if the domain matches directly `````` theos committed Sep 13, 2016 264 `````` # -> multiply the fields directly `````` Theo Steininger committed Feb 07, 2017 265 `````` if x.domain == self.domain: `````` theos committed Sep 13, 2016 266 267 268 269 270 271 272 `````` # here the actual multiplication takes place return operation(self.diagonal(copy=False))(x) # if the distribution_strategy of self is sub-slice compatible to # the one of x, reshape the local data of self and apply it directly active_axes = [] if spaces is None: `````` Theo Steininger committed Feb 15, 2017 273 `````` active_axes = range(len(x.shape)) `````` theos committed Sep 13, 2016 274 275 276 277 278 279 280 281 282 283 `````` else: for space_index in spaces: active_axes += x.domain_axes[space_index] axes_local_distribution_strategy = \ x.val.get_axes_local_distribution_strategy(active_axes) if axes_local_distribution_strategy == self.distribution_strategy: local_diagonal = self._diagonal.val.get_local_data(copy=False) else: # create an array that is sub-slice compatible `````` Theo Steininger committed Feb 15, 2017 284 285 `````` self.logger.warn("The input field is not sub-slice compatible to " "the distribution strategy of the operator.") `````` theos committed Sep 13, 2016 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 `````` redistr_diagonal_val = self._diagonal.val.copy( distribution_strategy=axes_local_distribution_strategy) local_diagonal = redistr_diagonal_val.get_local_data(copy=False) reshaper = [x.shape[i] if i in active_axes else 1 for i in xrange(len(x.shape))] reshaped_local_diagonal = np.reshape(local_diagonal, reshaper) # here the actual multiplication takes place local_result = operation(reshaped_local_diagonal)( x.val.get_local_data(copy=False)) result_field = x.copy_empty(dtype=local_result.dtype) result_field.val.set_local_data(local_result, copy=False) return result_field``````