Scheduled maintenance on Monday 2019-06-24 between 10:00-11:00 CEST

Value Insert InverseGamma at 0mode in MFCorrelatedField

parent eba1fd40
Pipeline #44103 failed with stages
in 8 minutes and 43 seconds
......@@ -93,10 +93,12 @@ if __name__ == '__main__':
# Apply the mapping
prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
print(prior_correlation_structure)
# Insert the result into the diagonal of an harmonic space operator
S = ift.DiagonalOperator(prior_correlation_structure)
# S is the prior field covariance
assert False
# Build instrument response consisting of a discretization, mask
# and harmonic transformaion
......
......@@ -17,12 +17,18 @@
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
def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
......@@ -76,7 +82,7 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
return ht(vol*A*ducktape(h_space, None, name))
def MfCorrelatedField(target, amplitudes, name='xi'):
def MfCorrelatedField(target, amplitudes, q=1e-14, 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.
......@@ -131,4 +137,15 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
A = pd @ a
# For `vol` see comment in `CorrelatedField`
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
return ht(vol*A*ducktape(hsp, None, name))
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))
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