Commit bf1d9f6e authored by Jakob Knollmueller's avatar Jakob Knollmueller
Browse files

getting_started_1 cconceptually finished, needs polishing

parent ab8fc3d0
import nifty5 as ift
import numpy as np
from global_newton.models_other.apply_data import ApplyData
from global_newton.models_energy.hamiltonian import Hamiltonian
from nifty5.library.unit_log_gauss import UnitLogGauss
def make_chess_mask():
mask = np.ones(position_space.shape)
for i in range(4):
for j in range(4):
if (i+j)%2 == 0:
mask[i*128/4:(i+1)*128/4, j*128/4:(j+1)*128/4] = 0
return mask
def make_random_mask():
mask = ift.from_random('pm1',position_space)
mask = (mask+1)/2
return mask.val
if __name__ == '__main__':
# s_space = ift.RGSpace([1024])
s_space = ift.RGSpace([128,128])
# s_space = ift.HPSpace(64)
## describtion of the tutorial ###
h_space = s_space.get_default_codomain()
total_domain = ift.MultiDomain.make({'xi': h_space})
HT = ift.HarmonicTransformOperator(h_space, s_space)
# Choose problem geometry and masking
def sqrtpspec(k):
return 16. / (20.+k**2)
# # One dimensional regular grid
# position_space = ift.RGSpace([1024])
# mask = np.ones(position_space.shape)
GR = ift.GeometryRemover(s_space)
# # Two dimensional regular grid with chess mask
# position_space = ift.RGSpace([128,128])
# mask = make_chess_mask()
d_space = GR.target
B = ift.FFTSmoothingOperator(s_space,0.1)
mask = np.ones(s_space.shape)
mask[64:89,76:100] = 0.
mask = ift.Field(s_space,val=mask)
Mask = ift.DiagonalOperator(mask)
R = GR * Mask * B
noise = 1.
N = ift.ScalingOperator(noise, d_space)
# # Sphere with half of its locations randomly masked
position_space = ift.HPSpace(128)
mask = make_random_mask()
p_space = ift.PowerSpace(h_space)
pd = ift.PowerDistributor(h_space, p_space)
position = ift.from_random('normal', total_domain)
xi = ift.Variable(position)['xi']
a = ift.Constant(position, ift.PS_field(p_space, sqrtpspec))
A = pd(a)
s_h = A * xi
s = HT(s_h)
Rs = R(s)
# set up corresponding harmonic transform and space
harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
# set correlation structure with a power spectrum and build prior correlation covariance
def power_spectrum(k):
return 100. / (20.+k**3)
power_space = ift.PowerSpace(harmonic_space)
PD = ift.PowerDistributor(harmonic_space, power_space)
prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
S = ift.DiagonalOperator(prior_correlation_structure)
MOCK_POSITION = ift.from_random('normal',total_domain)
data = Rs.at(MOCK_POSITION).value + N.draw_sample()
# build instrument response consisting of a discretization, mask and harmonic transformaion
GR = ift.GeometryRemover(position_space)
mask = ift.Field(position_space,val=mask)
Mask = ift.DiagonalOperator(mask)
R = GR * Mask * HT
data_space = GR.target
# setting the noise covariance
noise = 5.
N = ift.ScalingOperator(noise, data_space)
NWR = ApplyData(data, ift.Field(d_space,val=noise), Rs)
INITIAL_POSITION = ift.from_random('normal',total_domain)
likelihood = UnitLogGauss(INITIAL_POSITION, NWR)
# creating mock data
MOCK_SIGNAL = S.draw_sample()
MOCK_NOISE = N.draw_sample()
data = R(MOCK_SIGNAL) + MOCK_NOISE
# building propagator D and information source j
j = R.adjoint_times(N.inverse_times(data))
D_inv = R.adjoint * N.inverse * R + S.inverse
# make it invertible
IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
inverter = ift.ConjugateGradient(controller=IC)
IC2 = ift.GradientNormController(name='Newton', iteration_limit=15)
minimizer = ift.RelaxedNewton(IC2)
D = ift.InversionEnabler(D_inv,IC,approximation=S.inverse).inverse
# WIENER FILTER
m = D(j)
##PLOTTING
#Truth, data, reconstruction, residuals
H = Hamiltonian(likelihood, inverter)
H, convergence = minimizer(H)
result = s.at(H.position).value
import nifty5 as ift
import sys
import numpy as np
import global_newton as gn
from nifty5.library.nonlinearities import Exponential
# from nifty5.library.nonlinearities import Exponential
#DEFINE THE SIGNAL:
#Define signal space as a regular grid
#s_space = ift.RGSpace(10024)
s_space = ift.RGSpace([128,128])
#s_space = ift.RGSpace(1024)
# s_space = ift.RGSpace([128,128])
s_space = ift.HPSpace(128)
#Define the harmonic space
h_space = s_space.get_default_codomain()
......@@ -21,10 +19,10 @@ domain = ift.MultiDomain.make({'xi': h_space})
#Define positions from a Gaussian distribution
position = ift.from_random('normal', domain)
Nsamples = 5
#Define a power spectrum
def sqrtpspec(k):
return 16. / (20.+k**2)
return 10. / (20.+k**2)
#Define a power space
p_space = ift.PowerSpace(h_space)
......@@ -53,17 +51,17 @@ sky = ift.PointwiseExponential(logsky)
#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.
exposure = 1*np.ones(s_space.shape)
# exposure[int(s_space.shape[0]/3):int(s_space.shape[0]/3+10)] = 10.
#Convert the mask into a field
mask = ift.Field(s_space,val=mask)
exposure = ift.Field(s_space,val=exposure)
#Create a diagonal matrix corresponding to the mask
M = ift.DiagonalOperator(mask)
E = ift.DiagonalOperator(exposure)
#Create the response operator and apply the mask on it
R = ift.ScalingOperator(1., s_space) * M
R = ift.GeometryRemover(s_space) * E
#CREATE THE MOCK DATA:
......@@ -94,41 +92,20 @@ 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)
#Update the position
position = KL.position
#Obtain the value of the KL and the convergence at the new position
KL, convergence = minimizer(KL)
H = ift.Hamiltonian(likelihood, ic_cg)
H, convergence = minimizer(H)
#PLOT RESULTS:
#Evaluate lambda at final position
lamb_recontructed = lamb.at(KL.position).value
lamb_recontructed = lamb.at(H.position).value
sky_reconstructed = sky.at(H.position).value
#Evaluate lambda at the data position for comparison
lamb_mock = lamb.at(mock_position).value
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment