propagator_operator.py 5.12 KB
 Theo Steininger committed Apr 13, 2017 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 # NIFTy # Copyright (C) 2017 Theo Steininger # # Author: Theo Steininger # # 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 Jan 27, 2017 18   Theo Steininger committed Oct 14, 2016 19 from nifty.operators import EndomorphicOperator,\  Theo Steininger committed Jan 27, 2017 20 21  FFTOperator,\ InvertibleOperatorMixin  Theo Steininger committed Sep 17, 2016 22 23   Theo Steininger committed Jan 27, 2017 24 class PropagatorOperator(InvertibleOperatorMixin, EndomorphicOperator):  Theo Steininger committed May 13, 2017 25  """ NIFTY Propagator Operator D.  Theo Steininger committed Sep 17, 2016 26   Pumpe, Daniel (dpumpe) committed May 08, 2017 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52  The propagator operator D, is known from the Wiener Filter. Its inverse functional form might look like: D = (S^(-1) + M)^(-1) D = (S^(-1) + N^(-1))^(-1) D = (S^(-1) + R^(\dagger) N^(-1) R)^(-1) Parameters ---------- S : LinearOperator Covariance of the signal prior. M : LinearOperator Likelihood contribution. R : LinearOperator Response operator translating signal to (noiseless) data. N : LinearOperator Covariance of the noise prior or the likelihood, respectively. inverter : class to invert explicitly defined operators (default:ConjugateGradient) preconditioner : Field numerical preconditioner to speed up convergence Attributes ---------- Raises ------  Pumpe, Daniel (dpumpe) committed May 09, 2017 53 54 55  ValueError is raised if * neither N nor M is given  Pumpe, Daniel (dpumpe) committed May 08, 2017 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75  Notes ----- Examples -------- >>> x_space = RGSpace(4) >>> k_space = RGRGTransformation.get_codomain(x_space) >>> f = Field(x_space, val=[2, 4, 6, 8]) >>> S = create_power_operator(k_space, spec=1) >>> N = DiagonalOperaor(f.domain, diag=1) >>> D = PropagatorOperator(S=S, N=N) # D^{-1} = S^{-1} + N^{-1} >>> D(f).val array([ 1., 2., 3., 4.] See Also -------- Scientific reference https://arxiv.org/abs/0806.3474  Theo Steininger committed May 13, 2017 76   Pumpe, Daniel (dpumpe) committed May 08, 2017 77 78  """  Theo Steininger committed Sep 17, 2016 79 80  # ---Overwritten properties and methods---  Theo Steininger committed Oct 14, 2016 81  def __init__(self, S=None, M=None, R=None, N=None, inverter=None,  Theo Steininger committed May 10, 2017 82  preconditioner=None, default_spaces=None):  Theo Steininger committed Sep 17, 2016 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98  """ Sets the standard operator properties and codomain, _A1, _A2, and RN if required. Parameters ---------- S : operator Covariance of the signal prior. M : operator Likelihood contribution. R : operator Response operator translating signal to (noiseless) data. N : operator Covariance of the noise prior or the likelihood, respectively. """  Theo Steininger committed Oct 14, 2016 99 100 101 102 103 104 105 106 107 108 109  # infer domain, and target if M is not None: self._domain = M.domain self._likelihood_times = M.times elif N is None: raise ValueError("Either M or N must be given!") elif R is not None: self._domain = R.domain self._likelihood_times = \  110  lambda z: R.adjoint_times(N.inverse_times(R.times(z)))  Theo Steininger committed Oct 14, 2016 111  else:  112 113  self._domain = N.domain self._likelihood_times = lambda z: N.inverse_times(z)  Theo Steininger committed Sep 17, 2016 114   Theo Steininger committed Nov 02, 2016 115 116  self._S = S self._fft_S = FFTOperator(self._domain, target=self._S.domain)  Theo Steininger committed Oct 14, 2016 117 118 119 120  if preconditioner is None: preconditioner = self._S_times  Theo Steininger committed Jan 27, 2017 121  super(PropagatorOperator, self).__init__(inverter=inverter,  Theo Steininger committed May 10, 2017 122 123  preconditioner=preconditioner, default_spaces=default_spaces)  Theo Steininger committed Sep 17, 2016 124 125 126 127 128  # ---Mandatory properties and methods--- @property def domain(self):  Theo Steininger committed Oct 14, 2016 129  return self._domain  Theo Steininger committed Sep 17, 2016 130 131  @property  Martin Reinecke committed Apr 28, 2017 132  def self_adjoint(self):  Theo Steininger committed Sep 17, 2016 133 134 135 136 137 138 139 140  return True @property def unitary(self): return False # ---Added properties and methods---  Theo Steininger committed Feb 07, 2017 141 142 143 144  def _S_times(self, x, spaces=None): transformed_x = self._fft_S(x, spaces=spaces) y = self._S(transformed_x, spaces=spaces) transformed_y = self._fft_S.inverse_times(y, spaces=spaces)  Theo Steininger committed Nov 02, 2016 145 146 147 148  result = x.copy_empty() result.set_val(transformed_y, copy=False) return result  Theo Steininger committed Feb 07, 2017 149 150 151 152  def _S_inverse_times(self, x, spaces=None): transformed_x = self._fft_S(x, spaces=spaces) y = self._S.inverse_times(transformed_x, spaces=spaces) transformed_y = self._fft_S.inverse_times(y, spaces=spaces)  Theo Steininger committed Nov 02, 2016 153 154 155 156  result = x.copy_empty() result.set_val(transformed_y, copy=False) return result  Theo Steininger committed Feb 07, 2017 157 158  def _inverse_times(self, x, spaces): pre_result = self._S_inverse_times(x, spaces)  Theo Steininger committed Nov 02, 2016 159 160 161  pre_result += self._likelihood_times(x) result = x.copy_empty() result.set_val(pre_result, copy=False)  Theo Steininger committed Sep 17, 2016 162  return result