...
 
Commits (4)
......@@ -55,7 +55,7 @@ if __name__ == '__main__':
if not mode in (0, 1):
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
N = 128
......@@ -76,47 +76,48 @@ if __name__ == '__main__':
power_space_p = ift.PowerSpace(harmonic_space_p)
power_space_e = ift.PowerSpace(harmonic_space_e)
amp_p = ift.SLAmplitude(
**{
'target': power_space_p,
'n_pix': 64, # 64 spectral bins
'keys': ['tau_p', 'phi_p', 'zm_p'],
# Spectral smoothness (affects Gaussian process part)
'a': 3, # relatively high variance of spectral curvature
'k0': 0.4, # quefrency mode below which cepstrum flattens
# Power-law part of spectrum
'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 -||-
# 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', 'zm_e'],
# Spectral smoothness (affects Gaussian process part)
'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
'sv': 1., # 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
'zq': 0.002 # IG q parameter
})
amp_p_dict = {
'target': power_space_p,
'n_pix': 64, # 64 spectral bins
'keys': ['tau_p', 'phi_p', 'zm_p'],
# Spectral smoothness (affects Gaussian process part)
'a': 3, # relatively high variance of spectral curvature
'k0': 0.4, # quefrency mode below which cepstrum flattens
# Power-law part of spectrum
'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 -||-
# zero-mode
'za': 5., # IG alpha parameter
'zq': 0.001 # IG q parameter
}
amp_e_dict = {
'target': power_space_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
'k0': 0.4, # quefrency mode below which cepstrum flattens
# Power-law part of spectrum
'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 -||-
# zero-mode
'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
# za=4, zq=5 leads to an amplitude median of one
......@@ -125,10 +126,9 @@ if __name__ == '__main__':
elif mode == 0:
amplitudes = [None, amp_e]
log_diffuse = ift.MultiCorrelatedField(sky_target,
amplitudes,
gs_a=11.,
gs_q=50.)
diff_zm_params = {'gs_a': 11., 'gs_q': 50.}
log_diffuse = ift.MultiCorrelatedField(sky_target, amplitudes,
**diff_zm_params)
sky = log_diffuse.exp()
M = ift.DiagonalOperator(exposure)
......@@ -142,7 +142,6 @@ if __name__ == '__main__':
data = lamb(mock_position)
signal = sky(mock_position)
# DEBUG
mock_log_diffuse_vals = log_diffuse(mock_position).to_global_data()
print("std(mock_log_diffuse):", mock_log_diffuse_vals.std())
......@@ -165,7 +164,7 @@ if __name__ == '__main__':
iteration = 0
# Minimize
for i in range(3):
for i in range(10):
KL = ift.MetricGaussianKL(pos, H, 10)
KL, convergence = minimizer(KL)
pos = KL.position
......@@ -182,7 +181,9 @@ if __name__ == '__main__':
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}'))
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_e = ift.PowerDistributor(harmonic_space_e, power_space_e) @ amp_e
......@@ -197,13 +198,26 @@ if __name__ == '__main__':
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}'))
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'{iteration:02d}-power_{n}'))
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):
np.random.seed(seed)
domain = ift.PowerSpace(space.get_default_codomain())
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))
S = ift.ScalingOperator(1., model.domain)
pos = S.draw_sample()
......@@ -101,11 +101,11 @@ def testModelLibrary(space, seed):
pos = S.draw_sample()
ift.extra.check_jacobian_consistency(model2, pos, ntries=20)
#domtup = ift.DomainTuple.make((space, space))
#model3 = ift.MfCorrelatedField(domtup, [model, model])
#S = ift.ScalingOperator(1., model3.domain)
#pos = S.draw_sample()
#ift.extra.check_jacobian_consistency(model3, pos, ntries=20)
domtup = ift.DomainTuple.make((space, space))
model3 = ift.MultiCorrelatedField(domtup, [model, model])
S = ift.ScalingOperator(1., model3.domain)
pos = S.draw_sample()
ift.extra.check_jacobian_consistency(model3, pos, ntries=20)
def testPointModel(space, seed):
......