...
 
Commits (2)
...@@ -129,7 +129,7 @@ if __name__ == '__main__': ...@@ -129,7 +129,7 @@ if __name__ == '__main__':
amplitudes = [None, amp_e] amplitudes = [None, amp_e]
# Diffuse global scaling prior parameters # Diffuse global scaling prior parameters
diff_gs_params = {'gs_mean': 5., 'gs_var': 2.5} diff_gs_params = {'gs_alpha': 0.5, 'gs_q': 2.0}
log_diffuse = ift.MultiCorrelatedField(sky_target, amplitudes, log_diffuse = ift.MultiCorrelatedField(sky_target, amplitudes,
**diff_gs_params) **diff_gs_params)
...@@ -218,15 +218,18 @@ if __name__ == '__main__': ...@@ -218,15 +218,18 @@ if __name__ == '__main__':
keys = ['zm_e', 'zm_p', 'xi_pspec_scaling'] keys = ['zm_e', 'zm_p', 'xi_pspec_scaling']
hyperparams = [(amp_e_dict['zm'], amp_e_dict['zv']), hyperparams = [(amp_e_dict['zm'], amp_e_dict['zv']),
(amp_p_dict['zm'], amp_p_dict['zv']), (amp_p_dict['zm'], amp_p_dict['zv']),
(diff_gs_params['gs_mean'], diff_gs_params['gs_var'])] (diff_gs_params['gs_alpha'], diff_gs_params['gs_q'])]
elif mode == 0: elif mode == 0:
keys = ['zm_e', 'xi_pspec_scaling'] keys = ['zm_e', 'xi_pspec_scaling']
hyperparams = [(amp_e_dict['zm'], amp_e_dict['zv']), hyperparams = [(amp_e_dict['zm'], amp_e_dict['zv']),
(diff_gs_params['gs_mean'], diff_gs_params['gs_var'])] (diff_gs_params['gs_alpha'], diff_gs_params['gs_q'])]
print("InverseGamma fits:") print("InverseGamma fits:")
for key, hpar in zip(keys, hyperparams): for key, hpar in zip(keys, hyperparams):
op = InvGammaFromMeanVar(sky.domain[key], *hpar) if key in ['zm_e', 'zm_p']:
op = InvGammaFromMeanVar(sky.domain[key], *hpar)
else:
op = ift.InverseGammaOperator(sky.domain[key], *hpar)
print(f"{key}:") print(f"{key}:")
print(f"\tmock: {float(op(mock_position[key]).local_data):+1.2e}") 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") print(f"\t rec: {float(op(pos[key]).local_data):+1.2e}\n")
...@@ -82,7 +82,7 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None): ...@@ -82,7 +82,7 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
return ht(vol * A * ducktape(h_space, None, name)) return ht(vol * A * ducktape(h_space, None, name))
def MultiCorrelatedField(target, amplitudes, gs_mean, gs_var, name='xi'): def MultiCorrelatedField(target, amplitudes, gs_alpha, gs_q, name='xi'):
"""Constructs an operator which turns white Gaussian excitation fields """Constructs an operator which turns white Gaussian excitation fields
into a (partially) correlated field defined on a DomainTuple with separate into a (partially) correlated field defined on a DomainTuple with separate
correlation structures per subdomain. correlation structures per subdomain.
...@@ -111,14 +111,12 @@ def MultiCorrelatedField(target, amplitudes, gs_mean, gs_var, name='xi'): ...@@ -111,14 +111,12 @@ def MultiCorrelatedField(target, amplitudes, gs_mean, gs_var, name='xi'):
pass `None` for this subdomain. pass `None` for this subdomain.
name : string name : string
:class:`MultiField` key for xi-field. :class:`MultiField` key for xi-field.
gs_mean : float gs_alpha : float
Prior mean of global_scaling**2. `alpha` parameter of the global_scaling prior
Used to calculate appropriate parameters for See :class:`InverseGammaOperator` for interpretation.
a :class:`InverseGammaOperator`. gs_q : float
gs_var : float `q` parameter of the global_scaling prior
Prior variance of global_scaling**2. See :class:`InverseGammaOperator` for interpretation.
Used to calculate appropriate parameters for
a :class:`InverseGammaOperator`.
Returns Returns
------- -------
...@@ -171,13 +169,14 @@ def MultiCorrelatedField(target, amplitudes, gs_mean, gs_var, name='xi'): ...@@ -171,13 +169,14 @@ def MultiCorrelatedField(target, amplitudes, gs_mean, gs_var, name='xi'):
# Normalize amplitudes and apply spreader stack # Normalize amplitudes and apply spreader stack
if amplitudes[i] is not None: if amplitudes[i] is not None:
# Dimension with correlation structure # Dimension with correlation structure
normed_amplitude = OperatorNormalizer(amplitudes[i]) normed_amplitude = OperatorNormalizer(amplitudes[i], p=2)
a_i.append(spr @ normed_amplitude) a_i.append(spr @ normed_amplitude)
else: else:
# No correlation structure in this dimension # No correlation structure in this dimension
# -> constant amplitude # -> constant amplitude
# still got to normalize, though. # Still got to normalize, tough:
amplitude_integral = Field.full(spr.domain, 1.).integrate() # Since the field contains only ones, 2-norm = sqrt(1-norm)
amplitude_integral = Field.full(spr.domain, 1.).integrate().sqrt()
fct /= amplitude_integral fct /= amplitude_integral
fct = ScalingOperator(fct, pd.domain) fct = ScalingOperator(fct, pd.domain)
...@@ -193,9 +192,6 @@ def MultiCorrelatedField(target, amplitudes, gs_mean, gs_var, name='xi'): ...@@ -193,9 +192,6 @@ def MultiCorrelatedField(target, amplitudes, gs_mean, gs_var, name='xi'):
# Generate global power spectrum scaling operator # Generate global power spectrum scaling operator
vdot_ig = VdotOperator(Field.full(A.target, 1.)) vdot_ig = VdotOperator(Field.full(A.target, 1.))
gs_base = gs_mean**2 / gs_var + 1.
gs_alpha = gs_base + 1.
gs_q = gs_mean * gs_base
inv_gamma_op = InverseGammaOperator(vdot_ig.adjoint.domain, gs_alpha, gs_q) inv_gamma_op = InverseGammaOperator(vdot_ig.adjoint.domain, gs_alpha, gs_q)
global_amplitude = vdot_ig.adjoint @ inv_gamma_op.sqrt() global_amplitude = vdot_ig.adjoint @ inv_gamma_op.sqrt()
......
...@@ -19,25 +19,27 @@ from ..field import Field ...@@ -19,25 +19,27 @@ from ..field import Field
from .simple_linear_operators import VdotOperator, WeightApplier from .simple_linear_operators import VdotOperator, WeightApplier
def OperatorNormalizer(op): def OperatorNormalizer(op, p=1):
"""Wraps `op` into a shell that ensures the fields returned by """Wraps `op` into a shell that ensures the fields returned by
(shell @ op) will have an integral of one. (shell @ op) will have a p-norm of one.
Implements an operator analogously to Implements an operator analogously to
OperatorNormalizer(op) = op / op.integrate() OperatorNormalizer(op) = op / |op|_p
where the integration is implemented as an operator chain. where the integration is implemented as an operator chain.
Parameters Parameters
---------- ----------
op : Operator op : Operator
p : int
Which p-norm to normalize under
""" """
ones_field = Field.full(op.target, 1.) ones_field = Field.full(op.target, 1.)
vdot = VdotOperator(ones_field) vdot = VdotOperator(ones_field)
wa = WeightApplier(op.target, spaces=None, power=1) wa = WeightApplier(op.target, spaces=None, power=1)
# integrate as an operator # integrate as an operator
op_integrated = vdot.adjoint @ vdot @ wa @ op op_integrated_inv = vdot.adjoint @ (vdot @ wa @ (op)**p)**(-1/p)
return op * op_integrated**(-1) return op * op_integrated_inv
...@@ -103,8 +103,8 @@ def testModelLibrary(space, seed): ...@@ -103,8 +103,8 @@ def testModelLibrary(space, seed):
domtup = ift.DomainTuple.make((space, space)) domtup = ift.DomainTuple.make((space, space))
model3 = ift.MultiCorrelatedField(domtup, [model, model], model3 = ift.MultiCorrelatedField(domtup, [model, model],
gs_mean=5, gs_alpha=0.5,
gs_var=2.5) gs_q=2.0)
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)
......