diff --git a/2_critical_filter_solution.py b/2_critical_filter_solution.py
index 8e78e3eae010a1aeea331995e27a53e94d877a22..351e084151883ff6bdd5366cfb1821043a504680 100644
--- a/2_critical_filter_solution.py
+++ b/2_critical_filter_solution.py
@@ -48,7 +48,8 @@ signal_response = R @ signal
N = ift.ScalingOperator(data_space, 0.1)
data, ground_truth = generate_gaussian_data(signal_response, N)
-plot_prior_samples_2d(5, signal, R, signal, signal.amplitude, '2_gauss', N=N)
+plot_prior_samples_2d(5, signal, R, signal, signal.power_spectrum, 'gauss',
+ N=N)
likelihood = ift.GaussianEnergy(
mean=data, inverse_covariance=N.inverse)(signal_response)
@@ -72,5 +73,5 @@ for ii in range(n_rounds):
mean = KL.position
# Plot posterior samples
-plot_reconstruction_2d(data, ground_truth, KL, signal, R, signal.amplitude,
- '2_criticalfilter')
+plot_reconstruction_2d(data, ground_truth, KL, signal, R,
+ signal.power_spectrum, '2_criticalfilter')
diff --git a/3_more_examples.py b/3_more_examples.py
index a82507b71a9175dce8aed7a807548c5b0c52a56f..ca9c0646c3b45d3da5e54149ffde69c96611de52 100644
--- a/3_more_examples.py
+++ b/3_more_examples.py
@@ -45,7 +45,7 @@ for mode in [0, 1]:
'asperity': (0.5, 0.4) # 0.1, 0.5
}
correlated_field = ift.SimpleCorrelatedField(position_space, **args)
- A = correlated_field.amplitude
+ pspec = correlated_field.power_spectrum
dct = {}
if mode == 0:
@@ -54,7 +54,7 @@ for mode in [0, 1]:
else:
signal = correlated_field.exp()
R = h.exposure_response(position_space)
- h.plot_prior_samples_2d(5, signal, R, correlated_field, A, name[mode],
+ h.plot_prior_samples_2d(5, signal, R, correlated_field, pspec, name[mode],
**dct)
signal_response = R @ signal
if mode == 0:
@@ -79,4 +79,5 @@ for mode in [0, 1]:
KL = ift.MetricGaussianKL(mean, H, 30 if ii == n_rounds-1 else 5, True)
KL, convergence = minimizer(KL)
mean = KL.position
- h.plot_reconstruction_2d(data, ground_truth, KL, signal, R, A, name[mode])
+ h.plot_reconstruction_2d(data, ground_truth, KL, signal, R, pspec,
+ name[mode])
diff --git a/helpers/plot.py b/helpers/plot.py
index 2f2c3da676f5eaec5602943cc81b9561acac2697..03c452c7f68e10b6d828037d1b1288b28b43a233 100644
--- a/helpers/plot.py
+++ b/helpers/plot.py
@@ -85,10 +85,9 @@ def power_plot(name, s, m, samples=None):
plt.close('all')
-def plot_prior_samples_2d(n_samps, signal, R, correlated_field, A, likelihood,
- N=None):
+def plot_prior_samples_2d(n_samps, signal, R, correlated_field, pspec,
+ likelihood, N=None):
samples, pspecmin, pspecmax = [], np.inf, 0
- pspec = A*A
for _ in range(n_samps):
ss = ift.from_random(signal.domain)
samples.append(ss)
@@ -118,8 +117,7 @@ def plot_prior_samples_2d(n_samps, signal, R, correlated_field, A, likelihood,
raise ValueError('likelihood type not implemented')
data = R.adjoint(data + 0.)
- As = pspec.force(sample)
- ax[ii, 0].plot(As.domain[0].k_lengths, As.val)
+ ax[ii, 0].plot(pspec.target[0].k_lengths, pspec.force(sample))
ax[ii, 0].set_yscale('log')
ax[ii, 0].set_xscale('log')
ax[ii, 0].set_ylim(pspecmin, pspecmax)
@@ -155,14 +153,14 @@ def plot_prior_samples_2d(n_samps, signal, R, correlated_field, A, likelihood,
plt.close('all')
-def plot_reconstruction_2d(data, ground_truth, KL, signal, R, A, name):
+def plot_reconstruction_2d(data, ground_truth, KL, signal, R, pspec, name):
sc = ift.StatCalculator()
sky_samples, pspec_samples = [], []
for sample in KL.samples:
tmp = signal(sample + KL.position)
sc.add(tmp)
sky_samples.append(tmp)
- pspec_samples.append(A.force(sample)**2)
+ pspec_samples.append(pspec.force(sample))
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(4*3, 4*2))
im = []
@@ -198,7 +196,7 @@ def plot_reconstruction_2d(data, ground_truth, KL, signal, R, A, name):
label='reconstruction')
ax[1, 2].plot(
amp_mean.domain[0].k_lengths,
- A.force(ground_truth).val**2,
+ pspec.force(ground_truth).val,
color='b',
label='ground truth')
ax[1, 2].legend()
diff --git a/teaser_critical_filter.py b/teaser_critical_filter.py
index 468caebf8eb9808adf365679b81ab0af711ac6b9..1c13b7772aa6252cb25ed0b310e1c87bd4fd4e2f 100644
--- a/teaser_critical_filter.py
+++ b/teaser_critical_filter.py
@@ -15,8 +15,12 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+from functools import reduce
+from operator import add
+
import nifty7 as ift
-from helpers import plot_WF, power_plot, generate_mysterious_data
+
+from helpers import generate_mysterious_data, plot_WF, power_plot
position_space = ift.RGSpace(256)
harmonic_space = position_space.get_default_codomain()
@@ -25,8 +29,8 @@ HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
power_space = ift.PowerSpace(harmonic_space)
-# Set up an amplitude operator for the field
-# We want to set up a model for the amplitude spectrum with some magic numbers
+# Set up prior model for the field for the field
+# We want to set up a model for the power spectrum with some magic numbers
args = {
'offset_mean': 0,
'offset_std': (1e-3, 1e-6),
@@ -43,10 +47,8 @@ args = {
# How ragged the integrated Wiener process component is
'asperity': (0.5, 0.4) # 0.1, 0.5
}
-
correlated_field = ift.SimpleCorrelatedField(position_space, **args)
-A = correlated_field.amplitude
-# FIXME amplitude -> power spectrum (global)
+pspec = correlated_field.power_spectrum
# Set up specific scenario
R = ift.GeometryRemover(position_space)
@@ -95,10 +97,10 @@ plot_WF('teaser_unknown_power', ground_truth, data, m=mean,
# Plot the reconstruction of the power spectrum
mysterious_spectrum = lambda k: 5/((7**2 - k**2)**2 + 3**2*k**2)
-ground_truth_spectrum = ift.makeField(power_space, mysterious_spectrum(power_space.k_lengths))
-posterior_power_samples = [A.force(KL.position+samp)**2 for samp in KL.samples]
-power_mean = 0.*posterior_power_samples[0]
-for p in posterior_power_samples:
- power_mean = power_mean + p/len(posterior_power_samples)
+# FIXME There is a simpler way, isn't it?
+ground_truth_spectrum = ift.makeField(power_space,
+ mysterious_spectrum(power_space.k_lengths))
+posterior_power_samples = [pspec.force(KL.position+samp) for samp in KL.samples]
+power_mean = reduce(add, posterior_power_samples)
power_plot('teaser_power_reconstruction', ground_truth_spectrum, power_mean,
posterior_power_samples)