Commit 0a688d2b authored by Martin Reinecke's avatar Martin Reinecke
Browse files

extensive cleanups

parent d0c9a8ab
Pipeline #21555 passed with stage
in 4 minutes and 21 seconds
...@@ -106,7 +106,9 @@ if __name__ == "__main__": ...@@ -106,7 +106,9 @@ if __name__ == "__main__":
def ps0(k): def ps0(k):
return (1./(1.+k)**2) return (1./(1.+k)**2)
t0 = ift.Field(p_space, val=np.log(1./(1+p_space.k_lengths)**2))
t0 = ift.Field(p_space,
val=ift.dobj.from_global_data(np.log(1./(1+p_space.k_lengths)**2)))
for i in range(500): for i in range(500):
S0 = ift.create_power_operator(h_space, power_spectrum=ps0) S0 = ift.create_power_operator(h_space, power_spectrum=ps0)
...@@ -138,6 +140,6 @@ if __name__ == "__main__": ...@@ -138,6 +140,6 @@ if __name__ == "__main__":
t0 = power_energy.position.real t0 = power_energy.position.real
# Plot current estimate # Plot current estimate
print(i) ift.dobj.mprint(i)
if i % 5 == 0: if i % 5 == 0:
plot_parameters(m0, t0, ift.log(sp), data_power) plot_parameters(m0, t0, ift.log(sp), data_power)
...@@ -3,35 +3,38 @@ import numpy as np ...@@ -3,35 +3,38 @@ import numpy as np
if __name__ == "__main__": if __name__ == "__main__":
np.random.seed(42) np.random.seed(42)
# Setting up parameters |\label{code:wf_parameters}| # Setting up parameters
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
a = 4 * correlation_length * field_variance**2 a = 4 * correlation_length * field_variance**2
return a / (1 + k * correlation_length) ** 4 return a / (1 + k * correlation_length) ** 4
# Setting up the geometry |\label{code:wf_geometry}| # Setting up the geometry
L = 2. # Total side-length of the domain L = 2. # Total side-length of the domain
N_pixels = 128 # Grid resolution (pixels per axis) N_pixels = 128 # Grid resolution (pixels per axis)
#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) fft = ift.FFTOperator(harmonic_space, target=signal_space)
power_space = ift.PowerSpace(harmonic_space) power_space = ift.PowerSpace(harmonic_space)
# Creating the mock signal |\label{code:wf_mock_signal}| # Creating the mock signal
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum) S = ift.create_power_operator(harmonic_space,
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 = fft(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,), exposure=(mask,)) #|\label{code:wf_response}| R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
exposure=(mask,))
data_domain = R.target[0] data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R]) R_harmonic = ift.ComposedOperator([fft, R])
...@@ -39,40 +42,52 @@ if __name__ == "__main__": ...@@ -39,40 +42,52 @@ if __name__ == "__main__":
ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise) ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag.weight(1)) N = ift.DiagonalOperator(ndiag.weight(1))
noise = ift.Field.from_random(domain=data_domain, random_type='normal', noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(ift.exp(mock_signal)) + noise #|\label{code:wf_mock_data}| data = R(ift.exp(mock_signal)) + noise
# Wiener filter # Wiener filter
m0 = ift.Field.zeros(harmonic_space) m0 = ift.Field.zeros(harmonic_space)
ctrl = ift.GradientNormController(verbose=False,tol_abs_gradnorm=1) ctrl = ift.GradientNormController(verbose=False, tol_abs_gradnorm=1)
ctrl2 = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1, name="outer") ctrl2 = ift.GradientNormController(verbose=True, 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, N, S, inverter=inverter) energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic,
minimizer1 = ift.VL_BFGS(controller=ctrl2,max_history_length=20) N, S, inverter=inverter)
# minimizer1 = ift.VL_BFGS(controller=ctrl2, max_history_length=20)
minimizer2 = ift.RelaxedNewton(controller=ctrl2) minimizer2 = ift.RelaxedNewton(controller=ctrl2)
minimizer3 = ift.SteepestDescent(controller=ctrl2) # minimizer3 = ift.SteepestDescent(controller=ctrl2)
#me1 = minimizer1(energy) # me1 = minimizer1(energy)
me2 = minimizer2(energy) me2 = minimizer2(energy)
#me3 = minimizer3(energy) # me3 = minimizer3(energy)
#m1 = fft(me1[0].position) # m1 = fft(me1[0].position)
m2 = fft(me2[0].position) m2 = fft(me2[0].position)
#m3 = fft(me3[0].position) # m3 = fft(me3[0].position)
#Plotting #|\label{code:wf_plotting}| # Plotting
ift.plotting.plot(mock_signal.real,name='mock_signal.pdf', colormap="plasma",xlabel="Pixel Index",ylabel="Pixel Index") ift.plotting.plot(mock_signal.real, name='mock_signal.pdf',
ift.plotting.plot(ift.Field(signal_space, val=np.log(data.val.real).reshape(signal_space.shape)),name="log_of_data.pdf", colormap="plasma",xlabel="Pixel Index",ylabel="Pixel Index") colormap="plasma", xlabel="Pixel Index",
#ift.plotting.plot(m1.real,name='m_LBFGS.pdf', colormap="plasma",xlabel="Pixel Index",ylabel="Pixel Index") ylabel="Pixel Index")
ift.plotting.plot(m2.real,name='m_Newton.pdf', colormap="plasma",xlabel="Pixel Index",ylabel="Pixel Index") logdata = np.log(ift.dobj.to_global_data(data.val.real)).reshape(signal_space.shape)
#ift.plotting.plot(m3.real,name='m_SteepestDescent.pdf', colormap="plasma",xlabel="Pixel Index",ylabel="Pixel Index") ift.plotting.plot(ift.Field(signal_space,
val=ift.dobj.from_global_data(logdata)),
name="log_of_data.pdf", colormap="plasma",
xlabel="Pixel Index", ylabel="Pixel Index")
# ift.plotting.plot(m1.real,name='m_LBFGS.pdf', colormap="plasma",
# xlabel="Pixel Index", ylabel="Pixel Index")
ift.plotting.plot(m2.real, name='m_Newton.pdf', colormap="plasma",
xlabel="Pixel Index", ylabel="Pixel Index")
# ift.plotting.plot(m3.real, name='m_SteepestDescent.pdf',
# colormap="plasma", xlabel="Pixel Index",
# ylabel="Pixel Index")
# Probing the variance # Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober): pass class Proby(ift.DiagonalProberMixin, ift.Prober):
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: fft(me2[0].curvature.inverse_times(fft.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))
ift.plotting.plot(variance, name = 'variance.pdf') ift.plotting.plot(variance, name='variance.pdf')
...@@ -2,82 +2,82 @@ import numpy as np ...@@ -2,82 +2,82 @@ import numpy as np
import nifty2go as ift import nifty2go as ift
if __name__ == "__main__": if __name__ == "__main__":
signal_to_noise = 1.5 # The signal to noise ratio signal_to_noise = 1.5 # The signal to noise ratio
# Setting up parameters |\label{code:wf_parameters}| # Setting up parameters
correlation_length_1 = 1. # Typical distance over which the field is correlated correlation_length_1 = 1. # Typical distance over which the field is correlated
field_variance_1 = 2. # Variance of field in position space field_variance_1 = 2. # Variance of field in position space
response_sigma_1 = 0.05 # Smoothing length of response (in same unit as L) response_sigma_1 = 0.05 # Smoothing length of response (in same unit as L)
def power_spectrum_1(k): # note: field_variance**2 = a*k_0/4. def power_spectrum_1(k): # note: field_variance**2 = a*k_0/4.
a = 4 * correlation_length_1 * field_variance_1**2 a = 4 * correlation_length_1 * field_variance_1**2
return a / (1 + k * correlation_length_1) ** 4. return a / (1 + k * correlation_length_1) ** 4.
# Setting up the geometry |\label{code:wf_geometry}| # Setting up the geometry
L_1 = 2. # Total side-length of the domain L_1 = 2. # Total side-length of the domain
N_pixels_1 = 512 # Grid resolution (pixels per axis) N_pixels_1 = 512 # Grid resolution (pixels per axis)
signal_space_1 = ift.RGSpace([N_pixels_1], distances=L_1/N_pixels_1) signal_space_1 = ift.RGSpace([N_pixels_1], distances=L_1/N_pixels_1)
harmonic_space_1 = signal_space_1.get_default_codomain() harmonic_space_1 = signal_space_1.get_default_codomain()
# Setting up the geometry |\label{code:wf_geometry}| # Setting up the geometry
L_2 = 2. # Total side-length of the domain L_2 = 2. # Total side-length of the domain
N_pixels_2 = 512 # Grid resolution (pixels per axis) N_pixels_2 = 512 # Grid resolution (pixels per axis)
signal_space_2 = ift.RGSpace([N_pixels_2], distances=L_2/N_pixels_2) signal_space_2 = ift.RGSpace([N_pixels_2], distances=L_2/N_pixels_2)
harmonic_space_2 = signal_space_2.get_default_codomain() harmonic_space_2 = signal_space_2.get_default_codomain()
signal_domain = ift.DomainTuple.make((signal_space_1, signal_space_2)) signal_domain = ift.DomainTuple.make((signal_space_1, signal_space_2))
mid_domain = ift.DomainTuple.make((signal_space_1, harmonic_space_2)) mid_domain = ift.DomainTuple.make((signal_space_1, harmonic_space_2))
harmonic_domain = ift.DomainTuple.make((harmonic_space_1, harmonic_space_2)) harmonic_domain = ift.DomainTuple.make((harmonic_space_1,
harmonic_space_2))
fft_1 = ift.FFTOperator(harmonic_domain, space=0) fft_1 = ift.FFTOperator(harmonic_domain, space=0)
power_space_1 = ift.PowerSpace(harmonic_space_1) power_space_1 = ift.PowerSpace(harmonic_space_1)
mock_power_1 = ift.Field(power_space_1, val=power_spectrum_1(power_space_1.k_lengths)) mock_power_1 = ift.PS_field(power_space_1, power_spectrum_1)
# Setting up parameters
correlation_length_2 = 1. # Typical distance over which the field is correlated
field_variance_2 = 2. # Variance of field in position space
response_sigma_2 = 0.01 # Smoothing length of response (in same unit as L)
# Setting up parameters |\label{code:wf_parameters}| def power_spectrum_2(k): # note: field_variance**2 = a*k_0/4.
correlation_length_2 = 1. # Typical distance over which the field is correlated
field_variance_2 = 2. # Variance of field in position space
response_sigma_2 = 0.01 # Smoothing length of response (in same unit as L)
def power_spectrum_2(k): # note: field_variance**2 = a*k_0/4.
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) fft_2 = ift.FFTOperator(mid_domain, space=1)
power_space_2 = ift.PowerSpace(harmonic_space_2) power_space_2 = ift.PowerSpace(harmonic_space_2)
mock_power_2 = ift.Field(power_space_2, val=power_spectrum_2(power_space_2.k_lengths)) mock_power_2 = ift.PS_field(power_space_2, power_spectrum_2)
fft = ift.ComposedOperator((fft_1, fft_2)) fft = ift.ComposedOperator((fft_1, fft_2))
mock_power = ift.Field(domain=(power_space_1, power_space_2), mock_power = ift.Field(domain=(power_space_1, power_space_2),
val=np.outer(mock_power_1.val, mock_power_2.val)) val=ift.dobj.from_global_data(np.outer(ift.dobj.to_global_data(mock_power_1.val), ift.dobj.to_global_data(mock_power_2.val))))
diagonal = ift.power_synthesize_special(mock_power, spaces=(0, 1))**2 diagonal = ift.power_synthesize_special(mock_power, spaces=(0, 1))**2
diagonal = diagonal.real diagonal = diagonal.real
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 = fft(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)
mask_1 = ift.Field.ones(signal_space_1) mask_1 = np.ones(signal_space_1.shape)
mask_1.val[N1_10*7:N1_10*9] = 0. mask_1[N1_10*7:N1_10*9] = 0.
mask_1 = ift.Field(signal_space_1, ift.dobj.from_global_data(mask_1))
N2_10 = int(N_pixels_2/10) N2_10 = int(N_pixels_2/10)
mask_2 = ift.Field.ones(signal_space_2) mask_2 = np.ones(signal_space_2.shape)
mask_2.val[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))
R = ift.ResponseOperator(signal_domain,spaces=(0,1), R = ift.ResponseOperator(signal_domain, spaces=(0, 1),
sigma=(response_sigma_1, response_sigma_2), sigma=(response_sigma_1, response_sigma_2),
exposure=(mask_1, mask_2)) #|\label{code:wf_response}| exposure=(mask_1, mask_2))
data_domain = R.target data_domain = R.target
R_harmonic = ift.ComposedOperator([fft, R]) R_harmonic = ift.ComposedOperator([fft, R])
...@@ -87,21 +87,23 @@ if __name__ == "__main__": ...@@ -87,21 +87,23 @@ if __name__ == "__main__":
noise = ift.Field.from_random(domain=data_domain, random_type='normal', noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), std=mock_signal.std()/np.sqrt(signal_to_noise),
mean=0) mean=0)
data = R(mock_signal) + noise #|\label{code:wf_mock_data}| data = R(mock_signal) + noise
# Wiener filter # Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data)) j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1) ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic) wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N,
R=R_harmonic)
wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter) wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter)
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}| m_k = wiener_curvature.inverse_times(j)
m = fft(m_k) m = fft(m_k)
# Probing the variance # Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober): pass class Proby(ift.DiagonalProberMixin, ift.Prober):
proby = Proby((signal_space_1, signal_space_2), probe_count=1,ncpu=1) 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)))) proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
# sm = SmoothingOperator(signal_space, sigma=0.02) # sm = SmoothingOperator(signal_space, sigma=0.02)
# variance = sm(proby.diagonal.weight(-1)) # variance = sm(proby.diagonal.weight(-1))
......
...@@ -3,34 +3,38 @@ import numpy as np ...@@ -3,34 +3,38 @@ import numpy as np
if __name__ == "__main__": if __name__ == "__main__":
# Setting up parameters |\label{code:wf_parameters}| # Setting up parameters
correlation_length_scale = 1. # Typical distance over which the field is correlated correlation_length_scale = 1. # Typical distance over which the field is 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 (in same unit as L) response_sigma = 0.05 # Smoothing length of response (in same unit as L)
signal_to_noise = 1.5 # The signal to noise ratio signal_to_noise = 1.5 # 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
a = 4 * correlation_length_scale * fluctuation_scale**2 a = 4 * correlation_length_scale * fluctuation_scale**2
return a / (1 + (k * correlation_length_scale)**2) ** 2 return a / (1 + (k * correlation_length_scale)**2) ** 2
# Setting up the geometry |\label{code:wf_geometry}| # Setting up the geometry
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)
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.FFTOperator(harmonic_space, target=signal_space) fft = ift.FFTOperator(harmonic_space, target=signal_space)
power_space = ift.PowerSpace(harmonic_space,binbounds=ift.PowerSpace.useful_binbounds(harmonic_space,logarithmic=True)) power_space = ift.PowerSpace(harmonic_space, binbounds=ift.PowerSpace.useful_binbounds(harmonic_space,logarithmic=True))
# Creating the mock signal |\label{code:wf_mock_signal}| # Creating the mock signal
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum) S = ift.create_power_operator(harmonic_space,
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 = fft(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 = np.ones(signal_space.shape)
N10 = int(N_pixels/10) N10 = int(N_pixels/10)
mask.val[N10*5:N10*9, N10*5:N10*9] = 0. mask[N10*5:N10*9, N10*5:N10*9] = 0.
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}| mask = ift.Field(signal_space, ift.dobj.from_global_data(mask))
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
exposure=(mask,))
data_domain = R.target[0] data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R]) R_harmonic = ift.ComposedOperator([fft, R])
...@@ -39,28 +43,29 @@ if __name__ == "__main__": ...@@ -39,28 +43,29 @@ if __name__ == "__main__":
N = ift.DiagonalOperator(ndiag.weight(1)) N = ift.DiagonalOperator(ndiag.weight(1))
noise = ift.Field.from_random(domain=data_domain, random_type='normal', noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise #|\label{code:wf_mock_data}| data = R(mock_signal) + noise
# Wiener filter # Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data)) j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1) ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic) wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N,
wiener_curvature =ift.InversionEnabler(wiener_curvature, inverter) R=R_harmonic)
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}| wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter)
m_k = wiener_curvature.inverse_times(j)
m = fft(m_k) m = fft(m_k)
# Probing the uncertainty |\label{code:wf_uncertainty_probing}| # Probing the uncertainty
class Proby(ift.DiagonalProberMixin, ift.Prober): pass class Proby(ift.DiagonalProberMixin, ift.Prober):
proby = Proby(signal_space, probe_count=1,ncpu=1) pass
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z)))) #|\label{code:wf_variance_fft_wrap}| proby = Proby(signal_space, probe_count=1, ncpu=1)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03) sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03)
variance = ift.sqrt(sm(proby.diagonal.weight(-1))) #|\label{code:wf_variance_weighting}| variance = ift.sqrt(sm(proby.diagonal.weight(-1)))
# Plotting #|\label{code:wf_plotting}| # Plotting
ift.plotting.plot(variance,name="uncertainty.pdf",xlabel='Pixel index', ylabel='Pixel index') ift.plotting.plot(variance,name="uncertainty.pdf",xlabel='Pixel index', ylabel='Pixel index')
ift.plotting.plot(mock_signal,name="mock_signal.pdf",xlabel='Pixel index', ylabel='Pixel index') ift.plotting.plot(mock_signal,name="mock_signal.pdf",xlabel='Pixel index', ylabel='Pixel index')
ift.plotting.plot(ift.Field(signal_space, val=data.val),name="data.pdf",xlabel='Pixel index', ylabel='Pixel index') ift.plotting.plot(ift.Field(signal_space, val=data.val),name="data.pdf",xlabel='Pixel index', ylabel='Pixel index')
ift.plotting.plot(m,name="map.pdf",xlabel='Pixel index', ylabel='Pixel index') ift.plotting.plot(m,name="map.pdf",xlabel='Pixel index', ylabel='Pixel index')
from __future__ import print_function
import nifty2go as ift import nifty2go as ift
import numpy as np import numpy as np
...@@ -20,9 +19,9 @@ diagOp = ift.DiagonalOperator(f) ...@@ -20,9 +19,9 @@ diagOp = ift.DiagonalOperator(f)
diagProber = DiagonalProber(domain=x) diagProber = DiagonalProber(domain=x)
diagProber(diagOp) diagProber(diagOp)
print((f - diagProber.diagonal).norm()) ift.dobj.mprint((f - diagProber.diagonal).norm())
multiProber = MultiProber(domain=x) multiProber = MultiProber(domain=x)
multiProber(diagOp) multiProber(diagOp)
print((f - multiProber.diagonal).norm()) ift.dobj.mprint((f - multiProber.diagonal).norm())
print(f.sum() - multiProber.trace) ift.dobj.mprint(f.sum() - multiProber.trace)
...@@ -3,7 +3,8 @@ import nifty2go as ift ...@@ -3,7 +3,8 @@ import nifty2go as ift
import numericalunits as nu import numericalunits as nu
if __name__ == "__main__": if __name__ == "__main__":
nu.reset_units("SI") # In MPI mode, the random seed for numericalunits must be set by hand
nu.reset_units(43)
dimensionality = 2 dimensionality = 2
np.random.seed(43) np.random.seed(43)
...@@ -44,8 +45,7 @@ if __name__ == "__main__": ...@@ -44,8 +45,7 @@ if __name__ == "__main__":
power_spectrum=power_spectrum) power_spectrum=power_spectrum)
np.random.seed(43) np.random.seed(43)
mock_power = ift.Field(power_space, mock_power = ift.PS_field(power_space, power_spectrum)
val=ift.dobj.from_global_data(power_spectrum(power_space.k_lengths)))
mock_harmonic = ift.power_synthesize(mock_power, real_signal=True) mock_harmonic = ift.power_synthesize(mock_power, real_signal=True)
mock_harmonic = mock_harmonic.real mock_harmonic = mock_harmonic.real
mock_signal = fft(mock_harmonic) mock_signal = fft(mock_harmonic)
...@@ -81,7 +81,7 @@ if __name__ == "__main__": ...@@ -81,7 +81,7 @@ if __name__ == "__main__":
ift.plotting.plot(ift.Field(sspace2, mock_signal.real.val)/nu.K, ift.plotting.plot(ift.Field(sspace2, mock_signal.real.val)/nu.K,
name="mock_signal.pdf") name="mock_signal.pdf")
ift.plotting.plot(ift.Field( data = ift.dobj.to_global_data(data.val.real).reshape(sspace2.shape)/nu.K
sspace2, val=ift.dobj.from_global_data(ift.dobj.to_global_data(data.val.real).reshape(signal_space.shape)))/nu.K, data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))/nu.K
name="data.pdf") ift.plotting.plot(ift.Field(sspace2, val=data), name="data.pdf")
ift.plotting.plot(ift.Field(sspace2, m_s.real.val)/nu.K, name="map.pdf") ift.plotting.plot(ift.Field(sspace2, m_s.real.val)/nu.K, name="map.pdf")
...@@ -50,7 +50,7 @@ if __name__ == "__main__": ...@@ -50,7 +50,7 @@ if __name__ == "__main__":
S = ift.create_power_operator(h_space, power_spectrum=p_spec) S = ift.create_power_operator(h_space, power_spectrum=p_spec)
# Drawing a sample sh from the prior distribution in harmonic space # Drawing a sample sh from the prior distribution in harmonic space
sp = ift.Field(p_space, ift.dobj.from_global_data(p_spec(p_space.k_lengths))) 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.adjoint_times(sh) ss = fft.adjoint_times(sh)