...
 
Commits (2)
......@@ -26,6 +26,11 @@ import numpy as np
import nifty5 as ift
def InvGammaFromMeanVar(domain, mean, variance):
base = mean**2 / variance + 1.
return ift.InverseGammaOperator(domain, alpha=base + 1., q=mean * base)
def exposure_2d():
# Structured exposure for 2D mode
x_shape, y_shape = sky_target.shape
......@@ -88,9 +93,9 @@ if __name__ == '__main__':
'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
# zero-mode prior
'zm': 2.5e-4, # Prior mean
'zv': 1.0e-4 # Prior variance
}
amp_e_dict = {
......@@ -108,9 +113,9 @@ if __name__ == '__main__':
'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
# zero-mode prior
'zm': 3.0e-3, # Prior mean
'zv': 2.0e-3 # Prior variance
}
amp_p = ift.SLAmplitude(**amp_p_dict)
......@@ -124,7 +129,7 @@ if __name__ == '__main__':
amplitudes = [None, amp_e]
# Diffuse global scaling prior parameters
diff_gs_params = {'gs_a': 11., 'gs_q': 50.}
diff_gs_params = {'gs_mean': 5., 'gs_var': 2.5}
log_diffuse = ift.MultiCorrelatedField(sky_target, amplitudes,
**diff_gs_params)
......@@ -211,17 +216,17 @@ if __name__ == '__main__':
# 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'])]
hyperparams = [(amp_e_dict['zm'], amp_e_dict['zv']),
(amp_p_dict['zm'], amp_p_dict['zv']),
(diff_gs_params['gs_mean'], diff_gs_params['gs_var'])]
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'])]
hyperparams = [(amp_e_dict['zm'], amp_e_dict['zv']),
(diff_gs_params['gs_mean'], diff_gs_params['gs_var'])]
print("InverseGamma fits:")
for key, hpar in zip(keys, hyperparams):
op = ift.InverseGammaOperator(sky.domain[key], *hpar)
op = InvGammaFromMeanVar(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")
......@@ -75,9 +75,9 @@ if __name__ == '__main__':
'im': 0, # y-intercept mean, in-/decrease for more/less contrast
'iv': .3, # y-intercept variance
# Zero mode prior
'za': 5./2, # alpha of InverseGammaOperator
'zq': 3./2 # q of InverseGammaOperator
# Zero mode prior (InverseGammaOperator)
'zm': 2., # Prior mean
'zq': 0.5 # Prior variance
}
A = ift.SLAmplitude(**dct)
......
......@@ -82,7 +82,7 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
return ht(vol * A * ducktape(h_space, None, name))
def MultiCorrelatedField(target, amplitudes, gs_a, gs_q, name='xi'):
def MultiCorrelatedField(target, amplitudes, gs_mean, gs_var, name='xi'):
"""Constructs an operator which turns white Gaussian excitation fields
into a (partially) correlated field defined on a DomainTuple with separate
correlation structures per subdomain.
......@@ -111,12 +111,14 @@ def MultiCorrelatedField(target, amplitudes, gs_a, gs_q, name='xi'):
pass `None` for this subdomain.
name : string
:class:`MultiField` key for xi-field.
gs_a : float
`alpha` parameter of the global scaling prior.
See :class:`InverseGammaOperator` for interpretaition.
gs_q : float
`q` parameter of the global scaling prior.
See :class:`InverseGammaOperator` for interpretaition.
gs_mean : float
Prior mean of global_scaling**2.
Used to calculate appropriate parameters for
a :class:`InverseGammaOperator`.
gs_var : float
Prior variance of global_scaling**2.
Used to calculate appropriate parameters for
a :class:`InverseGammaOperator`.
Returns
-------
......@@ -191,7 +193,10 @@ def MultiCorrelatedField(target, amplitudes, gs_a, gs_q, name='xi'):
# Generate global power spectrum scaling operator
vdot_ig = VdotOperator(Field.full(A.target, 1.))
inv_gamma_op = InverseGammaOperator(vdot_ig.adjoint.domain, gs_a, gs_q)
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)
global_amplitude = vdot_ig.adjoint @ inv_gamma_op.sqrt()
A_times_xi = vol * A * ducktape(hsp, None, name) * \
......
......@@ -180,7 +180,7 @@ def OneDCepstrumOperator(target, a, k0):
return sym @ qht @ makeOp(cepstrum.sqrt())
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, zm, zv,
keys=['tau', 'phi', 'zm']):
'''Operator for parametrizing smooth amplitudes (square roots of power
spectra).
......@@ -231,10 +231,10 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
LogRGSpace (see :class:`ExpTransform`).
iv : float
Prior standard deviation of y-intercept of power law.
za : float
zero-mode prior parameter `alpha`
zq : float
zero-mode prior parameter `q`
zm : float
Prior mean value of the zero mode.
zv : float
Prior variance of the zero mode.
Returns
-------
......@@ -244,11 +244,11 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
smooth and a linear part.
'''
return LinearSLAmplitude(target=target, n_pix=n_pix, a=a, k0=k0,
sm=sm, sv=sv, im=im, iv=iv, za=za, zq=zq,
sm=sm, sv=sv, im=im, iv=iv, zm=zm, zv=zv,
keys=keys).exp()
def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, zm, zv,
keys=['tau', 'phi', 'zm']):
'''LinearOperator for parametrizing smooth log-amplitudes (square roots of
power spectra).
......@@ -263,8 +263,8 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
sm, sv, im, iv = float(sm), float(sv), float(im), float(iv)
if sv <= 0 or iv < 0:
raise ValueError
za, zq = float(za), float(zq)
if za <= 0 or zq <= 0:
zm, zv = float(zm), float(zv)
if zm <= 0 or zv <= 0:
raise ValueError
et = ExpTransform(target, n_pix)
......@@ -289,7 +289,8 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
# Zero Mode
vdot = VdotOperator(zm_mask_inv)
ig = InverseGammaOperator(vdot.target, za, zq)
zb = zm**2 / zv + 1.
ig = InverseGammaOperator(vdot.target, alpha=zb+1., q=zm*zb)
zero_mode = vdot.adjoint @ ig.ducktape(keys[2])
# Combine components
......
......@@ -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, za=5, zq=4)
im=1.75, iv=1.3, zm=1., zv=0.5)
assert_(isinstance(model, ift.Operator))
S = ift.ScalingOperator(1., model.domain)
pos = S.draw_sample()
......@@ -102,7 +102,9 @@ def testModelLibrary(space, seed):
ift.extra.check_jacobian_consistency(model2, pos, ntries=20)
domtup = ift.DomainTuple.make((space, space))
model3 = ift.MultiCorrelatedField(domtup, [model, model], gs_a=5, gs_q=4)
model3 = ift.MultiCorrelatedField(domtup, [model, model],
gs_mean=5,
gs_var=2.5)
S = ift.ScalingOperator(1., model3.domain)
pos = S.draw_sample()
ift.extra.check_jacobian_consistency(model3, pos, ntries=20)
......