Commit 583ae212 authored by Lukas Platz's avatar Lukas Platz

MultiCorrelatedField with new ansatz

parent e1e4be1b
......@@ -27,6 +27,7 @@ from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import ducktape, VdotOperator
from ..operators.scaling_operator import ScalingOperator
from ..operators.contraction_operator import ContractionOperator
from ..operators.operator_normalizer import OperatorNormalizer
from ..library.inverse_gamma_operator import InverseGammaOperator
......@@ -81,36 +82,40 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
return ht(vol * A * ducktape(h_space, None, name))
def MultiCorrelatedField(target, amplitudes, va, vq, name='xi'):
def MultiCorrelatedField(target, amplitudes, gs_a, gs_q, 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.
The signal variance is set by an InverseGamma prior parametrized
by `va` and `vq`.
Uses of this operator include modeling diffuse emissions with correlation
structure in the spatial and energy dimensions and point sources with
correlation structure in the energy and time dimensions, but not in the
spatial dimension. Generalizations to arbitray combinations of domains
with and without correlation structure are possible.
The field is defined as
`MultiCorrelatedField` = HT( A * xi * volume_factor * global_scaling )
and takes a list of Amplitude operators from which to construct A
and the parameters of the global scaling prior, provided by a
:class:`InverseGammaOperator`.
Parameters
----------
target : Domain, DomainTuple or tuple of Domain
Target of the operator.
amplitudes: Iterable of Operator
amplitudes: List of Operators
List of amplitude operators. Must contain one per target subdomain.
If a subdomain without correlations is wanted,
pass `None` for this subdomain.
name : string
:class:`MultiField` key for xi-field.
va : float
`alpha` parameter of the variance scaling prior.
See :class:`InverseGammaOperator` for interpretaition.
vq : float
`q` parameter of the variance scaling prior.
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.
Returns
......@@ -136,11 +141,11 @@ def MultiCorrelatedField(target, amplitudes, va, vq, name='xi'):
hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
# Assemble PowerDistributor stack
def gen_pd(domain, i):
def gen_pd(pd_target, i):
if amplitudes[i] is not None:
return PowerDistributor(domain, amplitudes[i].target[0], i)
return PowerDistributor(pd_target, amplitudes[i].target[0], i)
else:
return PowerDistributor(domain, PowerSpace(domain[i]), i)
return PowerDistributor(pd_target, PowerSpace(pd_target[i]), i)
pd = gen_pd(hsp, 0)
for i in range(1, n_subdoms - 1):
......@@ -150,10 +155,8 @@ def MultiCorrelatedField(target, amplitudes, va, vq, name='xi'):
# Spread amplitudes into their respective constant dimensions
# and calculate their outer product
a_i = []
for i in range(n_subdoms):
if amplitudes[i] is None:
# No correlation structure in this dimension
continue
# Assemble spreader stack for this amplitude
spr = ScalingOperator(1., pd.domain)
for j in range(n_subdoms):
......@@ -161,34 +164,38 @@ def MultiCorrelatedField(target, amplitudes, va, vq, name='xi'):
spr = spr @ ContractionOperator(spr.domain, 0).adjoint
if j > i:
spr = spr @ ContractionOperator(spr.domain, 1).adjoint
a_i.append(spr @ amplitudes[i])
a = reduce(mul, a_i)
# Assemble global power spectrum
p = pd @ a
# Normalize amplitudes and apply spreader stack
if amplitudes[i] is not None:
# Dimension with correlation structure
normed_amplitude = OperatorNormalizer(amplitudes[i])
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)
# Fix signal variance to one by dividing the global power spectrum
# by its integral. The volume factor of the integration is
# included later on.
vdot = VdotOperator(Field.full(p.target, 1.))
normed_pspec = p * ((vdot.adjoint @ vdot)(p))**(-1)
a = reduce(mul, a_i)
# Multiply by phenomenological factor to get to
# a signal variance of one
normed_pspec = normed_pspec * 10.
# Assemble global power spectrum
A = pd @ a
# Volume factors
# One (pixel volume)**(-0.5) included, for explanaition see
# respective comment in `CorrelatedField`.
# From the normalization, there is an additional factor of
# (pixel volume)**(-1.0)
vol = reduce(mul, [sp.scalar_dvol**-1.5 for sp in hsp])
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
# Generate global power spectrum scaling operator
# Used to control the signal variance
vdot_ig = VdotOperator(Field.full(p.target, 1.))
inv_gamma_op = InverseGammaOperator(vdot_ig.adjoint.domain, va, vq)
global_amplitude = vdot_ig.adjoint @ inv_gamma_op.log()
vdot_ig = VdotOperator(Field.full(A.target, 1.))
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')
# Assemble global HT
ht = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
......@@ -200,8 +207,4 @@ def MultiCorrelatedField(target, amplitudes, va, vq, name='xi'):
space=n_subdoms - 1)
ht = htn @ ht
A_times_xi = vol * normed_pspec * ducktape(hsp, None, name)
#A_times_xi = vol * (normed_pspec * ducktape(hsp, None, name)) * \
# global_amplitude.ducktape(name + '_pspec_scaling')
return ht(A_times_xi)
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