Commit a4217c4e authored by Lukas Platz's avatar Lukas Platz

demo with new ansatz (still needs improvement)

parent 583ae212
Pipeline #61562 failed with stages
in 4 minutes and 11 seconds
......@@ -55,8 +55,10 @@ if __name__ == '__main__':
if not mode in (0, 1):
raise ValueError("Only modes 0 and 1 defined")
filename = f"test-MultiCorrelatedField-{mode}" + "-{}.png"
# Two-dimensional signal on RGs with inhomogeneous exposure
N = 512
N = 128
position_space = ift.RGSpace(N)
energy_space = ift.RGSpace(N)
sky_target = ift.DomainTuple.make((position_space, energy_space))
......@@ -78,7 +80,7 @@ if __name__ == '__main__':
**{
'target': power_space_p,
'n_pix': 64, # 64 spectral bins
'keys': ['tau_p', 'phi_p'],
'keys': ['tau_p', 'phi_p', 'zm_p'],
# Spectral smoothness (affects Gaussian process part)
'a': 3, # relatively high variance of spectral curvature
......@@ -88,14 +90,18 @@ if __name__ == '__main__':
'sm': -4., # Expected exponent of power law
'sv': 1., # Prior standard deviation of -||-
'im': 0., # Expected y-intercept of power law
'iv': 0. # Prior standard deviation of -||-
'iv': 0., # Prior standard deviation of -||-
# zero-mode
'za': 5., # IG alpha parameter
'zq': 0.001 # IG q parameter
})
amp_e = ift.SLAmplitude(
**{
'target': power_space_e,
'n_pix': 64, # 64 spectral bins
'keys': ['tau_e', 'phi_e'],
'keys': ['tau_e', 'phi_e', 'zm_e'],
# Spectral smoothness (affects Gaussian process part)
'a': 3, # relatively high variance of spectral curvature
......@@ -105,7 +111,11 @@ if __name__ == '__main__':
'sm': -3., # Expected exponent of power law
'sv': 1., # Prior standard deviation of -||-
'im': 0., # Expected y-intercept of power law
'iv': 0. # Prior standard deviation of -||-
'iv': 0., # Prior standard deviation of -||-
# zero-mode
'za': 6., # IG alpha parameter
'zq': 0.002 # IG q parameter
})
# Define sky operator
......@@ -117,8 +127,8 @@ if __name__ == '__main__':
log_diffuse = ift.MultiCorrelatedField(sky_target,
amplitudes,
va=5.,
vq=1.)
gs_a=11.,
gs_q=50.)
sky = log_diffuse.exp()
M = ift.DiagonalOperator(exposure)
......@@ -130,6 +140,8 @@ if __name__ == '__main__':
lamb = R(sky)
mock_position = ift.from_random('normal', sky.domain)
data = lamb(mock_position)
signal = sky(mock_position)
# DEBUG
mock_log_diffuse_vals = log_diffuse(mock_position).to_global_data()
......@@ -140,47 +152,58 @@ if __name__ == '__main__':
likelihood = ift.PoissonianEnergy(data)(lamb)
# Settings for minimization
ic_newton = ift.DeltaEnergyController(name='Newton',
iteration_limit=100,
tol_rel_deltaE=1e-8)
ic_sampling = ift.AbsDeltaEnergyController(name='Sampling',
iteration_limit=200,
deltaE=0.5)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
ic_newton = ift.AbsDeltaEnergyController(name='Newton',
iteration_limit=30,
deltaE=0.5)
minimizer = ift.NewtonCG(ic_newton)
pos = ift.from_random('normal', sky.domain) * 0.01
iteration = 0
# Minimize
for i in range(3):
KL = ift.MetricGaussianKL(pos, H, 10)
KL, convergence = minimizer(KL)
pos = KL.position
iteration += 1
# Plotting
sc = ift.StatCalculator()
for sample in KL.samples:
sc.add(sky(pos + sample))
reconst = sc.mean
# Compute MAP solution by minimizing the information Hamiltonian
H = ift.StandardHamiltonian(likelihood)
initial_position = ift.from_random('normal', sky.domain) * 0.01
H = ift.EnergyAdapter(initial_position, H, want_metric=True)
H, convergence = minimizer(H)
# Plotting
filename = f"test-MultiCorrelatedField-{mode}" + "-{}.png"
signal = sky(mock_position)
reconst = sky(H.position)
plot = ift.Plot()
plot.add(cd(signal), title='Signal')
plot.add(cd(GR.adjoint(data)), title='Data')
plot.add(cd(reconst), title='Reconstruction')
plot.add(cd(reconst - signal), title='Error')
plot.output(xsize=12, ysize=10, name=filename.format('signal'))
p_p = ift.PowerDistributor(harmonic_space_p, power_space_p) @ amp_p
p_e = ift.PowerDistributor(harmonic_space_e, power_space_e) @ amp_e
if mode == 1:
amplitude_specs = [(p_p, 'position'), (p_e, 'energy')]
elif mode == 0:
amplitude_specs = [(p_e, 'energy')]
for p, n in amplitude_specs:
plot = ift.Plot()
plot.add([
ift.power_analyze(p.force(H.position)),
ift.power_analyze(p.force(mock_position))
],
linewidth=[1., 3.],
label=['Reconstruction', 'Ground Truth'],
title=f"Posterior {n} spectrum")
plot.output(xsize=12, ysize=10, name=filename.format('power_' + n))
print("Saved results as '{}'.".format(filename.format('*')))
plot.add(cd(signal), title='Signal')
plot.add(cd(GR.adjoint(data)), title='Data')
plot.add(cd(reconst), title='Posterior Mean')
plot.add(cd(reconst - signal), title='Error')
plot.output(xsize=12, ysize=10, name=filename.format(f'reconstr-{iteration:02d}'))
p_p = ift.PowerDistributor(harmonic_space_p, power_space_p) @ amp_p
p_e = ift.PowerDistributor(harmonic_space_e, power_space_e) @ amp_e
if mode == 1:
amplitude_specs = [(p_p, 'position'), (p_e, 'energy')]
elif mode == 0:
amplitude_specs = [(p_e, 'energy')]
for p, n in amplitude_specs:
sc = ift.StatCalculator()
for sample in KL.samples:
sc.add(ift.power_analyze(p.force(pos + sample)))
plot = ift.Plot()
plot.add([
sc.mean,
ift.power_analyze(p.force(mock_position))
],
linewidth=[1., 3.],
label=['Posterior Mean', 'Ground Truth'],
title=f"Posterior {n} spectrum")
plot.output(xsize=12, ysize=10, name=filename.format(f'power_{n}-{iteration:02d}'))
print("Saved results as '{}'.".format(filename.format('*')))
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