diff --git a/2_critical_filter_solution.py b/2_critical_filter_solution.py index e2f422538ef020bcf0c8af3e76c102a033c059a4..0fe20015649a8f58621ed34d0e852b75106eae42 100644 --- a/2_critical_filter_solution.py +++ b/2_critical_filter_solution.py @@ -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') diff --git a/teaser_critical_filter.py b/teaser_critical_filter.py index 55c4441686ce74e90433d99e2222c43018266bf8..a0cee0ee7f0add374974aa981f484addf170ad99 100644 --- a/teaser_critical_filter.py +++ b/teaser_critical_filter.py @@ -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)