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

nifty5 -> nifty7 (6/n)

parent 9b3d6274
No related branches found
No related tags found
1 merge request!2Draft: Nifty5 to nifty7
Pipeline #107685 failed
......@@ -48,14 +48,15 @@ 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, 'gauss', N=N)
plot_prior_samples_2d(5, signal, R, signal, signal.amplitude, '2_gauss', N=N)
likelihood = ift.GaussianEnergy(
mean=data, inverse_covariance=N.inverse)(signal_response)
# Solve inference problem
ic_sampling = ift.GradientNormController(iteration_limit=100)
ic_newton = ift.GradInfNormController(name='Newton', tol=1e-6, iteration_limit=10)
ic_newton = ift.GradInfNormController(name='Newton', tol=1e-6,
iteration_limit=10)
minimizer = ift.NewtonCG(ic_newton)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
......@@ -71,4 +72,5 @@ for _ in range(3):
# Draw posterior samples and plot
N_posterior_samples = 30
KL = ift.MetricGaussianKL(mean, H, N_posterior_samples, True)
plot_reconstruction_2d(data, ground_truth, KL, signal, R, signal.amplitude, 'criticalfilter')
plot_reconstruction_2d(data, ground_truth, KL, signal, R, signal.amplitude,
'2_criticalfilter')
......@@ -15,8 +15,6 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
import nifty7 as ift
from helpers import plot_WF, power_plot, generate_mysterious_data
......@@ -29,44 +27,46 @@ 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
dct = {
'target': power_space,
'n_pix': 64, # 64 spectral bins
# Spectral smoothness (affects Gaussian process part)
'a': 10, # relatively high variance of spectral curvature
'k0': .2, # quefrency mode below which cepstrum flattens
# Power-law part of spectrum:
'sm': -4, # preferred power-law slope
'sv': .6, # low variance of power-law slope
'im': -6, # y-intercept mean, in-/decrease for more/less contrast
'iv': 2. # y-intercept variance
}
A = ift.SLAmplitude(**dct)
args = {
'offset_mean': 0,
'offset_std': (1e-3, 1e-6),
# Amplitude of field fluctuations
'fluctuations': (1., 0.8), # 1.0, 1e-2
correlated_field = ift.CorrelatedField(position_space, A)
# Exponent of power law power spectrum component
'loglogavgslope': (-3., 1), # -6.0, 1
### SETTING UP SPECIFIC SCENARIO ####
# Amplitude of integrated Wiener process power spectrum component
'flexibility': (2, 1.), # 1.0, 0.5
# 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)
# Set up specific scenario
R = ift.GeometryRemover(position_space)
data_space = R.target
signal_response = R(correlated_field)
# Set up likelihood and load data
N = ift.ScalingOperator(data_space, 0.1)
data, ground_truth = generate_mysterious_data(position_space)
data = ift.makeField(data_space, data)
likelihood = ift.GaussianEnergy(mean=data,
inverse_covariance=N.inverse)(signal_response)
likelihood = (ift.GaussianEnergy(mean=data, inverse_covariance=N.inverse)
@ signal_response)
#### SOLVING PROBLEM ####
# Solve problem
ic_sampling = ift.GradientNormController(iteration_limit=100)
ic_newton = ift.GradInfNormController(
name='Newton', tol=1e-6, iteration_limit=30)
ic_newton = ift.GradInfNormController(name='Newton', tol=1e-6,
iteration_limit=10)
minimizer = ift.NewtonCG(ic_newton)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
......@@ -74,34 +74,36 @@ H = ift.StandardHamiltonian(likelihood, ic_sampling)
initial_mean = ift.MultiField.full(H.domain, 0.)
mean = initial_mean
# number of samples used to estimate the KL
# Number of samples used to estimate the KL
N_samples = 10
# Draw new samples to approximate the KL ten times
for i in range(10):
# Draw new samples and minimize KL
KL = ift.MetricGaussianKL(mean, H, N_samples)
KL = ift.MetricGaussianKL(mean, H, N_samples, True)
KL, convergence = minimizer(KL)
mean = KL.position
# FIXME Minimize posterior samples (global)
# Draw posterior samples and plotting
N_posterior_samples = 10
KL = ift.MetricGaussianKL(mean, H, N_posterior_samples)
KL = ift.MetricGaussianKL(mean, H, N_posterior_samples, True)
# Plotting the reconstruction result
# Plot the reconstruction result
ground_truth = ift.makeField(position_space, ground_truth)
posterior_samples = [correlated_field(KL.position+samp) for samp in KL.samples]
mean = 0.*posterior_samples[0]
for p in posterior_samples:
mean = mean + p/len(posterior_samples)
plot_WF('unknown_power', ground_truth, data, m=mean, samples=posterior_samples)
plot_WF('teaser_unknown_power', ground_truth, data, m=mean,
samples=posterior_samples)
# Plotting 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)
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)
power_plot('power_reconstruction', ground_truth_spectrum, power_mean, posterior_power_samples)
power_plot('teaser_power_reconstruction', ground_truth_spectrum, power_mean,
posterior_power_samples)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment