...
 
Commits (2)
......@@ -64,10 +64,7 @@ if __name__ == '__main__':
sky_target = ift.DomainTuple.make((position_space, energy_space))
plot_domain = ift.RGSpace((N, N))
exposure = exposure_2d()
if mode == 0:
# need a much higher exposure be able to reconstruct the spectrum
exposure = exposure * 100
exposure = 10. * exposure_2d()
# Define harmonic space and harmonic transform
harmonic_space_p = position_space.get_default_codomain()
......@@ -126,9 +123,11 @@ if __name__ == '__main__':
elif mode == 0:
amplitudes = [None, amp_e]
diff_zm_params = {'gs_a': 11., 'gs_q': 50.}
# Diffuse global scaling prior parameters
diff_gs_params = {'gs_a': 11., 'gs_q': 50.}
log_diffuse = ift.MultiCorrelatedField(sky_target, amplitudes,
**diff_zm_params)
**diff_gs_params)
sky = log_diffuse.exp()
M = ift.DiagonalOperator(exposure)
......@@ -209,13 +208,18 @@ if __name__ == '__main__':
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 prior fits
if mode == 1:
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_gs_params['gs_a'], diff_gs_params['gs_q'])]
elif mode == 0:
keys = ['zm_e', 'xi_pspec_scaling']
hyperparams = [(amp_e_dict['za'], amp_e_dict['zq']),
(diff_gs_params['gs_a'], diff_gs_params['gs_q'])]
print("InverseGamma missfits:")
print("InverseGamma fits:")
for key, hpar in zip(keys, hyperparams):
op = ift.InverseGammaOperator(sky.domain[key], *hpar)
print(f"{key}:")
......
......@@ -155,6 +155,7 @@ def MultiCorrelatedField(target, amplitudes, gs_a, gs_q, name='xi'):
# Spread amplitudes into their respective constant dimensions
# and calculate their outer product
a_i = []
fct = 1.
for i in range(n_subdoms):
# Assemble spreader stack for this amplitude
......@@ -169,16 +170,16 @@ def MultiCorrelatedField(target, amplitudes, gs_a, gs_q, name='xi'):
if amplitudes[i] is not None:
# Dimension with correlation structure
normed_amplitude = OperatorNormalizer(amplitudes[i])
a_i.append(spr @ normed_amplitude)
else:
# No correlation structure in this dimension
# -> constant amplitude
# still got to normalize, though.
amplitude_integral = Field.full(spr.domain, 1.).integrate()
normed_amplitude = ScalingOperator(spr.domain,
1. / amplitude_integral)
a_i.append(spr @ normed_amplitude)
fct /= amplitude_integral
a = reduce(mul, a_i)
fct = ScalingOperator(fct, pd.domain)
a = fct @ reduce(mul, a_i)
# Assemble global power spectrum
A = pd @ a
......@@ -193,7 +194,6 @@ def MultiCorrelatedField(target, amplitudes, gs_a, gs_q, name='xi'):
inv_gamma_op = InverseGammaOperator(vdot_ig.adjoint.domain, gs_a, gs_q)
global_amplitude = vdot_ig.adjoint @ inv_gamma_op.sqrt()
#A_times_xi = vol * A * ducktape(hsp, None, name)
A_times_xi = vol * A * ducktape(hsp, None, name) * \
global_amplitude.ducktape(name + '_pspec_scaling')
......