Skip to content
Snippets Groups Projects
Commit d25aa77f authored by Philipp Arras's avatar Philipp Arras
Browse files

8/n

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