Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • 420-tracerboolconversion-error-in-lognormal_moments-py
  • 423-minisanity-re-improve-likelihood-readability
  • NIFTy_1
  • NIFTy_2
  • NIFTy_3
  • NIFTy_4
  • NIFTy_5
  • NIFTy_6
  • NIFTy_7
  • NIFTy_8
  • TransitionFunction
  • adaptions_cfm_veberle
  • altVI
  • cf_matern_N_copies
  • cfm_demo_sliders
  • change_ncg_default
  • code_formatter
  • config_multifrequency_trans
  • cupy_backend
  • cupy_backend_stable
  • feature_nifty_re_static_newton_cg
  • fft_interpolator
  • fix_cfm_amplitude_normalization
  • fix_cfm_amplitude_normalization_nifty8
  • fix_color_mf_plotting
  • fix_invgamma_prior_re
  • fix_nonlinearity_gradients
  • fix_saved_pickle_files
  • flexible_FieldZeroPadder
  • force_precompute_dists
  • format_to_fstring
  • frequency_model
  • gauss_markov_processes
  • general_sky
  • geomap_pf
  • hammer_projection
  • hmc
  • improve_getting_started_0
  • independent_SLAmplitudes
  • inline_issue
  • introduce_qol_operators
  • jaxopt_vi
  • lanczos
  • likelihood_constant_terms
  • likelihoods
  • main
  • mf_plotting_refactor
  • mf_sky
  • more_derivatives_pf
  • more_restrictive_scaling_operator
  • mpcdf_runners
  • msc
  • new_sampling
  • nifty
  • nifty2jax_david_vmap
  • nifty8_philipps_unmerged_patches
  • nifty_jr
  • options_for_zeropadder
  • paper_lsp
  • paper_matteo
  • paper_pf
  • paper_ve
  • parametric_kl
  • playground_iwp_mf
  • projection_operators_maxim
  • qpo_model
  • restructuring_proposal
  • rmmap
  • symbolicDifferentiation
  • test_add_fluctuations
  • tests_jhk
  • tmp_mf_plotting
  • total_variation_prior
  • yango_minimizer
  • yapf_01
  • 0.2.0
  • 0.3.0
  • 0.4.0
  • 0.4.2
  • 0.6.0
  • 0.7.0
  • 0.8.0
  • 0.8.4
  • 0.9.0
  • 1.0.6
  • 1.0.7
  • 1.0.8
  • 4.0.0
  • 4.1.0
  • 4.2.0
  • 5.0.0
  • 6.0_kern
  • historical/head_RXTE_compatible
  • historical/head_field_testing
  • historical/head_master
  • historical/head_plotting
  • historical/head_theo_master
  • historical/nifty2go
  • historical/theo_master
  • point_source_separator
  • v5.0.1
  • v7.0
  • v7.1
  • v7.2
  • v7.3
  • v7.4
  • v7.5
  • v8.0
  • v8.1
  • v8.2
  • v8.3
  • v8.4
  • v8.5
  • v8.5.1
  • v8.5.2
  • v8.5.3
  • v8.5.4
  • v8.5.5
  • v8.5.6
  • v8.5.7
120 results

Target

Select target project
  • ift/nifty
  • g-philipparras/NIFTy
  • tpeters/NIFTy
  • g-neelshah/nifty
  • philiar/nifty
5 results
Select Git revision
  • 420-tracerboolconversion-error-in-lognormal_moments-py
  • 423-minisanity-re-improve-likelihood-readability
  • NIFTy_1
  • NIFTy_2
  • NIFTy_3
  • NIFTy_4
  • NIFTy_5
  • NIFTy_6
  • NIFTy_7
  • NIFTy_8
  • TransitionFunction
  • adaptions_cfm_veberle
  • altVI
  • cf_matern_N_copies
  • cfm_demo_sliders
  • change_ncg_default
  • code_formatter
  • config_multifrequency_trans
  • cupy_backend
  • cupy_backend_stable
  • feature_nifty_re_static_newton_cg
  • fft_interpolator
  • fix_cfm_amplitude_normalization
  • fix_cfm_amplitude_normalization_nifty8
  • fix_color_mf_plotting
  • fix_invgamma_prior_re
  • fix_nonlinearity_gradients
  • fix_saved_pickle_files
  • flexible_FieldZeroPadder
  • force_precompute_dists
  • format_to_fstring
  • frequency_model
  • gauss_markov_processes
  • general_sky
  • geomap_pf
  • hammer_projection
  • hmc
  • improve_getting_started_0
  • independent_SLAmplitudes
  • inline_issue
  • introduce_qol_operators
  • jaxopt_vi
  • lanczos
  • likelihood_constant_terms
  • likelihoods
  • main
  • mf_plotting_refactor
  • mf_sky
  • more_derivatives_pf
  • more_restrictive_scaling_operator
  • mpcdf_runners
  • msc
  • new_sampling
  • nifty
  • nifty2jax_david_vmap
  • nifty8_philipps_unmerged_patches
  • nifty_jr
  • options_for_zeropadder
  • paper_lsp
  • paper_matteo
  • paper_pf
  • paper_ve
  • parametric_kl
  • playground_iwp_mf
  • projection_operators_maxim
  • qpo_model
  • restructuring_proposal
  • rmmap
  • symbolicDifferentiation
  • test_add_fluctuations
  • tests_jhk
  • tmp_mf_plotting
  • total_variation_prior
  • yango_minimizer
  • yapf_01
  • 0.2.0
  • 0.3.0
  • 0.4.0
  • 0.4.2
  • 0.6.0
  • 0.7.0
  • 0.8.0
  • 0.8.4
  • 0.9.0
  • 1.0.6
  • 1.0.7
  • 1.0.8
  • 4.0.0
  • 4.1.0
  • 4.2.0
  • 5.0.0
  • 6.0_kern
  • historical/head_RXTE_compatible
  • historical/head_field_testing
  • historical/head_master
  • historical/head_plotting
  • historical/head_theo_master
  • historical/nifty2go
  • historical/theo_master
  • point_source_separator
  • v5.0.1
  • v7.0
  • v7.1
  • v7.2
  • v7.3
  • v7.4
  • v7.5
  • v8.0
  • v8.1
  • v8.2
  • v8.3
  • v8.4
  • v8.5
  • v8.5.1
  • v8.5.2
  • v8.5.3
  • v8.5.4
  • v8.5.5
  • v8.5.6
  • v8.5.7
