diff --git a/demos/.DS_Store b/demos/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6
Binary files /dev/null and b/demos/.DS_Store differ
diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py
index df291340f2e75aba4b98ccadd7a807ee716ca67e..6d157fc1a6ca03790114d10797f4dfbb026dce78 100644
--- a/demos/critical_filtering.py
+++ b/demos/critical_filtering.py
@@ -8,33 +8,33 @@ from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
-np.random.seed(62)
+np.random.seed(232)
-class NonlinearResponse(LinearOperator):
- def __init__(self, FFT, Instrument, function, derivative, default_spaces=None):
- super(NonlinearResponse, self).__init__(default_spaces)
+
+def plot_parameters(m,t,t_true, t_real, t_d):
+ m = fft.adjoint_times(m)
+ m_data = m.val.get_full_data().real
+ t_data = t.val.get_full_data().real
+ t_d_data = t_d.val.get_full_data().real
+ t_true_data = t_true.val.get_full_data().real
+ t_real_data = t_real.val.get_full_data().real
+ pl.plot([go.Heatmap(z=m_data)], filename='map.html')
+ pl.plot([go.Scatter(y=t_data), go.Scatter(y=t_true_data),
+ go.Scatter(y=t_real_data), go.Scatter(y=t_d_data)], filename="t.html")
+
+class AdjointFFTResponse(LinearOperator):
+ def __init__(self, FFT, R, default_spaces=None):
+ super(AdjointFFTResponse, self).__init__(default_spaces)
self._domain = FFT.target
- self._target = Instrument.target
- self.function = function
- self.derivative = derivative
- self.I = Instrument
+ self._target = R.target
+ self.R = R
self.FFT = FFT
-
def _times(self, x, spaces=None):
- return self.I(self.function(self.FFT.adjoint_times(x)))
+ return self.R(self.FFT.adjoint_times(x))
def _adjoint_times(self, x, spaces=None):
- return self.FFT(self.function(self.I.adjoint_times(x)))
-
- def derived_times(self, x, position):
- position_derivative = self.derivative(self.FFT.adjoint_times(position))
- return self.I(position_derivative * self.FFT.adjoint_times(x))
-
- def derived_adjoint_times(self, x, position):
- position_derivative = self.derivative(self.FFT.adjoint_times(position))
- return self.FFT(position_derivative * self.I.adjoint_times(x))
-
+ return self.FFT(self.R.adjoint_times(x))
@property
def domain(self):
return self._domain
@@ -47,16 +47,6 @@ class NonlinearResponse(LinearOperator):
def unitary(self):
return False
-def plot_parameters(m,t,t_true, t_real):
- m = fft.adjoint_times(m)
- m_data = m.val.get_full_data().real
- t_data = t.val.get_full_data().real
- t_true_data = t_true.val.get_full_data().real
- t_real_data = t_real.val.get_full_data().real
- pl.plot([go.Heatmap(z=m_data)], filename='map.html')
- pl.plot([go.Scatter(y=t_data), go.Scatter(y=t_true_data),
- go.Scatter(y=t_real_data)], filename="t.html")
-
if __name__ == "__main__":
distribution_strategy = 'not'
@@ -70,11 +60,11 @@ if __name__ == "__main__":
h_space = fft.target[0]
# Setting up power space
- p_space = PowerSpace(h_space, logarithmic = True,
- distribution_strategy=distribution_strategy)
+ p_space = PowerSpace(h_space, logarithmic=False,
+ distribution_strategy=distribution_strategy, nbin=128)
# Choosing the prior correlation structure and defining correlation operator
- pow_spec = (lambda k: (.05 / (k + 1) ** 3))
+ pow_spec = (lambda k: (.05 / (k + 1) ** 2))
S = create_power_operator(h_space, power_spectrum=pow_spec,
distribution_strategy=distribution_strategy)
@@ -89,43 +79,11 @@ if __name__ == "__main__":
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
- # Choosing nonlinearity
-
- # log-normal model:
-
- function = exp
- derivative = exp
-
- # tan-normal model
-
- # def function(x):
- # return 0.5 * tanh(x) + 0.5
- # def derivative(x):
- # return 0.5*(1 - tanh(x)**2)
-
- # no nonlinearity, Wiener Filter
-
- # def function(x):
- # return x
- # def derivative(x):
- # return 1
-
- # small quadratic pertubarion
-
- # def function(x):
- # return 0.5*x**2 + x
- # def derivative(x):
- # return x + 1
-
- # def function(x):
- # return 0.9*x**4 +0.2*x**2 + x
- # def derivative(x):
- # return 0.9*4*x**3 + 0.4*x +1
- #
#Adding a harmonic transformation to the instrument
- R = NonlinearResponse(fft, Instrument, function, derivative)
- noise = .01
+ R = AdjointFFTResponse(fft, Instrument)
+
+ noise = 1.
N = DiagonalOperator(s_space, diagonal=noise, bare=True)
n = Field.from_random(domain=s_space,
random_type='normal',
@@ -134,66 +92,64 @@ if __name__ == "__main__":
# Creating the mock data
d = R(sh) + n
- realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"])**2)
+
+ realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"],
+ nbin=p_space.config["nbin"])**2)
+ data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
+ nbin=p_space.config["nbin"])**2)
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html')
- # Choosing the minimization strategy
+ # minimization strategy
def convergence_measure(a_energy, iteration): # returns current energy
x = a_energy.value
print (x, iteration)
- # minimizer1 = SteepestDescent(convergence_tolerance=0,
- # iteration_limit=50,
- # callback=convergence_measure)
minimizer1 = RelaxedNewton(convergence_tolerance=0,
convergence_level=1,
- iteration_limit=5,
+ iteration_limit=2,
callback=convergence_measure)
- # minimizer2 = RelaxedNewton(convergence_tolerance=0,
- # convergence_level=1,
- # iteration_limit=2,
- # callback=convergence_measure)
- #
- # minimizer1 = VL_BFGS(convergence_tolerance=0,
- # iteration_limit=5,
- # callback=convergence_measure,
- # max_history_length=3)
-
-
+ minimizer2 = VL_BFGS(convergence_tolerance=0,
+ iteration_limit=50,
+ callback=convergence_measure,
+ max_history_length=3)
# Setting starting position
flat_power = Field(p_space,val=10e-8)
m0 = flat_power.power_synthesize(real_signal=True)
t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
- # t0 = Field(p_space,val=-8)
- # t0 = log(sp.copy()**2)
+ #t0 = data_power
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
-
for i in range(100):
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
# Initializing the nonlinear Wiener Filter energy
- map_energy = NonlinearWienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
+ map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
# Minimization with chosen minimizer
- (map_energy, convergence) = minimizer1(map_energy)
+ map_energy = map_energy.analytic_solution()
+
# Updating parameters for correlation structure reconstruction
m0 = map_energy.position
D0 = map_energy.curvature
# Initializing the power energy with updated parameters
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=10., samples=3)
- (power_energy, convergence) = minimizer1(power_energy)
+ if i > 0:
+ (power_energy, convergence) = minimizer1(power_energy)
+ else:
+ (power_energy, convergence) = minimizer2(power_energy)
# Setting new power spectrum
t0 = power_energy.position
- plot_parameters(m0,t0,log(sp**2),realized_power)
+ t0.val[-1] = t0.val[-2]
+ # Plotting current estimate
+ plot_parameters(m0,t0,log(sp**2),realized_power, data_power)
# Transforming fields to position space for plotting
@@ -223,10 +179,3 @@ if __name__ == "__main__":
m_data = m.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=m_data)], filename='map.html')
-
- f_m_data = function(m).val.get_full_data().real
- if rank == 0:
- pl.plot([go.Heatmap(z=f_m_data)], filename='f_map.html')
- f_ss_data = function(ss).val.get_full_data().real
- if rank == 0:
- pl.plot([go.Heatmap(z=f_ss_data)], filename='f_ss.html')
diff --git a/demos/laplace_testing.py b/demos/laplace_testing.py
new file mode 100644
index 0000000000000000000000000000000000000000..be878e2feddcdf991bfef9f325ee3c098bea124b
--- /dev/null
+++ b/demos/laplace_testing.py
@@ -0,0 +1,105 @@
+from nifty import *
+
+import numpy as np
+from nifty import Field,\
+ EndomorphicOperator,\
+ PowerSpace
+import plotly.offline as pl
+import plotly.graph_objs as go
+
+import numpy as np
+from nifty import Field, \
+ EndomorphicOperator, \
+ PowerSpace
+
+
+class TestEnergy(Energy):
+ def __init__(self, position, Op):
+ super(TestEnergy, self).__init__(position)
+ self.Op = Op
+
+ def at(self, position):
+ return self.__class__(position=position, Op=self.Op)
+
+ @property
+ def value(self):
+ return 0.5 * self.position.dot(self.Op(self.position))
+
+ @property
+ def gradient(self):
+ return self.Op(self.position)
+
+ @property
+ def curvature(self):
+ curv = CurvOp(self.Op)
+ return curv
+
+class CurvOp(InvertibleOperatorMixin, EndomorphicOperator):
+ def __init__(self, Op,inverter=None, preconditioner=None):
+ self.Op = Op
+ self._domain = Op.domain
+ super(CurvOp, self).__init__(inverter=inverter,
+ preconditioner=preconditioner)
+
+ def _times(self, x, spaces):
+ return self.Op(x)
+
+if __name__ == "__main__":
+
+ distribution_strategy = 'not'
+
+ # Set up position space
+ s_space = RGSpace([128,128])
+ # s_space = HPSpace(32)
+
+ # Define harmonic transformation and associated harmonic space
+ fft = FFTOperator(s_space)
+ h_space = fft.target[0]
+
+ # Setting up power space
+ p_space = PowerSpace(h_space, logarithmic=False,
+ distribution_strategy=distribution_strategy, nbin=16)
+
+ # Choosing the prior correlation structure and defining correlation operator
+ pow_spec = (lambda k: (.05 / (k + 1) ** 2))
+ # t = Field(p_space, val=pow_spec)
+ t= Field.from_random("normal", domain=p_space)
+ lap = LogLaplaceOperator(p_space)
+ T = SmoothnessOperator(p_space,sigma=1.)
+ test_energy = TestEnergy(t,T)
+
+ def convergence_measure(a_energy, iteration): # returns current energy
+ x = a_energy.value
+ print (x, iteration)
+ minimizer1 = VL_BFGS(convergence_tolerance=0,
+ iteration_limit=1000,
+ callback=convergence_measure,
+ max_history_length=3)
+
+
+ def explicify(op, domain):
+ space = domain
+ d = space.dim
+ res = np.zeros((d, d))
+ for i in range(d):
+ x = np.zeros(d)
+ x[i] = 1.
+ f = Field(space, val=x)
+ res[:, i] = op(f).val
+ return res
+ A = explicify(lap.times, p_space)
+ B = explicify(lap.adjoint_times, p_space)
+ test_energy,convergence = minimizer1(test_energy)
+ data = test_energy.position.val.get_full_data()
+ pl.plot([go.Scatter(x=log(p_space.kindex)[1:], y=data[1:])], filename="t.html")
+ tt = Field.from_random("normal", domain=t.domain)
+ print "adjointness"
+ print t.dot(lap(tt))
+ print tt.dot(lap.adjoint_times(t))
+ print "log kindex"
+ aa = Field(p_space, val=p_space.kindex.copy())
+ aa.val[0] = 1
+
+ print lap(log(aa)).val
+ print "######################"
+ print test_energy.position.val
\ No newline at end of file
diff --git a/demos/nonlinear_wiener_filter.py b/demos/nonlinear_wiener_filter.py
deleted file mode 100644
index 355d791257027f485eefc79e8c8216f4c4f80e46..0000000000000000000000000000000000000000
--- a/demos/nonlinear_wiener_filter.py
+++ /dev/null
@@ -1,182 +0,0 @@
-
-from nifty import *
-
-import plotly.offline as pl
-import plotly.graph_objs as go
-
-from mpi4py import MPI
-comm = MPI.COMM_WORLD
-rank = comm.rank
-
-np.random.seed(42)
-
-class NonlinearResponse(LinearOperator):
- def __init__(self, FFT, Instrument, function, derivative, default_spaces=None):
- super(NonlinearResponse, self).__init__(default_spaces)
- self._domain = FFT.target
- self._target = Instrument.target
- self.function = function
- self.derivative = derivative
- self.I = Instrument
- self.FFT = FFT
-
-
- def _times(self, x, spaces=None):
- return self.I(self.function(self.FFT.adjoint_times(x)))
-
- def _adjoint_times(self, x, spaces=None):
- return self.FFT(self.function(self.I.adjoint_times(x)))
-
- def derived_times(self, x, position):
- position_derivative = self.derivative(self.FFT.adjoint_times(position))
- return self.I(position_derivative * self.FFT.adjoint_times(x))
-
- def derived_adjoint_times(self, x, position):
- position_derivative = self.derivative(self.FFT.adjoint_times(position))
- return self.FFT(position_derivative * self.I.adjoint_times(x))
-
- @property
- def domain(self):
- return self._domain
-
- @property
- def target(self):
- return self._target
-
- @property
- def unitary(self):
- return False
-
-
-
-if __name__ == "__main__":
-
- distribution_strategy = 'not'
-
- # Set up position space
- s_space = RGSpace([100,100])
- # s_space = HPSpace(32)
-
- # Define harmonic transformation and associated harmonic space
- fft = FFTOperator(s_space)
- h_space = fft.target[0]
-
- # Setting up power space
- p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
-
- # Choosing the prior correlation structure and defining correlation operator
- pow_spec = (lambda k: (1. / (k + 1) ** 3))
- S = create_power_operator(h_space, power_spectrum=pow_spec,
- distribution_strategy=distribution_strategy)
-
- # Drawing a sample sh from the prior distribution in harmonic space
- sp = Field(p_space, val=lambda z: pow_spec(z)**(1./2),
- distribution_strategy=distribution_strategy)
- sh = sp.power_synthesize(real_signal=True)
-
-
- # Choosing the measurement instrument
-# Instrument = SmoothingOperator(s_space, sigma=0.01)
- Instrument = DiagonalOperator(s_space, diagonal=1.)
-# Instrument._diagonal.val[200:400, 200:400] = 0
-
- # Choosing nonlinearity
- #
- # function = exp
- # derivative = exp
- def function(x):
- return 0.5 * tanh(x) + 0.5
- def derivative(x):
- return 0.5*(1. - tanh(x)**2)
- # def function(x):
- # return 20*x
- # def derivative(x):
- # return 20.
-
- # def function(x):
- # return 0.1*0.5*x**2 + x
- # def derivative(x):
- # return 0.1*x + 1
-
- # def function(x):
- # return 0.9*x**4 +0.2*x**2 + x
- # def derivative(x):
- # return 0.9*4*x**3 + 0.4*x +1
- #
-
- #Adding a harmonic transformation to the instrument
-
- noise = 10.
- N = DiagonalOperator(s_space, diagonal=noise, bare=True)
- n = Field.from_random(domain=s_space,
- random_type='normal',
- std=sqrt(noise),
- mean=0)
- R = NonlinearResponse(fft, Instrument, function, derivative)
- # Creating the mock data
- d = R(sh) + n
-
- # Choosing the minimization strategy
-
- def convergence_measure(energy, iteration): # returns current energy
- x = energy.value
- print (x, iteration)
-
-# minimizer = SteepestDescent(convergence_tolerance=0,
-# iteration_limit=50,
-# callback=convergence_measure)
-#
- minimizer = RelaxedNewton(convergence_tolerance=0,
- convergence_level=1,
- iteration_limit=3,
- callback=convergence_measure)
-
- # minimizer = VL_BFGS(convergence_tolerance=0,
- # iteration_limit=50,
- # callback=convergence_measure,
- # max_history_length=3)
- #
-
-
- # Setting starting position
- flat_power = Field(p_space,val=10e-8)
- m0 = flat_power.power_synthesize(real_signal=True)
-
- # Initializing the Wiener Filter energy
- energy = NonlinearWienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S)
-
- # Solving the problem with chosen minimization strategy
- (energy, convergence) = minimizer(energy)
-
- # Transforming fields to position space for plotting
-
- ss = fft.adjoint_times(sh)
- m = fft.adjoint_times(energy.position)
-
-
- # Plotting
-
- d_data = d.val.get_full_data().real
- if rank == 0:
- pl.plot([go.Heatmap(z=d_data)], filename='data.html')
-
-
- ss_data = ss.val.get_full_data().real
- if rank == 0:
- pl.plot([go.Heatmap(z=ss_data)], filename='ss.html')
-
- sh_data = sh.val.get_full_data().real
- if rank == 0:
- pl.plot([go.Heatmap(z=sh_data)], filename='sh.html')
-
-
- m_data = m.val.get_full_data().real
- if rank == 0:
- pl.plot([go.Heatmap(z=m_data)], filename='map.html')
-
- f_m_data = function(m).val.get_full_data().real
- if rank == 0:
- pl.plot([go.Heatmap(z=f_m_data)], filename='f_map.html')
- f_ss_data = function(ss).val.get_full_data().real
- if rank == 0:
- pl.plot([go.Heatmap(z=f_ss_data)], filename='f_ss.html')
diff --git a/demos/train_example.py b/demos/train_example.py
index 2b44223d792f1fc18823301e0c9117284948395d..07ff6c2b482b433003bb781de723db2b7efed69d 100644
--- a/demos/train_example.py
+++ b/demos/train_example.py
@@ -54,12 +54,36 @@ def plot_parameters(m,t, t_real):
pl.plot([go.Scatter(y=m_data)], filename='map.html')
pl.plot([go.Scatter(y=t_data),
go.Scatter(y=t_real_data)], filename="t.html")
+class AdjointFFTResponse(LinearOperator):
+ def __init__(self, FFT, R, default_spaces=None):
+ super(AdjointFFTResponse, self).__init__(default_spaces)
+ self._domain = FFT.target
+ self._target = R.target
+ self.R = R
+ self.FFT = FFT
+
+ def _times(self, x, spaces=None):
+ return self.R(self.FFT.adjoint_times(x))
+
+ def _adjoint_times(self, x, spaces=None):
+ return self.FFT(self.R.adjoint_times(x))
+ @property
+ def domain(self):
+ return self._domain
+
+ @property
+ def target(self):
+ return self._target
+
+ @property
+ def unitary(self):
+ return False
if __name__ == "__main__":
distribution_strategy = 'not'
full_data = np.genfromtxt("train_data.csv", delimiter = ',')
- d = full_data.T[2]
+ d = full_data.T[3]
d[0] = 0.
d -= d.mean()
d[0] = 0.
@@ -80,42 +104,9 @@ if __name__ == "__main__":
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
- # Choosing nonlinearity
-
- # log-normal model:
-
- # function = exp
- # derivative = exp
-
- # tan-normal model
-
- # def function(x):
- # return 0.5 * tanh(x) + 0.5
- # def derivative(x):
- # return 0.5*(1 - tanh(x)**2)
-
- # no nonlinearity, Wiener Filter
-
- def function(x):
- return x
- def derivative(x):
- return 1
-
- # small quadratic pertubarion
-
- # def function(x):
- # return 0.5*x**2 + x
- # def derivative(x):
- # return x + 1
-
- # def function(x):
- # return 0.9*x**4 +0.2*x**2 + x
- # def derivative(x):
- # return 0.9*4*x**3 + 0.4*x +1
- #
#Adding a harmonic transformation to the instrument
- R = NonlinearResponse(fft, Instrument, function, derivative)
+ R = AdjointFFTResponse(fft, Instrument)
noise = .1
N = DiagonalOperator(s_space, diagonal=noise, bare=True)
@@ -139,13 +130,13 @@ if __name__ == "__main__":
callback=convergence_measure)
minimizer2 = RelaxedNewton(convergence_tolerance=0,
convergence_level=1,
- iteration_limit=10,
+ iteration_limit=30,
callback=convergence_measure)
- # minimizer1 = VL_BFGS(convergence_tolerance=0,
- # iteration_limit=5,
- # callback=convergence_measure,
- # max_history_length=3)
+ minimizer1 = VL_BFGS(convergence_tolerance=0,
+ iteration_limit=30,
+ callback=convergence_measure,
+ max_history_length=3)
@@ -159,25 +150,31 @@ if __name__ == "__main__":
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
-
- data_power = fft(d).power_analyze()
+ data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
+ nbin=p_space.config["nbin"])**2)
for i in range(100):
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
# Initializing the nonlinear Wiener Filter energy
- map_energy = NonlinearWienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
+ map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
# Minimization with chosen minimizer
- (map_energy, convergence) = minimizer1(map_energy)
+ map_energy = map_energy.analytic_solution()
+
# Updating parameters for correlation structure reconstruction
m0 = map_energy.position
D0 = map_energy.curvature
# Initializing the power energy with updated parameters
- power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=.5, samples=5)
- (power_energy, convergence) = minimizer1(power_energy)
+ power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=.5, samples=3)
+ if i > 0:
+ (power_energy, convergence) = minimizer1(power_energy)
+ else:
+ (power_energy, convergence) = minimizer1(power_energy)
# Setting new power spectrum
t0 = power_energy.position
- plot_parameters(m0,t0,log(data_power**2))
+ t0.val[-1] = t0.val[-2]
+ # Plotting current estimate
+ plot_parameters(m0,t0,data_power)
# Transforming fields to position space for plotting
@@ -185,6 +182,7 @@ if __name__ == "__main__":
m = fft.adjoint_times(map_energy.position)
+
# Plotting
d_data = d.val.get_full_data().real
diff --git a/demos/wiener_filter_advanced.py b/demos/wiener_filter_advanced.py
index e9660e682241b826cc8008e477be769a99d4b9a0..c98f2b58d668fb6a182ef1d3dfdf65a49cf5bd7b 100644
--- a/demos/wiener_filter_advanced.py
+++ b/demos/wiener_filter_advanced.py
@@ -61,7 +61,7 @@ if __name__ == "__main__":
sp = Field(p_space, val=lambda z: pow_spec(z)**(1./2),
distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True)
-
+ ss = fft.adjoint_times(sh)
# Choosing the measurement instrument
Instrument = SmoothingOperator(s_space, sigma=0.05)
diff --git a/nifty/field.py b/nifty/field.py
index 2f00b96565bcd930d09ff5c442baa22ee7ca2fa9..fa25faacc61472f91e1296527ae6608765d73c57 100644
--- a/nifty/field.py
+++ b/nifty/field.py
@@ -167,110 +167,20 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods---
- def power_analyze(self, spaces=None, logarithmic=False, nbin=None, binbounds=None,
- decompose_power=False):
- # check if all spaces in `self.domain` are either harmonic or
- # power_space instances
- for sp in self.domain:
- if not sp.harmonic and not isinstance(sp, PowerSpace):
- self.logger.info(
- "Field has a space in `domain` which is neither "
- "harmonic nor a PowerSpace.")
-
- # check if the `spaces` input is valid
- spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
- if spaces is None:
- if len(self.domain) == 1:
- spaces = (0,)
- else:
- raise ValueError(
- "Field has multiple spaces as domain "
- "but `spaces` is None.")
-
- if len(spaces) == 0:
- raise ValueError(
- "No space for analysis specified.")
- elif len(spaces) > 1:
- raise ValueError(
- "Conversion of only one space at a time is allowed.")
-
- space_index = spaces[0]
-
- if not self.domain[space_index].harmonic:
- raise ValueError(
- "The analyzed space must be harmonic.")
-
- # Create the target PowerSpace instance:
- # If the associated signal-space field was real, we extract the
- # hermitian and anti-hermitian parts of `self` and put them
- # into the real and imaginary parts of the power spectrum.
- # If it was complex, all the power is put into a real power spectrum.
-
- distribution_strategy = \
- self.val.get_axes_local_distribution_strategy(
- self.domain_axes[space_index])
-
- harmonic_partner = self.domain[space_index]
- power_domain = PowerSpace(harmonic_partner=harmonic_partner,
- distribution_strategy=distribution_strategy,
- logarithmic=logarithmic, nbin=nbin, binbounds=binbounds)
-
- # extract pindex and rho from power_domain
- pindex = power_domain.pindex
- rho = power_domain.rho
-
- if decompose_power:
- hermitian_part, anti_hermitian_part = \
- harmonic_partner.hermitian_decomposition(
- self.val,
- axes=self.domain_axes[space_index])
-
- [hermitian_power, anti_hermitian_power] = \
- [self._calculate_power_spectrum(
- x=part,
- pindex=pindex,
- rho=rho,
- axes=self.domain_axes[space_index])
- for part in [hermitian_part, anti_hermitian_part]]
-
- power_spectrum = hermitian_power + 1j * anti_hermitian_power
- else:
- power_spectrum = self._calculate_power_spectrum(
- x=self.val,
- pindex=pindex,
- rho=rho,
- axes=self.domain_axes[space_index])
-
- # create the result field and put power_spectrum into it
- result_domain = list(self.domain)
- result_domain[space_index] = power_domain
-
- if decompose_power:
- result_dtype = np.complex
- else:
- result_dtype = np.float
-
- result_field = self.copy_empty(
- domain=result_domain,
- dtype=result_dtype,
- distribution_strategy=power_spectrum.distribution_strategy)
- result_field.set_val(new_val=power_spectrum, copy=False)
-
- return result_field
-
- def project_power(self, spaces=None, logarithmic=False, nbin=None,
+ def power_analyze(self, spaces=None, logarithmic=False, nbin=None,
binbounds=None, decompose_power=True):
- """ Computes the quadratic power projection for a subspace of the Field.
+ """ Computes the square root power spectrum for a subspace of `self`.
Creates a PowerSpace for the space addressed by `spaces` with the given
- binning and computes the quadratic power projection as a Field over this
+ binning and computes the power spectrum as a Field over this
PowerSpace. This can only be done if the subspace to be analyzed is a
- harmonic space.
+ harmonic space. The resulting field has the same units as the initial
+ field, corresponding to the square root of the power spectrum.
Parameters
----------
spaces : int *optional*
- The subspace for which the power projection shall be computed
+ The subspace for which the powerspectrum shall be computed
(default : None).
if spaces==None : Tries to synthesize for the whole domain
logarithmic : boolean *optional*
@@ -285,7 +195,7 @@ class Field(Loggable, Versionable, object):
if binbounds==None : bins are inferred. Overwrites nbins and log
decompose_power : boolean, *optional*
Whether the analysed signal-space Field is intrinsically real or
- complex and if the projected power shall therefore be computed
+ complex and if the power spectrum shall therefore be computed
for the real and the imaginary part of the Field separately
(default : True).
@@ -301,11 +211,11 @@ class Field(Loggable, Versionable, object):
-------
out : Field
The output object. It's domain is a PowerSpace and it contains
- the quadratic power projection of 'self's field.
+ the power spectrum of 'self's field.
See Also
--------
- power_synthesize, PowerSpace, power_analyze
+ power_synthesize, PowerSpace
"""
@@ -343,8 +253,8 @@ class Field(Loggable, Versionable, object):
# Create the target PowerSpace instance:
# If the associated signal-space field was real, we extract the
# hermitian and anti-hermitian parts of `self` and put them
- # into the real and imaginary parts of the projected power.
- # If it was complex, all the power is put into a real power projection
+ # into the real and imaginary parts of the power spectrum.
+ # If it was complex, all the power is put into a real power spectrum.
distribution_strategy = \
self.val.get_axes_local_distribution_strategy(
@@ -358,6 +268,7 @@ class Field(Loggable, Versionable, object):
# extract pindex and rho from power_domain
pindex = power_domain.pindex
+ rho = power_domain.rho
if decompose_power:
hermitian_part, anti_hermitian_part = \
@@ -366,37 +277,36 @@ class Field(Loggable, Versionable, object):
axes=self.domain_axes[space_index])
[hermitian_power, anti_hermitian_power] = \
- [self._calculate_projected_power(
+ [self._calculate_power_spectrum(
x=part,
pindex=pindex,
+ rho=rho,
axes=self.domain_axes[space_index])
for part in [hermitian_part, anti_hermitian_part]]
- power_projection = hermitian_power + 1j * anti_hermitian_power
+ power_spectrum = hermitian_power + 1j * anti_hermitian_power
+
else:
- power_projection = self._calculate_projected_power(
+ power_spectrum = self._calculate_power_spectrum(
x=self.val,
pindex=pindex,
+ rho=rho,
axes=self.domain_axes[space_index])
# create the result field and put power_spectrum into it
result_domain = list(self.domain)
result_domain[space_index] = power_domain
-
- if decompose_power:
- result_dtype = np.complex
- else:
- result_dtype = np.float
+ result_dtype = power_spectrum.dtype
result_field = self.copy_empty(
domain=result_domain,
dtype=result_dtype,
- distribution_strategy=power_projection.distribution_strategy)
- result_field.set_val(new_val=power_projection, copy=False)
+ distribution_strategy=power_spectrum.distribution_strategy)
+ result_field.set_val(new_val=power_spectrum, copy=False)
return result_field
- def _calculate_projected_power(self, x, pindex, axes=None):
+ def _calculate_power_spectrum(self, x, pindex, rho, axes=None):
fieldabs = abs(x)
fieldabs **= 2
@@ -406,12 +316,8 @@ class Field(Loggable, Versionable, object):
target_shape=x.shape,
target_strategy=x.distribution_strategy,
axes=axes)
- projected_power = pindex.bincount(weights=fieldabs,
+ power_spectrum = pindex.bincount(weights=fieldabs,
axis=axes)
- return projected_power
-
- def _calculate_power_spectrum(self, x, pindex, rho, axes=None):
- power_spectrum = self._calculate_projected_power(x, pindex, axes)
if axes is not None:
new_rho_shape = [1, ] * len(power_spectrum.shape)
new_rho_shape[axes[0]] = len(rho)
diff --git a/nifty/library/.DS_Store b/nifty/library/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..b81c8dc7087a9ee1d969b9583e404e257c710f29
Binary files /dev/null and b/nifty/library/.DS_Store differ
diff --git a/nifty/library/energy_library/.DS_Store b/nifty/library/energy_library/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6
Binary files /dev/null and b/nifty/library/energy_library/.DS_Store differ
diff --git a/nifty/library/energy_library/__init__.py b/nifty/library/energy_library/__init__.py
index 2cc72dd3a77bee8514f6f3025924afeea037fc1c..381bf16434f1c2d8b194c34d4b7df82a8754fe8f 100644
--- a/nifty/library/energy_library/__init__.py
+++ b/nifty/library/energy_library/__init__.py
@@ -1,3 +1,2 @@
from wiener_filter_energy import WienerFilterEnergy
-from nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
from critical_power_energy import CriticalPowerEnergy
\ No newline at end of file
diff --git a/nifty/library/energy_library/critical_power_energy.py b/nifty/library/energy_library/critical_power_energy.py
index 2aee44df934edefdbf85e03aa63c25d6dcd60804..81c75c005d438efc2edc74f7c25022bf00438bf4 100644
--- a/nifty/library/energy_library/critical_power_energy.py
+++ b/nifty/library/energy_library/critical_power_energy.py
@@ -1,6 +1,7 @@
from nifty.energies.energy import Energy
from nifty.library.operator_library import CriticalPowerCurvature,\
- SmoothnessOperator
+ LaplaceOperator
+
from nifty.sugar import generate_posterior_sample
from nifty import Field, exp
@@ -30,7 +31,9 @@ class CriticalPowerEnergy(Energy):
self.sigma = sigma
self.alpha = Field(self.position.domain, val=alpha)
self.q = Field(self.position.domain, val=q)
- self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma)
+ #self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma)
+ self.Laplace = LaplaceOperator(self.position.domain[0])
+ self.roughness = self.Laplace(self.position)
self.rho = self.position.domain[0].rho
self.w = w
if self.w is None:
@@ -47,21 +50,21 @@ class CriticalPowerEnergy(Energy):
@property
def value(self):
energy = exp(-self.position).dot(self.q + self.w / 2.)
- energy += self.position.dot(self.alpha - 1 + self.rho / 2.)
- energy += 0.5 * self.position.dot(self.T(self.position))
+ energy += self.position.dot(self.alpha - 1. + self.rho / 2.)
+ energy += 0.5 * self.roughness.dot(self.roughness) / self.sigma**2
return energy.real
@property
def gradient(self):
gradient = - self.theta
- gradient += self.alpha - 1 + self.rho / 2.
- gradient += self.T(self.position)
- gradient.val[0] = 0.
+ gradient += self.alpha - 1. + self.rho / 2.
+ gradient += self.Laplace(self.roughness) / self.sigma **2
return gradient
@property
def curvature(self):
- curvature = CriticalPowerCurvature(theta=self.theta, T = self.T)
+ curvature = CriticalPowerCurvature(theta=self.theta, Laplace = self.Laplace,
+ sigma = self.sigma)
return curvature
def _calculate_w(self, m, D, samples):
@@ -69,15 +72,18 @@ class CriticalPowerEnergy(Energy):
if D is not None:
for i in range(samples):
posterior_sample = generate_posterior_sample(m, D)
- projected_sample = posterior_sample.project_power(
+ projected_sample = posterior_sample.power_analyze(
logarithmic=self.position.domain[0].config["logarithmic"],
+ nbin= self.position.domain[0].config["nbin"],
decompose_power=False)
- w += projected_sample
+ w += (projected_sample **2) * self.rho
w /= float(samples)
else:
- w = m.project_power(
+ w = m.power_analyze(
logarithmic=self.position.domain[0].config["logarithmic"],
decompose_power=False)
+ w **= 2
+ w *= self.rho
return w
diff --git a/nifty/library/energy_library/nonlinear_wiener_filter_energy.py b/nifty/library/energy_library/nonlinear_wiener_filter_energy.py
deleted file mode 100644
index 9cfd2c1ec6005622652abbaf2efd7da530bff26d..0000000000000000000000000000000000000000
--- a/nifty/library/energy_library/nonlinear_wiener_filter_energy.py
+++ /dev/null
@@ -1,53 +0,0 @@
-from nifty.energies.energy import Energy
-from nifty.library.operator_library import NonlinearWienerFilterCurvature
-from nifty import Field
-class NonlinearWienerFilterEnergy(Energy):
- """The Energy for the Gaussian lognormal case.
-
- It describes the situation of linear measurement of a
- lognormal signal with Gaussian noise and Gaussain signal prior.
-
- Parameters
- ----------
- d : Field,
- the data.
- R : Operator,
- The nonlinear response operator, describtion of the measurement process.
- N : EndomorphicOperator,
- The noise covariance in data space.
- S : EndomorphicOperator,
- The prior signal covariance in harmonic space.
- """
-
- def __init__(self, position, d, R, N, S):
- super(NonlinearWienerFilterEnergy, self).__init__(position = position)
- self.d = d
- self.R = R
- self.N = N
- self.S = S
-
-
- def at(self, position):
- return self.__class__(position, self.d, self.R, self.N, self.S)
-
- @property
- def value(self):
- energy = 0.5 * self.position.dot(self.S.inverse_times(self.position))
- energy += 0.5 * (self.d - self.R(self.position)).dot(
- self.N.inverse_times(self.d - self.R(self.position)))
- return energy.real
-
- @property
- def gradient(self):
- gradient = self.S.inverse_times(self.position)
- gradient -= self.R.derived_adjoint_times(
- self.N.inverse_times(self.d - self.R(self.position)), self.position)
- return gradient
-
- @property
- def curvature(self):
- curvature = NonlinearWienerFilterCurvature(R=self.R,
- N=self.N,
- S=self.S,
- position=self.position)
- return curvature
diff --git a/nifty/library/operator_library/.DS_Store b/nifty/library/operator_library/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6
Binary files /dev/null and b/nifty/library/operator_library/.DS_Store differ
diff --git a/nifty/library/operator_library/__init__.py b/nifty/library/operator_library/__init__.py
index b9a3cfa0a1796b883c581969094b0232311af618..64e43c39c6c21c4925e45a10fa295bfb8f91c3de 100644
--- a/nifty/library/operator_library/__init__.py
+++ b/nifty/library/operator_library/__init__.py
@@ -1,5 +1,5 @@
from wiener_filter_curvature import WienerFilterCurvature
-from nonlinear_wiener_filter_curvature import NonlinearWienerFilterCurvature
from critical_power_curvature import CriticalPowerCurvature
from laplace_operator import LaplaceOperator
+from laplace_operator import LogLaplaceOperator
from smoothness_operator import SmoothnessOperator
\ No newline at end of file
diff --git a/nifty/library/operator_library/critical_power_curvature.py b/nifty/library/operator_library/critical_power_curvature.py
index 5be1e2a6963defe6e2b4ba25dfc729ce9197e867..ece3cb58bb82a9fd6253b57998ea97a110a265f3 100644
--- a/nifty/library/operator_library/critical_power_curvature.py
+++ b/nifty/library/operator_library/critical_power_curvature.py
@@ -3,13 +3,14 @@ from nifty.operators import EndomorphicOperator,\
class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
- def __init__(self, theta, T, inverter=None, preconditioner=None):
+ def __init__(self, theta, Laplace, sigma, inverter=None, preconditioner=None):
self.theta = theta
- self.T = T
+ self.Laplace = Laplace
+ self.sigma = sigma
# if preconditioner is None:
# preconditioner = self.T.times
- self._domain = self.T.domain
+ self._domain = self.theta.domain
super(CriticalPowerCurvature, self).__init__(inverter=inverter,
preconditioner=preconditioner)
@property
@@ -27,4 +28,5 @@ class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
# ---Added properties and methods---
def _times(self, x, spaces):
- return self.T(x) + self.theta * x
+ return self.Laplace.adjoint_times(self.Laplace(x)) / self.sigma ** 2 \
+ + self.theta * x
diff --git a/nifty/library/operator_library/laplace_operator.py b/nifty/library/operator_library/laplace_operator.py
index 4ae03e2d599f6a119d75b99c4f210eac5ec7044c..98934775da199f8bc2a4a4a1028af20b687a0f8b 100644
--- a/nifty/library/operator_library/laplace_operator.py
+++ b/nifty/library/operator_library/laplace_operator.py
@@ -8,30 +8,37 @@ class LaplaceOperator(EndomorphicOperator):
default_spaces=None):
super(LaplaceOperator, self).__init__(default_spaces)
if (domain is not None):
- if(not isinstance(domain, PowerSpace)):
- raise TypeError("The domain has to live over a PowerSpace")
+ if (not isinstance(domain, PowerSpace)):
+ raise TypeError("The domain has to live over a PowerSpace")
self._domain = self._parse_domain(domain)
+
"""
input parameters:
domain- can only live over one domain
to do:
correct implementation of D20 object
"""
+
@property
def target(self):
return self._domain
+
@property
def domain(self):
return self._domain
+
@property
def unitary(self):
return False
+
@property
def symmetric(self):
return False
+
@property
def self_adjoint(self):
return False
+
def _times(self, t, spaces):
if t.val.distribution_strategy != 'not':
# self.logger.warn("distribution_strategy should be 'not' but is %s"
@@ -39,22 +46,106 @@ class LaplaceOperator(EndomorphicOperator):
array = t.val.get_full_data()
else:
array = t.val.get_local_data(copy=False)
- ret = 2 * array - np.append(0, array[:-1]) - np.append(array[1:], 0)
- ret[0] = 0
- ret[1] = 0
+ ret = 2 * array
+ ret -= np.roll(array, 1)
+ ret -= np.roll(array, -1)
+ ret[0:2] = 0
ret[-1] = 0
return Field(self.domain, val=ret)
+
def _adjoint_times(self, t, spaces):
+ t = t.copy().weight(1)
+ t.val[1:] /= self.domain[0].kindex[1:]
if t.val.distribution_strategy != 'not':
# self.logger.warn("distribution_strategy should be 'not' but is %s"
# % t.val.distribution_strategy)
array = t.val.get_full_data()
else:
array = t.val.get_local_data(copy=False)
- ret = 2 * array - np.append(0, array[:-1]) - np.append(array[1:], 0)
- ret[0] = 0
- ret[1] = -array[2]
- ret[2] = 2*array[2]-array[3]
- ret[-1] = -array[-2]
- ret[-2] = -array[-3]+2*array[-2]
- return Field(self.domain, val=ret)
\ No newline at end of file
+ ret = np.copy(array)
+ ret[0:2] = 0
+ ret[-1] = 0
+ ret = 2 * ret - np.roll(ret, 1) - np.roll(ret, -1)
+ result = Field(self.domain, val=ret).weight(-1)
+ return result
+
+
+
+def _irregular_nabla(x,k):
+ #applies forward differences and does nonesense at the edge. Thus needs cutting
+ y = -x
+ y[:-1] += x[1:]
+ y[1:-1] /= - k[1:-1] + k[2:]
+ return y
+
+def _irregular_adj_nabla(z, k):
+ #applies backwards differences*(-1) and does nonesense at the edge. Thus needs cutting
+ x = z.copy()
+ x[1:-1] /= - k[1:-1] + k[2:]
+ y = -x
+ y[1:] += x[:-1]
+ return y
+
+
+class LogLaplaceOperator(EndomorphicOperator):
+ def __init__(self, domain,
+ default_spaces=None):
+ super(LogLaplaceOperator, self).__init__(default_spaces)
+ if (domain is not None):
+ if (not isinstance(domain, PowerSpace)):
+ raise TypeError("The domain has to live over a PowerSpace")
+ self._domain = self._parse_domain(domain)
+
+ """
+ input parameters:
+ domain- can only live over one domain
+ to do:
+ correct implementation of D20 object
+ """
+
+ @property
+ def target(self):
+ return self._domain
+
+ @property
+ def domain(self):
+ return self._domain
+
+ @property
+ def unitary(self):
+ return False
+
+ @property
+ def symmetric(self):
+ return False
+
+ @property
+ def self_adjoint(self):
+ return False
+
+
+ def _times(self, t, spaces):
+ l_vec = self.domain[0].kindex.copy()
+ l_vec[1:] = np.log(l_vec[1:])
+ l_vec[0] = -1.
+ ret = t.val.copy()
+ # ret[2:] *= np.sqrt(np.log(k_vec)-np.log(np.roll(k_vec,1)))[2:]
+ ret = _irregular_nabla(ret,l_vec)
+ ret = _irregular_adj_nabla(ret,l_vec)
+ ret[0:2] = 0
+ ret[-1] = 0
+ ret[2:] *= np.sqrt(l_vec-np.roll(l_vec,1))[2:]
+ return Field(self.domain, val=ret).weight(power=-0.5)
+
+ def _adjoint_times(self, t, spaces):
+ ret = t.copy().weight(power=0.5).val
+ l_vec = self.domain[0].kindex.copy()
+ l_vec[1:] = np.log(l_vec[1:])
+ l_vec[0] = -1.
+ ret[2:] *= np.sqrt(l_vec-np.roll(l_vec,1))[2:]
+ ret[0:2] = 0
+ ret[-1] = 0
+ ret = _irregular_nabla(ret,l_vec)
+ ret = _irregular_adj_nabla(ret,l_vec)
+ # ret[2:] *= np.sqrt(np.log(k_vec)-np.log(np.roll(k_vec,1)))[2:]
+ return Field(self.domain, val=ret).weight(-1)
\ No newline at end of file
diff --git a/nifty/library/operator_library/nonlinear_wiener_filter_curvature.py b/nifty/library/operator_library/nonlinear_wiener_filter_curvature.py
deleted file mode 100644
index f42a35ae899476f284d06f1aa0fc1561cbe1887f..0000000000000000000000000000000000000000
--- a/nifty/library/operator_library/nonlinear_wiener_filter_curvature.py
+++ /dev/null
@@ -1,35 +0,0 @@
-from nifty.operators import EndomorphicOperator,\
- InvertibleOperatorMixin
-
-
-class NonlinearWienerFilterCurvature(InvertibleOperatorMixin, EndomorphicOperator):
- def __init__(self, R, N, S, position, inverter=None, preconditioner=None):
-
- self.R = R
- self.N = N
- self.S = S
- self.position = position
- # if preconditioner is None:
- # preconditioner = self.S.times
- self._domain = self.S.domain
- super(NonlinearWienerFilterCurvature, self).__init__(inverter=inverter,
- preconditioner=preconditioner)
- @property
- def domain(self):
- return self._domain
-
- @property
- def self_adjoint(self):
- return True
-
- @property
- def unitary(self):
- return False
-
- # ---Added properties and methods---
-
- def _times(self, x, spaces):
- return self.R.derived_adjoint_times(
- self.N.inverse_times(self.R.derived_times(
- x, self.position)), self.position)\
- + self.S.inverse_times(x)
diff --git a/nifty/library/operator_library/smoothness_operator.py b/nifty/library/operator_library/smoothness_operator.py
index fca9d91e2fa4fd32a5ff604e201ef096fd77eb1d..24e285a73b582086fbd00108767fe1bf745f44bf 100644
--- a/nifty/library/operator_library/smoothness_operator.py
+++ b/nifty/library/operator_library/smoothness_operator.py
@@ -1,7 +1,7 @@
from nifty import EndomorphicOperator,\
PowerSpace
-from laplace_operator import LaplaceOperator
+from laplace_operator import LogLaplaceOperator
class SmoothnessOperator(EndomorphicOperator):
@@ -21,7 +21,7 @@ class SmoothnessOperator(EndomorphicOperator):
self._sigma = sigma
- self._Laplace = LaplaceOperator(domain=self._domain[0])
+ self._Laplace = LogLaplaceOperator(domain=self._domain[0])
"""
SmoothnessOperator
diff --git a/nifty/minimization/relaxed_newton.py b/nifty/minimization/relaxed_newton.py
index 467029317bba8a780627db6bad8d57858370c2a2..e0e86e778e38bf104833ae78767fd2b7bc4f7589 100644
--- a/nifty/minimization/relaxed_newton.py
+++ b/nifty/minimization/relaxed_newton.py
@@ -37,4 +37,4 @@ class RelaxedNewton(DescentMinimizer):
gradient = energy.gradient
curvature = energy.curvature
descend_direction = curvature.inverse_times(gradient)
- return descend_direction * -1
+ return descend_direction * -1.
diff --git a/nifty/sugar.py b/nifty/sugar.py
index 5f29e9d04e36d8fc4eb0d5f8d95c2939716fedcb..7cd6c1ff18abefb4db0ceb79884b9c4dbe62e437 100644
--- a/nifty/sugar.py
+++ b/nifty/sugar.py
@@ -55,13 +55,13 @@ def generate_posterior_sample(mean, covariance):
mock_signal = power.power_synthesize(real_signal=True)
- noise = N.diagonal().val
+ noise = N.diagonal(bare=True).val
mock_noise = Field.from_random(random_type="normal", domain=N.domain,
std = sqrt(noise), dtype = noise.dtype)
- mock_data = R.derived_times(mock_signal, mean) + mock_noise
+ mock_data = R(mock_signal) + mock_noise
- mock_j = R.derived_adjoint_times(N.inverse_times(mock_data), mean)
+ mock_j = R.adjoint_times(N.inverse_times(mock_data))
mock_m = covariance.inverse_times(mock_j)
sample = mock_signal - mock_m + mean
return sample