Skip to content

GitLab

Commit 7a785d5e authored by Martin Reinecke's avatar Martin Reinecke
Browse files

make demos look nicer

parent 8e01fbd3
Pipeline #22457 passed with stage
in 4 minutes and 46 seconds
Hide whitespace changes
Inline Side-by-side
demos/critical_filtering.py
View file @ 7a785d5e
......@@ -5,49 +5,14 @@ np.random.seed(42)
# np.seterr(all="raise",under="ignore")
def plot_parameters(m, t, p, p_d):
m = fft.adjoint_times(m)
t = t.val.real
p = p.val.real
p_d = p_d.val.real
ift.plotting.plot(m.real, name='map.png')
class AdjointFFTResponse(ift.LinearOperator):
def __init__(self, FFT, R):
super(AdjointFFTResponse, self).__init__()
self._domain = FFT.target
self._target = R.target
self.R = R
self.FFT = FFT
def _times(self, x):
return self.R(self.FFT.adjoint_times(x))
def _adjoint_times(self, x):
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__":
# Set up position space
s_space = ift.RGSpace([128, 128])
# s_space = ift.HPSpace(32)
# Define harmonic transformation and associated harmonic space
fft = ift.FFTOperator(s_space)
h_space = fft.target[0]
h_space = s_space.get_default_codomain()
fft = ift.FFTOperator(h_space, s_space)
# Set up power space
p_space = ift.PowerSpace(h_space,
......@@ -69,14 +34,12 @@ if __name__ == "__main__":
# Instrument._diagonal.val[64:512-64, 64:512-64] = 0
# Add a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
R = ift.ComposedOperator([fft, Instrument])
noise = 1.
N = ift.DiagonalOperator(ift.Field.full(s_space, noise).weight(1))
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=np.sqrt(noise),
mean=0)
n = ift.Field.from_random(domain=s_space, random_type='normal',
std=np.sqrt(noise), mean=0)
# Create mock data
d = R(sh) + n
......@@ -85,59 +48,47 @@ if __name__ == "__main__":
j = R.adjoint_times(N.inverse_times(d))
realized_power = ift.log(ift.power_analyze(sh,
binbounds=p_space.binbounds))
data_power = ift.log(ift.power_analyze(fft(d),
data_power = ift.log(ift.power_analyze(fft.adjoint_times(d),
binbounds=p_space.binbounds))
d_data = d.val.real
ift.plotting.plot(d.real, name="data.png")
d_data = d.val
ift.plotting.plot(d, name="data.png")
IC1 = ift.GradientNormController(verbose=True, iteration_limit=100,
tol_abs_gradnorm=0.1)
minimizer1 = ift.RelaxedNewton(IC1)
IC2 = ift.GradientNormController(verbose=True, iteration_limit=100,
tol_abs_gradnorm=0.1)
minimizer2 = ift.VL_BFGS(IC2, max_history_length=20)
IC3 = ift.GradientNormController(verbose=True, iteration_limit=1000,
tol_abs_gradnorm=0.1)
minimizer3 = ift.SteepestDescent(IC3)
minimizer = ift.RelaxedNewton(IC1)
ICI = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=0.1)
map_inverter = ift.ConjugateGradient(controller=ICI)
ICI2 = ift.GradientNormController(iteration_limit=200,
tol_abs_gradnorm=1e-5)
power_inverter = ift.ConjugateGradient(controller=ICI2)
# Set starting position
flat_power = ift.Field.full(p_space, 1e-8)
m0 = ift.power_synthesize(flat_power, real_signal=True)
t0 = ift.Field(p_space,
val=ift.dobj.from_global_data(-7.))
t0 = ift.Field(p_space, val=-7.)
for i in range(500):
S0 = ift.create_power_operator(h_space, power_spectrum=ift.exp(t0))
# Initialize non-linear Wiener Filter energy
ICI = ift.GradientNormController(verbose=False, name="ICI",
iteration_limit=500,
tol_abs_gradnorm=0.1)
map_inverter = ift.ConjugateGradient(controller=ICI)
map_energy = ift.library.WienerFilterEnergy(position=m0,
d=d, R=R, N=N, S=S0,
inverter=map_inverter)
map_energy = ift.library.WienerFilterEnergy(
position=m0, d=d, R=R, N=N, S=S0, inverter=map_inverter)
# Solve the Wiener Filter analytically
D0 = map_energy.curvature
m0 = D0.inverse_times(j)
# Initialize power energy with updated parameters
ICI2 = ift.GradientNormController(name="powI",
verbose=False,
iteration_limit=200,
tol_abs_gradnorm=1e-5)
power_inverter = ift.ConjugateGradient(controller=ICI2)
power_energy = ift.library.CriticalPowerEnergy(
position=t0, m=m0, D=D0, smoothness_prior=10., samples=3,
inverter=power_inverter)
(power_energy, convergence) = minimizer1(power_energy)
power_energy = minimizer(power_energy)[0]
# Set new power spectrum
t0 = power_energy.position.real
t0 = power_energy.position
# Plot current estimate
ift.dobj.mprint(i)
if i % 5 == 0:
plot_parameters(m0, t0, ift.log(sp), data_power)
ift.plotting.plot(fft(m0), name='map.png')
demos/nonlinear_critical_filter.py
View file @ 7a785d5e
......@@ -34,12 +34,12 @@ if __name__ == "__main__":
s_spec = ift.Field(p_space, val=1.)
# Choosing the prior correlation structure and defining
# correlation operator
p = ift.Field(p_space, val=p_spec(p_space.k_lengths))
p = ift.PS_field(p_space, p_spec)
log_p = ift.log(p)
S = ift.create_power_operator(h_space, power_spectrum=s_spec)
# Drawing a sample sh from the prior distribution in harmonic space
sp = ift.Field(p_space, val=s_spec)
sp = ift.Field(p_space, val=s_spec)
sh = ift.power_synthesize(sp)
# Choosing the measurement instrument
......@@ -56,19 +56,16 @@ if __name__ == "__main__":
noise_covariance = ift.Field(d_space, val=noise_level**2).weight()
N = ift.DiagonalOperator(noise_covariance)
n = ift.Field.from_random(domain=d_space,
random_type='normal',
std=noise_level,
mean=0)
n = ift.Field.from_random(domain=d_space, random_type='normal',
std=noise_level)
Projection = ift.PowerProjectionOperator(domain=h_space,
power_space=p_space)
power = Projection.adjoint_times(ift.exp(0.5 * log_p))
power = Projection.adjoint_times(ift.exp(0.5*log_p))
# Creating the mock data
true_sky = nonlinearity(FFT.adjoint_times(power * sh))
true_sky = nonlinearity(FFT.adjoint_times(power*sh))
d = MeasurementOperator(true_sky) + n
flat_power = ift.Field(p_space, val=1e-7)
m0 = ift.power_synthesize(flat_power)
m0 = ift.power_synthesize(ift.Field(p_space, val=1e-7))
t0 = ift.Field(p_space, val=-4.)
power0 = Projection.adjoint_times(ift.exp(0.5 * t0))
......@@ -83,7 +80,6 @@ if __name__ == "__main__":
inverter = ift.ConjugateGradient(controller=ICI)
for i in range(20):
power0 = Projection.adjoint_times(ift.exp(0.5*t0))
map0_energy = ift.library.NonlinearWienerFilterEnergy(
m0, d, MeasurementOperator, nonlinearity, FFT, power0, N, S,
......@@ -102,10 +98,9 @@ if __name__ == "__main__":
Instrument=MeasurementOperator, nonlinearity=nonlinearity,
Projection=Projection, sigma=1., samples=2, inverter=inverter)
(power0_energy, convergence) = minimizer(power0_energy)
power0_energy = minimizer(power0_energy)[0]
# Setting new power spectrum
t0_ = t0.copy()
t0 = power0_energy.position
# break degeneracy between amplitude and excitation by setting
......
demos/paper_demos/cartesian_wiener_filter.py
View file @ 7a785d5e
......@@ -5,24 +5,23 @@ if __name__ == "__main__":
signal_to_noise = 1.5 # The signal to noise ratio
# Setting up parameters
correlation_length_1 = 1. # Typical distance over which the field is correlated
field_variance_1 = 2. # Variance of field in position space
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.
L_1 = 2. # Total side-length of the domain
N_pixels_1 = 512 # Grid resolution (pixels per axis)
L_2 = 2. # Total side-length of the domain
N_pixels_2 = 512 # Grid resolution (pixels per axis)
correlation_length_1 = 1.
field_variance_1 = 2. # Variance of field in position space
response_sigma_1 = 0.05 # Smoothing length of response
correlation_length_2 = 1.
field_variance_2 = 2. # Variance of field in position space
response_sigma_2 = 0.01 # Smoothing length of response
def power_spectrum_1(k): # note: field_variance**2 = a*k_0/4.
a = 4 * correlation_length_1 * field_variance_1**2
return a / (1 + k * correlation_length_1) ** 4.
# Setting up the geometry
L_1 = 2. # Total side-length of the domain
N_pixels_1 = 512 # Grid resolution (pixels per axis)
signal_space_1 = ift.RGSpace([N_pixels_1], distances=L_1/N_pixels_1)
harmonic_space_1 = signal_space_1.get_default_codomain()
# Setting up the geometry
L_2 = 2. # Total side-length of the domain
N_pixels_2 = 512 # Grid resolution (pixels per axis)
signal_space_2 = ift.RGSpace([N_pixels_2], distances=L_2/N_pixels_2)
harmonic_space_2 = signal_space_2.get_default_codomain()
......@@ -36,12 +35,6 @@ if __name__ == "__main__":
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)
def power_spectrum_2(k): # note: field_variance**2 = a*k_0/4.
a = 4 * correlation_length_2 * field_variance_2**2
return a / (1 + k * correlation_length_2) ** 2.5
......@@ -53,8 +46,11 @@ if __name__ == "__main__":
fft = ift.ComposedOperator((fft_1, fft_2))
mock_power = ift.Field(domain=(power_space_1, power_space_2),
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))))
mock_power = ift.Field(
(power_space_1, power_space_2),
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 = diagonal.real
......@@ -84,9 +80,9 @@ if __name__ == "__main__":
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag.weight(1))
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
mean=0)
noise = ift.Field.from_random(
domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise
# Wiener filter
......@@ -110,7 +106,13 @@ if __name__ == "__main__":
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
sm = ift.FFTSmoothingOperator(plot_space, sigma=0.03)
ift.plotting.plot(ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))), name='uncertainty.png',zmin=0.,zmax=3.,title="Uncertainty map",colormap="Planck-like")
ift.plotting.plot(ift.Field(plot_space, val=mock_signal.val.real), name='mock_signal.png',colormap="Planck-like")
ift.plotting.plot(ift.Field(plot_space, val=data.val.real), name='data.png',colormap="Planck-like")
ift.plotting.plot(ift.Field(plot_space, val=m.val.real), name='map.png',colormap="Planck-like")
ift.plotting.plot(
ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))),
name='uncertainty.png', zmin=0., zmax=3., title="Uncertainty map",
colormap="Planck-like")
ift.plotting.plot(ift.Field(plot_space, val=mock_signal.val.real),
name='mock_signal.png', colormap="Planck-like")
ift.plotting.plot(ift.Field(plot_space, val=data.val.real),
name='data.png', colormap="Planck-like")
ift.plotting.plot(ift.Field(plot_space, val=m.val.real),
name='map.png', colormap="Planck-like")
demos/paper_demos/wiener_filter.py
View file @ 7a785d5e
......@@ -4,9 +4,12 @@ import numpy as np
if __name__ == "__main__":
# Setting up parameters
correlation_length_scale = 1. # Typical distance over which the field is correlated
L = 2. # Total side-length of the domain
N_pixels = 512 # Grid resolution (pixels per axis)
correlation_length_scale = 1. # Typical distance over which the field is
# correlated
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
signal_to_noise = 1.5 # The signal to noise ratio
np.random.seed(43) # Fixing the random seed
......@@ -14,13 +17,12 @@ if __name__ == "__main__":
a = 4 * correlation_length_scale * fluctuation_scale**2
return a / (1 + (k * correlation_length_scale)**2) ** 2
# Setting up the geometry
L = 2. # Total side-length of the domain
N_pixels = 512 # Grid resolution (pixels per axis)
signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = signal_space.get_default_codomain()
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
S = ift.create_power_operator(harmonic_space,
......@@ -41,8 +43,9 @@ if __name__ == "__main__":
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag.weight(1))
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
noise = ift.Field.from_random(
domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise
# Wiener filter
......@@ -64,7 +67,12 @@ if __name__ == "__main__":
variance = ift.sqrt(sm(proby.diagonal.weight(-1)))
# Plotting
ift.plotting.plot(variance,name="uncertainty.png",xlabel='Pixel index', ylabel='Pixel index')
ift.plotting.plot(mock_signal,name="mock_signal.png",xlabel='Pixel index', ylabel='Pixel index')
ift.plotting.plot(ift.Field(signal_space, val=data.val),name="data.png",xlabel='Pixel index', ylabel='Pixel index')
ift.plotting.plot(m,name="map.png",xlabel='Pixel index', ylabel='Pixel index')
ift.plotting.plot(variance, name="uncertainty.png", xlabel='Pixel index',
ylabel='Pixel index')
ift.plotting.plot(mock_signal, name="mock_signal.png",
xlabel='Pixel index', ylabel='Pixel index')
ift.plotting.plot(ift.Field(signal_space, val=data.val),
name="data.png", xlabel='Pixel index',
ylabel='Pixel index')
ift.plotting.plot(m, name="map.png", xlabel='Pixel index',
ylabel='Pixel index')
nifty/operators/dof_projection_operator.py
View file @ 7a785d5e
......@@ -29,8 +29,8 @@ class DOFProjectionOperator(LinearOperator):
nbin = dofdex.max()
if partner.scalar_dvol() is not None:
wgt = np.bincount(dobj.local_data(dofdex.val).ravel(),
minlength=nbin).astype(np.float64)
wgt *= partner.scalar_dvol()
minlength=nbin)
wgt = wgt*partner.scalar_dvol()
else:
dvol = dobj.local_data(partner.dvol())
wgt = np.bincount(dobj.local_data(dofdex.val).ravel(),
......
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