Commit 848b1006 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

fix markdown

parent c7a7b50e
import nifty4 as ift
import numpy as np
if __name__ == "__main__":
np.random.seed(42)
# Setting up parameters
correlation_length = 1. # Typical distance over which the field is correlated
field_variance = 2. # Variance of field in position space
response_sigma = 0.02 # Smoothing length of response (in same unit as L)
signal_to_noise = 100 # The signal to noise ratio
np.random.seed(43) # Fixing the random seed
def power_spectrum(k): # Defining the power spectrum
a = 4 * correlation_length * field_variance**2
return a / (1 + k * correlation_length) ** 4
# Setting up the geometry
L = 2. # Total side-length of the domain
N_pixels = 128 # Grid resolution (pixels per axis)
# signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
signal_space = ift.HPSpace(16)
harmonic_space = signal_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, target=signal_space)
power_space = ift.PowerSpace(harmonic_space)
# Creating the mock signal
S = ift.create_power_operator(harmonic_space,
power_spectrum=power_spectrum)
mock_power = ift.PS_field(power_space, power_spectrum)
mock_signal = ift.power_synthesize(mock_power, real_signal=True)
# Setting up an exemplary response
mask = ift.Field.ones(signal_space)
N10 = int(N_pixels/10)
# mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ift.GeometryRemover(signal_space)
R = R*ift.DiagonalOperator(mask)
R = R*HT
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
data_domain = R.target[0]
# Setting up the noise covariance and drawing a random noise realization
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.val.std()/signal_to_noise
N = ift.ScalingOperator(noise_amplitude**2, data_domain)
noise = ift.Field.from_random(
domain=data_domain, random_type='normal',
std=noise_amplitude, mean=0)
data = noiseless_data + noise
# Wiener filter
m0 = ift.Field.zeros(harmonic_space)
ctrl = ift.GradientNormController(tol_abs_gradnorm=0.0001)
ctrl2 = ift.GradientNormController(tol_abs_gradnorm=0.1, name="outer")
inverter = ift.ConjugateGradient(controller=ctrl)
energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R,
N, S, inverter=inverter)
#minimizer = ift.VL_BFGS(controller=ctrl2, max_history_length=20)
minimizer = ift.RelaxedNewton(controller=ctrl2)
#minimizer = ift.SteepestDescent(controller=ctrl2)
me = minimizer(energy)
m = HT(me[0].position)
# Plotting
plotdict = {"colormap": "Planck-like"}
ift.plot(HT(mock_signal), name="mock_signal.png", **plotdict)
logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape)
ift.plot(ift.Field(signal_space, val=ift.dobj.from_global_data(logdata)),
name="log_of_data.png", **plotdict)
ift.plot(m, name='m.png', **plotdict)
# Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober):
pass
proby = Proby(signal_space, probe_count=1)
proby(lambda z: HT(me2[0].curvature.inverse_times(HT.adjoint_times(z))))
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.02)
variance = sm(proby.diagonal.weight(-1))
ift.plot(variance, name='variance.png', **plotdict)
......@@ -131,6 +131,7 @@ typical examples for this category are the :py:class:`ScalingOperator`, which si
Nifty4 allows simple and intuitive construction of combined operators.
As an example, if :math:`A`, :math:`B` and :math:`C` are of type :py:class:`LinearOperator` and :math:`f_1` and :math:`f_2` are fields, writing::
X = A*B.inverse*A.adjoint + C
f2 = X(f1)
......
# 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 <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..operators.inversion_enabler import InversionEnabler
def LogNormalWienerFilterCurvature(R, N, S, ht, expp_sspace, inverter):
LinResp = R * ht.adjoint * expp_sspace * ht
t_op = LinResp.adjoint * N.inverse * LinResp
return InversionEnabler(S.inverse + t_op, inverter)
# 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 <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..minimization.energy import Energy
from ..utilities import memo
from .log_normal_wiener_filter_curvature import LogNormalWienerFilterCurvature
from ..sugar import create_composed_ht_operator
from ..field import exp
class LogNormalWienerFilterEnergy(Energy):
"""The Energy for the log-normal Wiener filter.
It covers the case of linear measurement with
Gaussian noise and Gaussain signal prior with known covariance.
Parameters
----------
position: Field,
The current position.
d: Field,
the data.
R: Operator,
The response operator, description 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, inverter, ht=None):
super(LogNormalWienerFilterEnergy, self).__init__(position=position)
self.d = d
self.R = R
self.N = N
self.S = S
self._inverter = inverter
if ht is None:
self._ht = create_composed_ht_operator(self.S.domain)
else:
self._ht = ht
self._expp_sspace = exp(self._ht(self.position))
Sp = self.S.inverse_times(self.position)
expp = self._ht.adjoint_times(self._expp_sspace)
Rexppd = self.R(expp) - self.d
NRexppd = self.N.inverse_times(Rexppd)
self._value = 0.5*(self.position.vdot(Sp) + Rexppd.vdot(NRexppd))
exppRNRexppd = self._ht.adjoint_times(
self._expp_sspace * self._ht(self.R.adjoint_times(NRexppd)))
self._gradient = Sp + exppRNRexppd
def at(self, position):
return self.__class__(position, self.d, self.R, self.N, self.S,
self._inverter, self._ht)
@property
def value(self):
return self._value
@property
def gradient(self):
return self._gradient
@property
@memo
def curvature(self):
return LogNormalWienerFilterCurvature(
R=self.R, N=self.N, S=self.S, ht=self._ht,
expp_sspace=self._expp_sspace, inverter=self._inverter)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment