Commit 63eb16f4 authored by Martin Reinecke's avatar Martin Reinecke

demo tweaks

parent 7272888e
Pipeline #24245 passed with stage
in 5 minutes and 21 seconds
...@@ -7,7 +7,7 @@ if __name__ == "__main__": ...@@ -7,7 +7,7 @@ if __name__ == "__main__":
correlation_length = 1. # Typical distance over which the field is correlated correlation_length = 1. # Typical distance over which the field is correlated
field_variance = 2. # Variance of field in position space field_variance = 2. # Variance of field in position space
response_sigma = 0.02 # Smoothing length of response (in same unit as L) response_sigma = 0.02 # Smoothing length of response (in same unit as L)
signal_to_noise = 100 # The signal to noise ratio signal_to_noise = 100 # The signal to noise ratio
np.random.seed(43) # Fixing the random seed np.random.seed(43) # Fixing the random seed
def power_spectrum(k): # Defining the power spectrum def power_spectrum(k): # Defining the power spectrum
...@@ -20,66 +20,65 @@ if __name__ == "__main__": ...@@ -20,66 +20,65 @@ if __name__ == "__main__":
# signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels) # signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
signal_space = ift.HPSpace(16) signal_space = ift.HPSpace(16)
harmonic_space = signal_space.get_default_codomain() harmonic_space = signal_space.get_default_codomain()
fft = ift.FFTOperator(harmonic_space, target=signal_space) ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space)
power_space = ift.PowerSpace(harmonic_space) power_space = ift.PowerSpace(harmonic_space)
# Creating the mock signal # Creating the mock signal
S = ift.create_power_operator(harmonic_space, S = ift.create_power_operator(harmonic_space,
power_spectrum=power_spectrum) power_spectrum=power_spectrum)
mock_power = ift.PS_field(power_space, power_spectrum) mock_power = ift.PS_field(power_space, power_spectrum)
mock_signal = fft(ift.power_synthesize(mock_power, real_signal=True)) mock_signal = ift.power_synthesize(mock_power, real_signal=True)
# Setting up an exemplary response # Setting up an exemplary response
mask = ift.Field.ones(signal_space) mask = ift.Field.ones(signal_space)
N10 = int(N_pixels/10) N10 = int(N_pixels/10)
# mask.val[N10*5:N10*9, N10*5:N10*9] = 0. # mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), R = ift.GeometryRemover(signal_space)
exposure=(mask,)) R = R*ift.DiagonalOperator(mask)
R = R*ht
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
#R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
# sensitivity=(mask,))
data_domain = R.target[0] data_domain = R.target[0]
R_harmonic = R*fft
# Setting up the noise covariance and drawing a random noise realization # Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise) noiseless_data = R(mock_signal)
N = ift.DiagonalOperator(ndiag.weight(1)) noise_amplitude = noiseless_data.std()/signal_to_noise
noise = ift.Field.from_random(domain=data_domain, random_type='normal', N = ift.DiagonalOperator(
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) ift.Field.full(data_domain, noise_amplitude**2))
data = R(ift.exp(mock_signal)) + noise noise = ift.Field.from_random(
domain=data_domain, random_type='normal',
std=noise_amplitude, mean=0)
data = noiseless_data + noise
# Wiener filter # Wiener filter
m0 = ift.Field.zeros(harmonic_space) m0 = ift.Field.zeros(harmonic_space)
ctrl = ift.GradientNormController(tol_abs_gradnorm=0.0001) ctrl = ift.GradientNormController(tol_abs_gradnorm=0.0001)
ctrl2 = ift.GradientNormController(tol_abs_gradnorm=0.1, name="outer") ctrl2 = ift.GradientNormController(tol_abs_gradnorm=0.1, name="outer")
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R,
N, S, inverter=inverter) N, S, inverter=inverter)
# minimizer1 = ift.VL_BFGS(controller=ctrl2, max_history_length=20) # minimizer = ift.VL_BFGS(controller=ctrl2, max_history_length=20)
minimizer2 = ift.RelaxedNewton(controller=ctrl2) minimizer = ift.RelaxedNewton(controller=ctrl2)
# minimizer3 = ift.SteepestDescent(controller=ctrl2) # minimizer = ift.SteepestDescent(controller=ctrl2)
# me1 = minimizer1(energy) me = minimizer(energy)
me2 = minimizer2(energy) m = ht(me[0].position)
# me3 = minimizer3(energy)
# m1 = fft(me1[0].position)
m2 = fft(me2[0].position)
# m3 = fft(me3[0].position)
# Plotting # Plotting
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"} "colormap": "Planck-like"}
ift.plot(mock_signal, name="mock_signal.png", **plotdict) ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict)
logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape) 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)), ift.plot(ift.Field(signal_space, val=ift.dobj.from_global_data(logdata)),
name="log_of_data.png", **plotdict) name="log_of_data.png", **plotdict)
# ift.plot(m1,name='m_LBFGS.png', **plotdict) ift.plot(m, name='m.png', **plotdict)
ift.plot(m2, name='m_Newton.png', **plotdict)
# ift.plot(m3, name='m_SteepestDescent.png', **plotdict)
# Probing the variance # Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober): class Proby(ift.DiagonalProberMixin, ift.Prober):
pass pass
proby = Proby(signal_space, probe_count=1) proby = Proby(signal_space, probe_count=1)
proby(lambda z: fft(me2[0].curvature.inverse_times(fft.adjoint_times(z)))) proby(lambda z: ht(me2[0].curvature.inverse_times(ht.adjoint_times(z))))
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.02) sm = ift.FFTSmoothingOperator(signal_space, sigma=0.02)
variance = sm(proby.diagonal.weight(-1)) variance = sm(proby.diagonal.weight(-1))
......
...@@ -2,7 +2,7 @@ import numpy as np ...@@ -2,7 +2,7 @@ import numpy as np
import nifty4 as ift import nifty4 as ift
if __name__ == "__main__": if __name__ == "__main__":
signal_to_noise = 1.5 # The signal to noise ratio signal_to_noise = 0.5 # The signal to noise ratio
# Setting up parameters # Setting up parameters
L_1 = 2. # Total side-length of the domain L_1 = 2. # Total side-length of the domain
...@@ -30,7 +30,7 @@ if __name__ == "__main__": ...@@ -30,7 +30,7 @@ if __name__ == "__main__":
harmonic_domain = ift.DomainTuple.make((harmonic_space_1, harmonic_domain = ift.DomainTuple.make((harmonic_space_1,
harmonic_space_2)) harmonic_space_2))
fft_1 = ift.FFTOperator(harmonic_domain, space=0) ht_1 = ift.HarmonicTransformOperator(harmonic_domain, space=0)
power_space_1 = ift.PowerSpace(harmonic_space_1) power_space_1 = ift.PowerSpace(harmonic_space_1)
mock_power_1 = ift.PS_field(power_space_1, power_spectrum_1) mock_power_1 = ift.PS_field(power_space_1, power_spectrum_1)
...@@ -39,12 +39,12 @@ if __name__ == "__main__": ...@@ -39,12 +39,12 @@ if __name__ == "__main__":
a = 4 * correlation_length_2 * field_variance_2**2 a = 4 * correlation_length_2 * field_variance_2**2
return a / (1 + k * correlation_length_2) ** 2.5 return a / (1 + k * correlation_length_2) ** 2.5
fft_2 = ift.FFTOperator(mid_domain, space=1) ht_2 = ift.HarmonicTransformOperator(mid_domain, space=1)
power_space_2 = ift.PowerSpace(harmonic_space_2) power_space_2 = ift.PowerSpace(harmonic_space_2)
mock_power_2 = ift.PS_field(power_space_2, power_spectrum_2) mock_power_2 = ift.PS_field(power_space_2, power_spectrum_2)
fft = fft_2*fft_1 ht = ht_2*ht_1
mock_power = ift.Field( mock_power = ift.Field(
(power_space_1, power_space_2), (power_space_1, power_space_2),
...@@ -57,7 +57,7 @@ if __name__ == "__main__": ...@@ -57,7 +57,7 @@ if __name__ == "__main__":
S = ift.DiagonalOperator(diagonal) S = ift.DiagonalOperator(diagonal)
np.random.seed(10) np.random.seed(10)
mock_signal = fft(ift.power_synthesize(mock_power, real_signal=True)) mock_signal = ift.power_synthesize(mock_power, real_signal=True)
# Setting up a exemplary response # Setting up a exemplary response
N1_10 = int(N_pixels_1/10) N1_10 = int(N_pixels_1/10)
...@@ -70,49 +70,50 @@ if __name__ == "__main__": ...@@ -70,49 +70,50 @@ if __name__ == "__main__":
mask_2[N2_10*7:N2_10*9] = 0. mask_2[N2_10*7:N2_10*9] = 0.
mask_2 = ift.Field(signal_space_2, ift.dobj.from_global_data(mask_2)) mask_2 = ift.Field(signal_space_2, ift.dobj.from_global_data(mask_2))
R = ift.ResponseOperator(signal_domain, R = ift.GeometryRemover(signal_domain)
sigma=(response_sigma_1, response_sigma_2), R = R*ift.DiagonalOperator(mask_1, signal_domain,spaces=0)
exposure=(mask_1, mask_2)) R = R*ift.DiagonalOperator(mask_2, signal_domain,spaces=1)
R = R*ht
R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 0,
response_sigma_1)
R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 1,
response_sigma_2)
data_domain = R.target data_domain = R.target
R_harmonic = R*fft
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.std()/signal_to_noise
# Setting up the noise covariance and drawing a random noise realization # Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise) ndiag = ift.Field.full(data_domain, noise_amplitude**2)
N = ift.DiagonalOperator(ndiag.weight(1)) N = ift.DiagonalOperator(ndiag)
noise = ift.Field.from_random( noise = ift.Field.from_random(
domain=data_domain, random_type='normal', domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) std=noise_amplitude, mean=0)
data = R(mock_signal) + noise data = noiseless_data + noise
# Wiener filter # Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data)) j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1) ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature( wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R_harmonic, inverter=inverter) S=S, N=N, R=R, inverter=inverter)
m_k = wiener_curvature.inverse_times(j) m_k = wiener_curvature.inverse_times(j)
m = fft(m_k) m = ht(m_k)
# Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober):
pass
proby = Proby((signal_space_1, signal_space_2), probe_count=1, ncpu=1)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
# sm = SmoothingOperator(signal_space, sigma=0.02)
# variance = sm(proby.diagonal.weight(-1))
variance = proby.diagonal.weight(-1)
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
sm = ift.FFTSmoothingOperator(plot_space, sigma=0.03)
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"} "colormap": "Planck-like"}
ift.plot( plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))), ift.plot(ift.Field(plot_space,val=ht(mock_signal).val), name='mock_signal.png',
name='uncertainty.png', zmin=0., zmax=3., title="Uncertainty map", **plotdict)
**plotdict) ift.plot(ift.Field(plot_space,val=data.val), name='data.png', **plotdict)
ift.plot(ift.Field(plot_space, val=mock_signal.val.real), ift.plot(ift.Field(plot_space,val=m.val), name='map.png', **plotdict)
name='mock_signal.png', **plotdict) # sampling the uncertainty map
ift.plot(ift.Field(plot_space, val=data.val.real), sample_variance = ift.Field.zeros(signal_domain)
name='data.png', **plotdict) sample_mean = ift.Field.zeros(signal_domain)
ift.plot(ift.Field(plot_space, val=m.val.real), name='map.png', **plotdict) n_samples = 10
for i in range(n_samples):
sample = ht(wiener_curvature.generate_posterior_sample()) + m
sample_variance += sample**2
sample_mean += sample
variance = sample_variance/n_samples - (sample_mean/n_samples)**2
ift.plot(ift.Field(plot_space, val=ift.sqrt(variance).val), name="uncertainty.png", **plotdict)
...@@ -6,7 +6,7 @@ if __name__ == "__main__": ...@@ -6,7 +6,7 @@ if __name__ == "__main__":
# Setting up parameters # Setting up parameters
L = 2. # Total side-length of the domain L = 2. # Total side-length of the domain
N_pixels = 512 # Grid resolution (pixels per axis) N_pixels = 512 # Grid resolution (pixels per axis)
correlation_length_scale = 1. # Typical distance over which the field is correlation_length_scale = .2 # Typical distance over which the field is
# correlated # correlated
fluctuation_scale = 2. # Variance of field in position space fluctuation_scale = 2. # Variance of field in position space
response_sigma = 0.05 # Smoothing length of response response_sigma = 0.05 # Smoothing length of response
...@@ -19,7 +19,7 @@ if __name__ == "__main__": ...@@ -19,7 +19,7 @@ if __name__ == "__main__":
signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels) signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = signal_space.get_default_codomain() harmonic_space = signal_space.get_default_codomain()
fft = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space)
power_space = ift.PowerSpace( power_space = ift.PowerSpace(
harmonic_space, binbounds=ift.PowerSpace.useful_binbounds( harmonic_space, binbounds=ift.PowerSpace.useful_binbounds(
harmonic_space, logarithmic=True)) harmonic_space, logarithmic=True))
...@@ -37,47 +37,43 @@ if __name__ == "__main__": ...@@ -37,47 +37,43 @@ if __name__ == "__main__":
mask = ift.Field(signal_space, ift.dobj.from_global_data(mask)) mask = ift.Field(signal_space, ift.dobj.from_global_data(mask))
R = ift.GeometryRemover(signal_space) R = ift.GeometryRemover(signal_space)
R = R*ift.DiagonalOperator(mask) R = R*ift.DiagonalOperator(mask)
R = R*fft R = R*ht
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma) R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
data_domain = R.target[0] data_domain = R.target[0]
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.std()/signal_to_noise
# Setting up the noise covariance and drawing a random noise realization # Setting up the noise covariance and drawing a random noise realization
ndiag = 1e-8*ift.Field.full(data_domain, fft(mock_signal).var()/signal_to_noise) ndiag = ift.Field.full(data_domain, noise_amplitude**2)
N = ift.DiagonalOperator(ndiag) N = ift.DiagonalOperator(ndiag)
noise = ift.Field.from_random( noise = ift.Field.from_random(
domain=data_domain, random_type='normal', domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) std=noise_amplitude, mean=0)
data = R(mock_signal) + noise data = noiseless_data + noise
# Wiener filter # Wiener filter
j = R.adjoint_times(N.inverse_times(data)) j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-6) ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-2)
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature( wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, inverter=inverter) S=S, N=N, R=R, inverter=inverter)
m_k = wiener_curvature.inverse_times(j) m_k = wiener_curvature.inverse_times(j)
m = fft(m_k) m = ht(m_k)
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"} "colormap": "Planck-like"}
ift.plot(mock_signal, name="mock_signal.png", **plotdict) ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict)
ift.plot(ift.Field(signal_space, val=data.val), ift.plot(ift.Field(signal_space, val=data.val),
name="data.png", **plotdict) name="data.png", **plotdict)
ift.plot(m, name="map.png", **plotdict) ift.plot(m, name="map.png", **plotdict)
# Probing the uncertainty
class Proby(ift.DiagonalProberMixin, ift.Prober):
pass
proby = Proby(harmonic_space, probe_count=1, ncpu=1)
proby(lambda z: wiener_curvature.inverse_times(z))
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03) # sampling the uncertainty map
variance = ift.sqrt(sm(proby.diagonal.weight(-1))) sample_variance = ift.Field.zeros(signal_space)
sample_mean = ift.Field.zeros(signal_space)
# Plotting n_samples = 10
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", for i in range(n_samples):
"colormap": "Planck-like"} sample = ht(wiener_curvature.generate_posterior_sample()) + m
ift.plot(variance, name="uncertainty.png", **plotdict) sample_variance += sample**2
ift.plot(mock_signal, name="mock_signal.png", **plotdict) sample_mean += sample
ift.plot(ift.Field(signal_space, val=data.val), variance = sample_variance/n_samples - (sample_mean/n_samples)**2
name="data.png", **plotdict) ift.plot(ift.sqrt(variance), name="uncertainty.png", **plotdict)
ift.plot(m, name="map.png", **plotdict)
...@@ -27,8 +27,8 @@ if __name__ == "__main__": ...@@ -27,8 +27,8 @@ if __name__ == "__main__":
# Set up the geometry # Set up the geometry
s_space = ift.RGSpace([N_pixels, N_pixels], distances=pixel_width) s_space = ift.RGSpace([N_pixels, N_pixels], distances=pixel_width)
fft = ift.FFTOperator(s_space) h_space = s_space.get_default_codomain()
h_space = fft.target[0] ht = ift.HarmonicTransformOperator(h_space, s_space)
p_space = ift.PowerSpace(h_space) p_space = ift.PowerSpace(h_space)
# Create mock data # Create mock data
...@@ -37,19 +37,19 @@ if __name__ == "__main__": ...@@ -37,19 +37,19 @@ if __name__ == "__main__":
sp = ift.PS_field(p_space, pow_spec) sp = ift.PS_field(p_space, pow_spec)
sh = ift.power_synthesize(sp, real_signal=True) sh = ift.power_synthesize(sp, real_signal=True)
ss = fft.inverse_times(sh)
R = ift.FFTSmoothingOperator(s_space, sigma=response_sigma) R = ht*ift.create_harmonic_smoothing_operator((h_space,),0,response_sigma)
signal_to_noise = 1 noiseless_data = R(sh)
diag = ift.Field(s_space, ss.var()/signal_to_noise).weight(1) signal_to_noise = 1.
N = ift.DiagonalOperator(diag) noise_amplitude = noiseless_data.std()/signal_to_noise
N = ift.DiagonalOperator(ift.Field.full(s_space, noise_amplitude**2))
n = ift.Field.from_random(domain=s_space, n = ift.Field.from_random(domain=s_space,
random_type='normal', random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise), std=noise_amplitude,
mean=0) mean=0)
d = R(ss) + n d = noiseless_data + n
# Wiener filter # Wiener filter
...@@ -57,8 +57,14 @@ if __name__ == "__main__": ...@@ -57,8 +57,14 @@ if __name__ == "__main__":
IC = ift.GradientNormController(name="inverter", iteration_limit=500, IC = ift.GradientNormController(name="inverter", iteration_limit=500,
tol_abs_gradnorm=0.1) tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=IC) inverter = ift.ConjugateGradient(controller=IC)
S_inv = fft.adjoint*Sh.inverse*fft D = (R.adjoint*N.inverse*R + Sh.inverse).inverse
D = (R.adjoint*N.inverse*R + S_inv).inverse
# MR FIXME: we can/should provide a preconditioner here as well! # MR FIXME: we can/should provide a preconditioner here as well!
D = ift.InversionEnabler(D, inverter) D = ift.InversionEnabler(D, inverter)
m = D(j) m = D(j)
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
ift.plot(ht(sh), name="mock_signal.png", **plotdict)
ift.plot(ift.Field(s_space, val=d.val),
name="data.png", **plotdict)
ift.plot(ht(m), name="map.png", **plotdict)
...@@ -16,14 +16,14 @@ if __name__ == "__main__": ...@@ -16,14 +16,14 @@ if __name__ == "__main__":
# sigma of field in position space sqrt(<|s_x|^2>) # sigma of field in position space sqrt(<|s_x|^2>)
field_sigma = 2. * nu.K field_sigma = 2. * nu.K
# smoothing length of response # smoothing length of response
response_sigma = 0.01*nu.m response_sigma = 0.03*nu.m
# The signal to noise ratio ***CURRENTLY BROKEN*** # The signal to noise ratio
signal_to_noise = 70.7 signal_to_noise = 1
# note that field_variance**2 = a*k_0/4. for this analytic form of power # note that field_variance**2 = a*k_0/4. for this analytic form of power
# spectrum # spectrum
def power_spectrum(k): def power_spectrum(k):
#RL FIXME: signal_amplitude is not how much signal varies #RL FIXME: signal_amplitude is not how much signal varies
cldim = correlation_length**(2*dimensionality) cldim = correlation_length**(2*dimensionality)
a = 4/(2*np.pi) * cldim * field_sigma**2 a = 4/(2*np.pi) * cldim * field_sigma**2
# to be integrated over spherical shells later on # to be integrated over spherical shells later on
...@@ -39,7 +39,7 @@ if __name__ == "__main__": ...@@ -39,7 +39,7 @@ if __name__ == "__main__":
signal_space = ift.RGSpace(shape, distances=L/N_pixels) signal_space = ift.RGSpace(shape, distances=L/N_pixels)
harmonic_space = signal_space.get_default_codomain() harmonic_space = signal_space.get_default_codomain()
fft = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space)
power_space = ift.PowerSpace(harmonic_space) power_space = ift.PowerSpace(harmonic_space)
# Creating the mock data # Creating the mock data
...@@ -53,18 +53,18 @@ if __name__ == "__main__": ...@@ -53,18 +53,18 @@ if __name__ == "__main__":
sensitivity = (1./nu.m)**dimensionality/nu.K sensitivity = (1./nu.m)**dimensionality/nu.K
R = ift.GeometryRemover(signal_space) R = ift.GeometryRemover(signal_space)
R = R*ift.ScalingOperator(sensitivity, signal_space) R = R*ift.ScalingOperator(sensitivity, signal_space)
R = R*fft R = R*ht
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma) R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
data_domain = R.target[0] data_domain = R.target[0]
noise_amplitude = 1./signal_to_noise*field_sigma*sensitivity*((L/N_pixels)**dimensionality) noiseless_data = R(mock_signal)
print "noise amplitude:", noise_amplitude noise_amplitude = noiseless_data.std()/signal_to_noise
N = ift.DiagonalOperator( N = ift.DiagonalOperator(
ift.Field.full(data_domain, noise_amplitude**2)) ift.Field.full(data_domain, noise_amplitude**2))
noise = ift.Field.from_random( noise = ift.Field.from_random(
domain=data_domain, random_type='normal', domain=data_domain, random_type='normal',
std=noise_amplitude, mean=0) std=noise_amplitude, mean=0)
data = R(mock_signal) + noise data = noiseless_data + noise
# Wiener filter # Wiener filter
j = R.adjoint_times(N.inverse_times(data)) j = R.adjoint_times(N.inverse_times(data))
...@@ -75,14 +75,14 @@ if __name__ == "__main__": ...@@ -75,14 +75,14 @@ if __name__ == "__main__":
S=S, N=N, R=R, inverter=inverter) S=S, N=N, R=R, inverter=inverter)
m = wiener_curvature.inverse_times(j) m = wiener_curvature.inverse_times(j)
m_s = fft(m) m_s = ht(m)
sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m) sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m)
ift.plot(ift.Field(sspace2, fft(mock_signal).val)/nu.K, name="mock_signal.png") ift.plot(ift.Field(sspace2, ht(mock_signal).val)/nu.K, name="mock_signal.png")
#data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape) #data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape)
#data = ift.Field(sspace2, val=ift.dobj.from_global_data(data)) #data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))
ift.plot(ift.Field(sspace2, val=fft(R.adjoint_times(data)).val), name="data.png") ift.plot(ift.Field(sspace2, val=data.val), name="data.png")
print "msig",np.min(fft(mock_signal).val)/nu.K, np.max(fft(mock_signal).val)/nu.K print "msig",np.min(ht(mock_signal).val)/nu.K, np.max(ht(mock_signal).val)/nu.K
print "map",np.min(m_s.val)/nu.K, np.max(m_s.val)/nu.K print "map",np.min(m_s.val)/nu.K, np.max(m_s.val)/nu.K
ift.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png") ift.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png")
...@@ -11,7 +11,7 @@ if __name__ == "__main__": ...@@ -11,7 +11,7 @@ if __name__ == "__main__":
# Define associated harmonic space and harmonic transformation # Define associated harmonic space and harmonic transformation
h_space = s_space.get_default_codomain() h_space = s_space.get_default_codomain()
fft = ift.HarmonicTransformOperator(h_space, s_space) ht = ift.HarmonicTransformOperator(h_space, s_space)
# Setting up power space # Setting up power space
p_space = ift.PowerSpace(h_space) p_space = ift.PowerSpace(h_space)
...@@ -25,27 +25,26 @@ if __name__ == "__main__": ...@@ -25,27 +25,26 @@ if __name__ == "__main__":
# Drawing a sample sh from the prior distribution in harmonic space # Drawing a sample sh from the prior distribution in harmonic space
sp = ift.PS_field(p_space, p_spec) sp = ift.PS_field(p_space, p_spec)
sh = ift.power_synthesize(sp, real_signal=True) sh = ift.power_synthesize(sp, real_signal=True)
ss = fft.times(sh)
# Choosing the measurement instrument # Choosing the measurement instrument
# Instrument = ift.FFTSmoothingOperator(s_space, sigma=0.05)
diag = np.ones(s_space.shape) diag = np.ones(s_space.shape)
diag[20