Commit 010555de authored by Lukas Platz's avatar Lukas Platz

inverse gamma on zero mode of SLAmplitude

parent eaaed9cc
Pipeline #61491 failed with stages
in 4 minutes and 11 seconds
......@@ -25,7 +25,9 @@ from ..operators.exp_transform import ExpTransform
from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
from ..operators.simple_linear_operators import VdotOperator
from ..sugar import makeOp
from ..library.inverse_gamma_operator import InverseGammaOperator
def _ceps_kernel(k, a, k0):
......@@ -85,7 +87,7 @@ def CepstrumOperator(target, a, k0):
k0 = np.array(k0)
if len(k0) != len(target.shape):
raise ValueError
if np.any(np.array(k0) <= 0):
if np.any(k0 <= 0):
raise ValueError
qht = QHTOperator(target)
......@@ -112,7 +114,74 @@ def CepstrumOperator(target, a, k0):
return sym @ qht @ makeOp(cepstrum.sqrt())
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
def OneDCepstrumOperator(target, a, k0):
"""Turns a white Gaussian random field into a smooth field on a LogRGSpace.
The smooth field is composed out of three operators:
sym @ qht @ diag(sqrt_ceps),
where sym is a :class:`SymmetrizingOperator`, qht is a :class:`QHTOperator`
and ceps is the so-called cepstrum:
.. math::
\\mathrm{sqrt\\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
These operators are combined in this fashion in order to generate:
- A field which is smooth, i.e. second derivatives are punished (note
that the sqrt-cepstrum is essentially proportional to 1/k**2).
- A field which is symmetric around the pixel in the middle of the space.
This is result of the :class:`SymmetrizingOperator` and needed in order
to decouple the degrees of freedom at the beginning and the end of the
amplitude whenever :class:`CepstrumOperator` is used as in
:class:`SLAmplitude`.
Parameters
----------
target : LogRGSpace
Target domain of the operator, needs to be non-harmonic and 1d.
a : float
Cutoff of smoothness prior (positive only). Controls the
regularization of the inverse laplace operator to be finite at zero.
Larger values for the cutoff results in a weaker constraining prior.
k0 : float
Strength of smoothness prior in quefrency space (positive only) along
each axis. If float then the strength is the same along each axis.
Larger values result in a weaker constraining prior.
"""
a = float(a)
if a <= 0:
raise ValueError
target = DomainTuple.make(target)
if len(target.shape) != 1 or target[0].harmonic:
raise TypeError
if not isinstance(k0, (float, int)):
raise ValueError
if k0 <= 0:
raise ValueError
qht = QHTOperator(target)
sym = SymmetrizingOperator(target)
# Compute cepstrum field
dom = qht.domain[0]
shape = dom.shape
q_array = dom.get_k_array()
# Fill all non-zero modes
no_zero_modes = (slice(1, None),)
ks = q_array[(slice(None),) + no_zero_modes]
cepstrum_field = np.zeros(shape)
cepstrum_field[no_zero_modes] = _ceps_kernel(ks, a, k0)
# zero mode stays zero
cepstrum = Field.from_global_data(dom, cepstrum_field)
return sym @ qht @ makeOp(cepstrum.sqrt())
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
keys=['tau', 'phi', 'zm']):
'''Operator for parametrizing smooth amplitudes (square roots of power
spectra).
......@@ -130,6 +199,7 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
the output of the generated operator is:
AmplitudeOperator = 0.5*(smooth_component + linear_component)
+ zero_mode
This is then exponentiated and exponentially binned (in this order).
......@@ -161,6 +231,10 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
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`
Returns
-------
......@@ -169,12 +243,13 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
which returns on its target a power spectrum which consists out of a
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, keys=keys).exp()
return LinearSLAmplitude(target=target, n_pix=n_pix, a=a, k0=k0,
sm=sm, sv=sv, im=im, iv=iv, za=za, zq=zq,
keys=keys).exp()
def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
keys=['tau', 'phi']):
def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
keys=['tau', 'phi', 'zm']):
'''LinearOperator for parametrizing smooth log-amplitudes (square roots of
power spectra).
......@@ -188,6 +263,9 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
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:
raise ValueError
et = ExpTransform(target, n_pix)
dom = et.domain[0]
......@@ -202,10 +280,20 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
# Smooth deviations from linear component
dct = {'a': a, 'k0': k0}
smooth = CepstrumOperator(dom, **dct).ducktape(keys[0])
# Combine linear and smooth component
loglog_ampl = 0.5*(smooth + linear)
zm_mask = np.ones(dom.shape)
zm_mask[(0,) * len(dom.shape)] = 0.
zm_mask = Field.from_global_data(dom, zm_mask)
zm_mask_inv = Field.full(dom, 1.) - zm_mask
zm_mask = makeOp(zm_mask)
smooth = zm_mask @ CepstrumOperator(dom, **dct).ducktape(keys[0])
# Zero Mode
vdot = VdotOperator(zm_mask_inv)
ig = InverseGammaOperator(vdot.target, za, zq)
zero_mode = vdot.adjoint @ ig.ducktape(keys[2])
# Combine components
loglog_ampl = 0.5*(smooth + linear) + zero_mode
# Go from loglog-space to linear-linear-space
return et @ loglog_ampl
# ToDo Development
## In `SmoothLinearAmplitude`:
- exchange zero mode for Inverse-Gamma distributed random variable
## In `MultiCorrelatedField`:
- before outer product devide by integral
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment