Commit 79dacba1 authored by Ruestig, Julian (jruestig)'s avatar Ruestig, Julian (jruestig) 📡

Smooth Linear Amplitude takes zero mode now

parent e6eb8d06
......@@ -17,21 +17,22 @@
from functools import reduce
from operator import mul
import numpy as np
from ..domain_tuple import DomainTuple
from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import ducktape
from ..operators.value_inserter import ValueInserter
from ..operators.simple_linear_operators import FieldAdapter
from ..operators.diagonal_operator import DiagonalOperator
from ..sugar import from_global_data
from ..library.inverse_gamma_operator import InverseGammaOperator
# import numpy as np
# from ..operators.value_inserter import ValueInserter
# from ..operators.simple_linear_operators import FieldAdapter
# from ..operators.diagonal_operator import DiagonalOperator
# from ..sugar import from_global_data
# from ..library.inverse_gamma_operator import InverseGammaOperator
def CorrelatedField(target, amplitude_operator, q=1e-5, name='xi', codomain=None):
def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
"""Constructs an operator which turns a white Gaussian excitation field
into a correlated field.
......@@ -76,24 +77,14 @@ def CorrelatedField(target, amplitude_operator, q=1e-5, name='xi', codomain=None
A = power_distributor(amplitude_operator)
vol = h_space.scalar_dvol**-0.5
zero_mode = ValueInserter(A.target, (0,)*len(A.target.shape))
mask = np.ones(A.target.shape)
mask[(0,)*len(A.target.shape)] = 0.
mask = from_global_data(A.target, mask)
mask = DiagonalOperator(mask)
zero_mode = zero_mode @ InverseGammaOperator(
zero_mode.domain, 0.5, q) @ FieldAdapter(
zero_mode.domain, "zero_mode")
# When doubling the resolution of `tgt` the value of the highest k-mode
# will scale with a square root. `vol` cancels this effect such that the
# same power spectrum can be used for the spaces with the same volume,
# different resolutions and the same object in them.
return ht(vol*(mask@A+zero_mode)*ducktape(h_space, None, name))
return ht(vol*(A)*ducktape(h_space, None, name))
def MfCorrelatedField(target, amplitudes, q=1e-14, name='xi'):
def MfCorrelatedField(target, amplitudes, name='xi'):
"""Constructs an operator which turns white Gaussian excitation fields
into a correlated field defined on a DomainTuple with two entries and two
separate correlation structures.
......@@ -149,14 +140,4 @@ def MfCorrelatedField(target, amplitudes, q=1e-14, name='xi'):
# For `vol` see comment in `CorrelatedField`
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
zero_mode = ValueInserter(A.target, (0,)*len(A.target.shape))
mask = np.ones(A.target.shape)
mask[(0,)*len(A.target.shape)] = 0.
mask = from_global_data(A.target, mask)
mask = DiagonalOperator(mask)
zero_mode = zero_mode @ InverseGammaOperator(
zero_mode.domain, 0.5, q) @ FieldAdapter(
zero_mode.domain, "zero_mode")
return ht(vol*(mask @ A + zero_mode)*ducktape(hsp, None, name))
return ht(vol*A*ducktape(hsp, None, name))
......@@ -27,6 +27,12 @@ from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
from ..sugar import makeOp
from ..library.inverse_gamma_operator import InverseGammaOperator
from ..sugar import from_global_data
from ..operators.diagonal_operator import DiagonalOperator
from ..operators.simple_linear_operators import FieldAdapter
from ..operators.value_inserter import ValueInserter
def _ceps_kernel(k, a, k0):
return (a/(1 + np.sum((k.T/k0)**2, axis=-1).T))**2
......@@ -112,7 +118,8 @@ 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 SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
keys=['tau', 'phi', 'zero_mode']):
'''Operator for parametrizing smooth amplitudes (square roots of power
spectra).
......@@ -161,6 +168,11 @@ 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
The alpha-parameter of the inverse-gamma distribution.
Setting the a seperate prior on the zeroGmode of the amplitude model.
zq : float
The q-parameter of the inverse-gamma distribution.
Returns
-------
......@@ -174,6 +186,7 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
a, k0 = float(a), float(k0)
sm, sv, im, iv = float(sm), float(sv), float(im), float(iv)
za, zq = float(za), float(zq)
if sv <= 0 or iv <= 0:
raise ValueError
......@@ -195,5 +208,15 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
# Combine linear and smooth component
loglog_ampl = 0.5*(smooth + linear)
zero_mode = ValueInserter(et.target, (0,)*len(et.target.shape))
mask = np.ones(et.target.shape)
mask[(0,)*len(et.target.shape)] = 0.
mask = from_global_data(et.target, mask)
mask = DiagonalOperator(mask)
zero_mode = zero_mode @ InverseGammaOperator(
zero_mode.domain, za, zq) @ FieldAdapter(
zero_mode.domain, keys[2])
# Go from loglog-space to linear-linear-space
return et @ loglog_ampl.exp()
return (mask @ (et @ loglog_ampl.exp()) + zero_mode)
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