From ab8fc3d0b78fd10adcb2060c472e26f5cb558f47 Mon Sep 17 00:00:00 2001 From: natalia Date: Thu, 28 Jun 2018 12:01:17 +0200 Subject: [PATCH] Add comments --- demos/getting_started_2.py | 116 ++++++++++++++++++++++++++++--------- 1 file changed, 89 insertions(+), 27 deletions(-) diff --git a/demos/getting_started_2.py b/demos/getting_started_2.py index 124df9f0..1ddbbabc 100644 --- a/demos/getting_started_2.py +++ b/demos/getting_started_2.py @@ -4,80 +4,142 @@ import numpy as np import global_newton as gn from nifty5.library.nonlinearities import Exponential -np.random.seed(42) +#DEFINE THE SIGNAL: -N = 2 -Nsamples = 5 +#Define signal space as a regular grid +#s_space = ift.RGSpace(10024) +s_space = ift.RGSpace([128,128]) -s_space = ift.RGSpace(10024) +#Define the harmonic space h_space = s_space.get_default_codomain() +#Prepare Harmonic transformation between the two spaces +HT = ift.HarmonicTransformOperator(h_space, s_space) + +#Define domain domain = ift.MultiDomain.make({'xi': h_space}) + +#Define positions from a Gaussian distribution position = ift.from_random('normal', domain) -HT = ift.HarmonicTransformOperator(h_space, s_space) +#Define a power spectrum def sqrtpspec(k): return 16. / (20.+k**2) - -# Define amplitude model +#Define a power space p_space = ift.PowerSpace(h_space) + +#Define the power distribution between the harmonic and the power spaces pd = ift.PowerDistributor(h_space, p_space) + +#Create a field with the defined power spectrum a = ift.PS_field(p_space, sqrtpspec) + +#Define the amplitudes A = pd(a) -# Define sky model +#Unpack the positions xi from the Multifield xi = ift.Variable(position)['xi'] + +#Multiply the positions by the amplitudes in the harmonic domain logsky_h = xi * A + +#Transform to the real domain logsky = HT(logsky_h) -nonlin = Exponential() + +#Create a sky model by applying the exponential (Poisson) sky = ift.PointwiseExponential(logsky) -R = ift.ScalingOperator(1., s_space) + +#DEFINE THE RESPONSE OPERATOR: + +#Define a mask to cover a patch of the real space +mask = np.ones(s_space.shape) +mask[int(s_space.shape[0]/3):int(s_space.shape[0]/3+10)] = 0. + +#Convert the mask into a field +mask = ift.Field(s_space,val=mask) + +#Create a diagonal matrix corresponding to the mask +M = ift.DiagonalOperator(mask) + +#Create the response operator and apply the mask on it +R = ift.ScalingOperator(1., s_space) * M + +#CREATE THE MOCK DATA: + +#Define the data space d_space = R.target[0] + +#Apply the response operator to the signal +#lamb corresponds to the mean in the Poisson distribution lamb = R(sky) -# Generate mock data -MOCK_POSITION = ift.from_random('normal', lamb.position.domain) -data = np.random.poisson(lamb.at(MOCK_POSITION).value.val.astype(np.float64)) +#Draw coordinates of the mock data from a Gaussian distribution +mock_position = ift.from_random('normal', lamb.position.domain) + +#Generate mock data from a Poisson distribution using lamb as a mean +data = np.random.poisson(lamb.at(mock_position).value.val.astype(np.float64)) + +#Store the data as a field data = ift.Field.from_local_data(d_space, data) -# Define Hamiltonian +#RECONSTRUCT THE SIGNAL: + +#Define the positions where we perform the analysis from a Gaussian distribution position = ift.from_random('normal', lamb.position.domain) + +#Define the Poisson likelihood knowing the mean and the data likelihood = ift.library.PoissonLogLikelihood(lamb, data) +#Define a iteration controller with a maximum number of iterations ic_cg = ift.GradientNormController(iteration_limit=50) + +#Define a iteration controller with convergence criteria ic_samps = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-4) +#Define a iteration controller for the minimizer ic_newton = ift.GradientNormController(name='Newton', tol_abs_gradnorm=1e-3) minimizer = ift.RelaxedNewton(ic_newton) +#Build the Hamiltonian H = gn.Hamiltonian(likelihood, ic_cg, ic_samps) +#Iterate the reconstruction for _ in range(N): + #Draw samples from the curvature samples = [H.curvature.draw_sample(from_inverse=True) for _ in range(Nsamples)] + sc_samplesky = ift.StatCalculator() for s in samples: sc_samplesky.add(sky.at(s+position).value) ift.plot(sc_samplesky.mean, name='sample_mean.png') + #Compute the Kullbachh Leibler divergence KL = gn.SampledKullbachLeiblerDivergence(H, samples, ic_cg) - KL, convergence = minimizer(KL) + + #Update the position position = KL.position + + #Obtain the value of the KL and the convergence at the new position + KL, convergence = minimizer(KL) + +#PLOT RESULTS: + +#Evaluate lambda at final position +lamb_recontructed = lamb.at(KL.position).value + +#Evaluate lambda at the data position for comparison +lamb_mock = lamb.at(mock_position).value -# Plot results -E = KL -l1 = lamb.at(E.position).value -l2 = lamb.at(MOCK_POSITION).value -ift.plot([data, l2, l1], name="poisson.png", +#Plot the data, reconstruction and underlying signal +ift.plot([data, lamb_mock, lamb_recontructed], name="poisson.png", label=['Data', 'Mock signal', 'Reconstruction'], alpha=[.5, 1, 1]) -if power_spectrum_estimation: - a_mock = a.at(MOCK_POSITION).value - a_recon = a.at(E.position).value -else: - a_mock = a - a_recon = a -ift.plot([a_mock**2, a_recon**2, ift.power_analyze(logsky_h.at(E.position).value)], + +#Plot power spectrum for posterior test +a_mock = a +a_recon = a +ift.plot([a_mock**2, a_recon**2, ift.power_analyze(logsky_h.at(KL.position).value)], name='power_spectrum.png', label=['Mock', 'Reconstruction', 'power_analyze']) -- GitLab