120 results
Show changes
Commits on Source (1)
......@@ -37,10 +37,11 @@ D_inv = ift.SandwichOperator.make(R_p, N.inverse) + S.inverse
N_samps = 200
N_iter = 100
IC = ift.GradientNormController(tol_abs_gradnorm=1e-3, iteration_limit=N_iter)
m, samps = ift.library.generate_krylov_samples(D_inv, S, j, N_samps, IC)
m_x = sky(m)
samps = ift.library.generate_krylov_samples(D_inv, S, N_samps, IC)
inverter = ift.ConjugateGradient(IC)
curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p, inverter=inverter)
m = curv.inverse_times(j)
m_x = sky(m)
samps_old = [curv.draw_sample(from_inverse=True) for i in range(N_samps)]
plt.plot(d.to_global_data(), '+', label="data", alpha=.5)
......
......@@ -20,26 +20,23 @@ import numpy as np
from ..minimization.quadratic_energy import QuadraticEnergy
def generate_krylov_samples(D_inv, S, j, N_samps, controller):
def generate_krylov_samples(D_inv, S, N_samps, controller):
"""
Generates inverse samples from a curvature D.
This algorithm iteratively generates samples from
a curvature D by applying conjugate gradient steps
and resampling the curvature in search direction.
It is basically just a more stable version of
Wiener Filter samples
Parameters
----------
D_inv : EndomorphicOperator
The curvature which will be the inverse of the covarianc
D_inv : WienerFilterCurvature
The curvature which will be the inverse of the covariance
of the generated samples
S : EndomorphicOperator (from which one can sample)
A prior covariance operator which is used to generate prior
samples that are then iteratively updated
j : Field, optional
A Field to which the inverse of D_inv is applied. The solution
of this matrix inversion problem is a side product of generating
the samples.
If not supplied, it is sampled from the inverse prior.
N_samps : Int
How many samples to generate.
controller : IterationController
......@@ -47,56 +44,62 @@ def generate_krylov_samples(D_inv, S, j, N_samps, controller):
Returns
-------
(solution, samples) : A tuple of a field 'solution' and a list of fields
'samples'. The first entry of the tuple is the solution x to
D_inv(x) = j
and the second entry are a list of samples from D_inv.inverse
samples : a list of samples from D_inv.inverse
"""
# RL FIXME: make consistent with complex numbers
j = S.draw_sample(from_inverse=True) if j is None else j
energy = QuadraticEnergy(j.empty_copy().fill(0.), D_inv, j)
y = [S.draw_sample() for _ in range(N_samps)]
status = controller.start(energy)
if status != controller.CONTINUE:
return energy.position, y
r = energy.gradient
d = r.copy()
previous_gamma = r.vdot(r).real
if previous_gamma == 0:
return energy.position, y
while True:
q = energy.curvature(d)
ddotq = d.vdot(q).real
if ddotq == 0.:
logger.error("Error: ConjugateGradient: ddotq==0.")
return energy.position, y
alpha = previous_gamma/ddotq
if alpha < 0:
logger.error("Error: ConjugateGradient: alpha<0.")
return energy.position, y
for i in range(len(y)):
y[i] += (np.random.randn()*np.sqrt(ddotq) - y[i].vdot(q))/ddotq * d
q *= -alpha
r = r + q
energy = energy.at_with_grad(energy.position - alpha*d, r)
gamma = r.vdot(r).real
if gamma == 0:
return energy.position, y
status = controller.check(energy)
samples = []
for i in range(N_samps):
x0 = S.draw_sample()
y = x0*0
j = y*0
#j = y
energy = QuadraticEnergy(x0, D_inv, j)
status = controller.start(energy)
if status != controller.CONTINUE:
return energy.position, y
d *= max(0, gamma/previous_gamma)
d += r
previous_gamma = gamma
samples += [y]
break
r = energy.gradient
d = r.copy()
previous_gamma = r.vdot(r).real
if previous_gamma == 0:
samples += [y+energy.position]
break
while True:
q = energy.curvature(d)
ddotq = d.vdot(q).real
if ddotq == 0.:
logger.error("Error: ConjugateGradient: ddotq==0.")
samples += [y+energy.position]
break
alpha = previous_gamma/ddotq
if alpha < 0:
logger.error("Error: ConjugateGradient: alpha<0.")
samples += [y+energy.position]
break
y += (np.random.randn()*np.sqrt(ddotq) )/ddotq * d
q *= -alpha
r = r + q
energy = energy.at_with_grad(energy.position - alpha*d, r)
gamma = r.vdot(r).real
if gamma == 0:
samples += [y+energy.position]
break
status = controller.check(energy)
if status != controller.CONTINUE:
samples += [y+energy.position]
break
d *= max(0, gamma/previous_gamma)
d += r
previous_gamma = gamma
return samples