...
 
Commits (4)
...@@ -55,7 +55,7 @@ if __name__ == '__main__': ...@@ -55,7 +55,7 @@ if __name__ == '__main__':
if not mode in (0, 1): if not mode in (0, 1):
raise ValueError("Only modes 0 and 1 defined") raise ValueError("Only modes 0 and 1 defined")
filename = f"test-MultiCorrelatedField-{mode}" + "-{}.png" filename = f"demo-MultiCorrelatedField-mode_{mode}" + "-{}.png"
# Two-dimensional signal on RGs with inhomogeneous exposure # Two-dimensional signal on RGs with inhomogeneous exposure
N = 128 N = 128
...@@ -76,47 +76,48 @@ if __name__ == '__main__': ...@@ -76,47 +76,48 @@ if __name__ == '__main__':
power_space_p = ift.PowerSpace(harmonic_space_p) power_space_p = ift.PowerSpace(harmonic_space_p)
power_space_e = ift.PowerSpace(harmonic_space_e) power_space_e = ift.PowerSpace(harmonic_space_e)
amp_p = ift.SLAmplitude( amp_p_dict = {
**{ 'target': power_space_p,
'target': power_space_p, 'n_pix': 64, # 64 spectral bins
'n_pix': 64, # 64 spectral bins 'keys': ['tau_p', 'phi_p', 'zm_p'],
'keys': ['tau_p', 'phi_p', 'zm_p'],
# Spectral smoothness (affects Gaussian process part)
# Spectral smoothness (affects Gaussian process part) 'a': 3, # relatively high variance of spectral curvature
'a': 3, # relatively high variance of spectral curvature 'k0': 0.4, # quefrency mode below which cepstrum flattens
'k0': 0.4, # quefrency mode below which cepstrum flattens
# Power-law part of spectrum
# Power-law part of spectrum 'sm': -4., # Expected exponent of power law
'sm': -4., # Expected exponent of power law 'sv': 1., # Prior standard deviation of -||-
'sv': 1., # Prior standard deviation of -||- 'im': 0., # Expected y-intercept of power law
'im': 0., # Expected y-intercept of power law 'iv': 0., # Prior standard deviation of -||-
'iv': 0., # Prior standard deviation of -||-
# zero-mode
# zero-mode 'za': 5., # IG alpha parameter
'za': 5., # IG alpha parameter 'zq': 0.001 # IG q parameter
'zq': 0.001 # IG q parameter }
})
amp_e_dict = {
amp_e = ift.SLAmplitude( 'target': power_space_e,
**{ 'n_pix': 64, # 64 spectral bins
'target': power_space_e, 'keys': ['tau_e', 'phi_e', 'zm_e'],
'n_pix': 64, # 64 spectral bins
'keys': ['tau_e', 'phi_e', 'zm_e'], # Spectral smoothness (affects Gaussian process part)
'a': 3, # relatively high variance of spectral curvature
# Spectral smoothness (affects Gaussian process part) 'k0': 0.4, # quefrency mode below which cepstrum flattens
'a': 3, # relatively high variance of spectral curvature
'k0': 0.4, # quefrency mode below which cepstrum flattens # Power-law part of spectrum
'sm': -3., # Expected exponent of power law
# Power-law part of spectrum 'sv': 1., # Prior standard deviation of -||-
'sm': -3., # Expected exponent of power law 'im': 0., # Expected y-intercept of power law
'sv': 1., # Prior standard deviation of -||- 'iv': 0., # Prior standard deviation of -||-
'im': 0., # Expected y-intercept of power law
'iv': 0., # Prior standard deviation of -||- # zero-mode
'za': 6., # IG alpha parameter
# zero-mode 'zq': 0.002 # IG q parameter
'za': 6., # IG alpha parameter }
'zq': 0.002 # IG q parameter
}) amp_p = ift.SLAmplitude(**amp_p_dict)
amp_e = ift.SLAmplitude(**amp_e_dict)
# Define sky operator # Define sky operator
# za=4, zq=5 leads to an amplitude median of one # za=4, zq=5 leads to an amplitude median of one
...@@ -125,10 +126,9 @@ if __name__ == '__main__': ...@@ -125,10 +126,9 @@ if __name__ == '__main__':
elif mode == 0: elif mode == 0:
amplitudes = [None, amp_e] amplitudes = [None, amp_e]
log_diffuse = ift.MultiCorrelatedField(sky_target, diff_zm_params = {'gs_a': 11., 'gs_q': 50.}
amplitudes, log_diffuse = ift.MultiCorrelatedField(sky_target, amplitudes,
gs_a=11., **diff_zm_params)
gs_q=50.)
sky = log_diffuse.exp() sky = log_diffuse.exp()
M = ift.DiagonalOperator(exposure) M = ift.DiagonalOperator(exposure)
...@@ -142,7 +142,6 @@ if __name__ == '__main__': ...@@ -142,7 +142,6 @@ if __name__ == '__main__':
data = lamb(mock_position) data = lamb(mock_position)
signal = sky(mock_position) signal = sky(mock_position)
# DEBUG # DEBUG
mock_log_diffuse_vals = log_diffuse(mock_position).to_global_data() mock_log_diffuse_vals = log_diffuse(mock_position).to_global_data()
print("std(mock_log_diffuse):", mock_log_diffuse_vals.std()) print("std(mock_log_diffuse):", mock_log_diffuse_vals.std())
...@@ -165,7 +164,7 @@ if __name__ == '__main__': ...@@ -165,7 +164,7 @@ if __name__ == '__main__':
iteration = 0 iteration = 0
# Minimize # Minimize
for i in range(3): for i in range(10):
KL = ift.MetricGaussianKL(pos, H, 10) KL = ift.MetricGaussianKL(pos, H, 10)
KL, convergence = minimizer(KL) KL, convergence = minimizer(KL)
pos = KL.position pos = KL.position
...@@ -182,7 +181,9 @@ if __name__ == '__main__': ...@@ -182,7 +181,9 @@ if __name__ == '__main__':
plot.add(cd(GR.adjoint(data)), title='Data') plot.add(cd(GR.adjoint(data)), title='Data')
plot.add(cd(reconst), title='Posterior Mean') plot.add(cd(reconst), title='Posterior Mean')
plot.add(cd(reconst - signal), title='Error') plot.add(cd(reconst - signal), title='Error')
plot.output(xsize=12, ysize=10, name=filename.format(f'reconstr-{iteration:02d}')) plot.output(xsize=12,
ysize=10,
name=filename.format(f'{iteration:02d}-reconstr'))
p_p = ift.PowerDistributor(harmonic_space_p, power_space_p) @ amp_p p_p = ift.PowerDistributor(harmonic_space_p, power_space_p) @ amp_p
p_e = ift.PowerDistributor(harmonic_space_e, power_space_e) @ amp_e p_e = ift.PowerDistributor(harmonic_space_e, power_space_e) @ amp_e
...@@ -197,13 +198,26 @@ if __name__ == '__main__': ...@@ -197,13 +198,26 @@ if __name__ == '__main__':
for sample in KL.samples: for sample in KL.samples:
sc.add(ift.power_analyze(p.force(pos + sample))) sc.add(ift.power_analyze(p.force(pos + sample)))
plot = ift.Plot() plot = ift.Plot()
plot.add([ plot.add(
sc.mean, [sc.mean, ift.power_analyze(p.force(mock_position))],
ift.power_analyze(p.force(mock_position)) linewidth=[1., 3.],
], label=['Posterior Mean', 'Ground Truth'],
linewidth=[1., 3.], title=f"Posterior {n} spectrum")
label=['Posterior Mean', 'Ground Truth'], plot.output(xsize=12,
title=f"Posterior {n} spectrum") ysize=10,
plot.output(xsize=12, ysize=10, name=filename.format(f'power_{n}-{iteration:02d}')) name=filename.format(f'{iteration:02d}-power_{n}'))
print("Saved results as '{}'.".format(filename.format('*'))) print("Saved results as '{}'.".format(filename.format('*')))
# Calculate InverseGamma prior misfits
keys = ['zm_e', 'zm_p', 'xi_pspec_scaling']
hyperparams = [(amp_e_dict['za'], amp_e_dict['zq']),
(amp_p_dict['za'], amp_p_dict['zq']),
(diff_zm_params['gs_a'], diff_zm_params['gs_q'])]
print("InverseGamma missfits:")
for key, hpar in zip(keys, hyperparams):
op = ift.InverseGammaOperator(sky.domain[key], *hpar)
print(f"{key}:")
print(f"\tmock: {float(op(mock_position[key]).local_data):+1.2e}")
print(f"\t rec: {float(op(pos[key]).local_data):+1.2e}\n")
...@@ -90,7 +90,7 @@ def testModelLibrary(space, seed): ...@@ -90,7 +90,7 @@ def testModelLibrary(space, seed):
np.random.seed(seed) np.random.seed(seed)
domain = ift.PowerSpace(space.get_default_codomain()) domain = ift.PowerSpace(space.get_default_codomain())
model = ift.SLAmplitude(target=domain, n_pix=4, a=.5, k0=2, sm=3, sv=1.5, model = ift.SLAmplitude(target=domain, n_pix=4, a=.5, k0=2, sm=3, sv=1.5,
im=1.75, iv=1.3) im=1.75, iv=1.3, za=5, zq=4)
assert_(isinstance(model, ift.Operator)) assert_(isinstance(model, ift.Operator))
S = ift.ScalingOperator(1., model.domain) S = ift.ScalingOperator(1., model.domain)
pos = S.draw_sample() pos = S.draw_sample()
...@@ -101,11 +101,11 @@ def testModelLibrary(space, seed): ...@@ -101,11 +101,11 @@ def testModelLibrary(space, seed):
pos = S.draw_sample() pos = S.draw_sample()
ift.extra.check_jacobian_consistency(model2, pos, ntries=20) ift.extra.check_jacobian_consistency(model2, pos, ntries=20)
#domtup = ift.DomainTuple.make((space, space)) domtup = ift.DomainTuple.make((space, space))
#model3 = ift.MfCorrelatedField(domtup, [model, model]) model3 = ift.MultiCorrelatedField(domtup, [model, model])
#S = ift.ScalingOperator(1., model3.domain) S = ift.ScalingOperator(1., model3.domain)
#pos = S.draw_sample() pos = S.draw_sample()
#ift.extra.check_jacobian_consistency(model3, pos, ntries=20) ift.extra.check_jacobian_consistency(model3, pos, ntries=20)
def testPointModel(space, seed): def testPointModel(space, seed):
......