Commit 364eab62 authored by Philipp Arras's avatar Philipp Arras

Formatting

parent 6c984d32
......@@ -19,10 +19,12 @@ import nifty5 as ift
from functools import reduce
from operator import mul
def CorrelatedFieldNormAmplitude(target, amplitudes,
def CorrelatedFieldNormAmplitude(target,
amplitudes,
stdmean,
stdstd,
names=['xi','std']):
names=['xi', 'std']):
"""Constructs an operator which turns white Gaussian excitation fields
into a correlated field defined on a DomainTuple with n entries and n
separate correlation structures.
......@@ -56,43 +58,48 @@ def CorrelatedFieldNormAmplitude(target, amplitudes,
one needs to combine this operator with a :class:`FieldZeroPadder` or even
two (one for the energy and one for the spatial subdomain)
"""
amps = [amplitudes,] if isinstance(amplitudes,ift.Operator) else amplitudes
amps = [
amplitudes,
] if isinstance(amplitudes, ift.Operator) else amplitudes
tgt = ift.DomainTuple.make(target)
if len(tgt) != len(amps):
raise ValueError
stdmean, stdstd = float(stdmean), float(stdstd)
if stdstd <= 0:
raise ValueError
psp = [aa.target[0] for aa in amps]
hsp = ift.DomainTuple.make([tt.get_default_codomain() for tt in tgt])
ht = ift.HarmonicTransformOperator(hsp, target=tgt[0], space=0)
pd = ift.PowerDistributor(hsp, psp[0], 0)
for i in range(1,len(amps)):
ht = ift.HarmonicTransformOperator(ht.target,target=tgt[i], space=i) @ ht
pd = pd @ ift.PowerDistributor(pd.domain, psp[i], space = i)
for i in range(1, len(amps)):
ht = ift.HarmonicTransformOperator(ht.target, target=tgt[i],
space=i) @ ht
pd = pd @ ift.PowerDistributor(pd.domain, psp[i], space=i)
spaces = tuple(range(len(amps)))
a = ift.ContractionOperator(pd.domain, spaces[1:]).adjoint @ amps[0]
for i in range(1,len(amps)):
a = a * (ift.ContractionOperator(pd.domain,
spaces[:i] + spaces[(i+1):]).adjoint @ amps[i])
expander = ift.VdotOperator(ift.full(a.target,1.)).adjoint
Std = stdstd * ift.ducktape(expander.domain, None, names[1])
Std = expander @ (ift.Adder(ift.full(expander.domain,stdmean))@Std).exp()
A = pd @ (Std * a)
for i in range(1, len(amps)):
a = a*(ift.ContractionOperator(pd.domain, spaces[:i] + spaces[
(i + 1):]).adjoint @ amps[i])
expander = ift.VdotOperator(ift.full(a.target, 1.)).adjoint
Std = stdstd*ift.ducktape(expander.domain, None, names[1])
Std = expander @ (
ift.Adder(ift.full(expander.domain, stdmean)) @ Std).exp()
A = pd @ (Std*a)
# For `vol` see comment in `CorrelatedField`
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
return ht(vol*A*ift.ducktape(hsp, None, names[0]))
if __name__ == '__main__':
import numpy as np
from normalized_amplitude import NormalizedAmplitude
......@@ -110,7 +117,6 @@ if __name__ == '__main__':
target = ift.PowerSpace(hspace)
a.append(NormalizedAmplitude(target, 16, 1, 1, -3, 1, 0, 1, 0, 1))
tgt = ift.makeDomain(tuple(spaces))
op = CorrelatedFieldNormAmplitude(tgt,a,0.,1.)
op = CorrelatedFieldNormAmplitude(tgt, a, 0., 1.)
fld = ift.from_random('normal', op.domain)
ift.extra.check_jacobian_consistency(op, fld)
\ No newline at end of file
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