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

...
 
Commits (4)
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
import nifty5 as ift
print("N1,N2,std", flush=True)
Ns = 2 ** np.arange(4, 11, dtype=np.int32)
for N1, N2 in [(i, j) for i in Ns for j in Ns]:
np.random.seed(42)
# Two-dimensional signal on RGs with inhomogeneous exposure
position_space = ift.RGSpace(N1)
energy_space = ift.RGSpace(N2)
sky_target = ift.DomainTuple.make((position_space, energy_space))
# Define harmonic space and harmonic transform
harmonic_space_p = position_space.get_default_codomain()
harmonic_space_e = energy_space.get_default_codomain()
# Shared configuration parameters of Amplitude Models
conf_dict = {
'n_pix': 64, # 64 spectral bins
# Spectral smoothness (affects Gaussian process part)
'a': 3, # relatively high variance of spectral curvature
'k0': 0.4, # quefrency mode below which cepstrum flattens
# Power-law part of spectrum
'sv': 1., # Prior standard deviation of -||-
'im': 0., # Expected y-intercept of power law
'iv': 0. # Prior standard deviation of -||-
}
amp_p = ift.SLAmplitude(**conf_dict,
target=ift.PowerSpace(harmonic_space_p),
keys=['tau_p', 'phi_p'],
sm=-4.)
amp_e = ift.SLAmplitude(**conf_dict,
target=ift.PowerSpace(harmonic_space_e),
keys=['tau_e', 'phi_e'],
sm=-3.)
# Define sky operator
# za=4, zq=5 leads to an amplitude median of one
log_diffuse = ift.MultiDimCorrelatedField(sky_target, [amp_p, amp_e],
za=5., zq=1.)
N_stds = 10**3 * max(10, int(Ns.max() / np.sqrt(N1 * N2)))
stds = np.empty(N_stds)
for i in range(N_stds):
mock_pos = ift.from_random('normal', log_diffuse.domain)
stds[i] = log_diffuse(mock_pos).to_global_data().std()
print(f"{N1},{N2},{stds.mean()}", flush=True)
......@@ -24,9 +24,11 @@ from ..domain_tuple import DomainTuple
from ..operators.distributors import PowerDistributor
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 ..library.inverse_gamma_operator import InverseGammaOperator
def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
"""Constructs an operator which turns a white Gaussian excitation field
into a correlated field.
......@@ -130,27 +132,48 @@ def MultiDimCorrelatedField(target, amplitudes, za, zq, name='xi'):
n_subdoms = len(tgt)
hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
print("hsp:", hsp)
psp = [aa.target[0] for aa in amplitudes]
# Assemble PowerDistributor stack
pd = PowerDistributor(hsp, psp[0], 0)
for i in range(1, n_subdoms - 1):
pd = pd @ PowerDistributor(pd.domain, psp[i], i)
pd = pd @ PowerDistributor(pd.domain, psp[-1], n_subdoms - 1)
# Calculate individual power spectra and spread them
# into their respective constant dimensions
p_i = []
print("pd.domain:", pd.domain)
# 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
pd = PowerDistributor(hsp[i], amplitudes[i].target[0])
spreader = ContractionOperator(hsp, i).adjoint
p_i.append(spreader @ pd @ amplitudes[i])
# Calculate global power spectrum as the
# outer product of individual power spectra
p = reduce(mul, p_i)
# Assemble spreader stack for this amplitude
spr = ScalingOperator(1., pd.domain)
for j in range(n_subdoms):
if j < i:
spr = spr @ ContractionOperator(spr.domain, 0).adjoint
if j > i:
spr = spr @ ContractionOperator(spr.domain, 1).adjoint
print("spreader:", spr)
print("amp.target:", amplitudes[i].target)
print("spreader.domain:", spr.domain)
a_i.append(spr @ amplitudes[i])
# Calculate global spectrum amplitude as the
# outer product of individual spectra amplitudes
a = reduce(mul, a_i)
# Power spectrum
p = pd @ a
# Fix pixel variance to one by dividing the global power spectrum
# by the integral of itself.
# Integration volume factor included later on
vdot = VdotOperator(Field.full(p.target, 1.))
normed_pspec = p * ((vdot.adjoint@vdot)(p))**(-1)
normed_pspec = p * ((vdot.adjoint @ vdot)(p))**(-1)
# Generate global power spectrum amplitude operator
vdot_ig = VdotOperator(Field.full(p.target, 1.))
......@@ -171,11 +194,10 @@ def MultiDimCorrelatedField(target, amplitudes, za, zq, name='xi'):
htn = HarmonicTransformOperator(ht.target, target=tgt[i], space=i)
ht = htn @ ht
htn = HarmonicTransformOperator(ht.target,
target=tgt[n_subdoms - 1],
target=tgt[-1],
space=n_subdoms - 1)
ht = htn @ ht
# FIXME remove all zero-mode DoF in SLA _OR_ reintroduce normed_pspec
return ht(vol * (normed_pspec * ducktape(hsp, None, name)))
# return ht(vol * (normed_pspec * ducktape(hsp, None, name)) *
# global_amplitude.ducktape(name + '_global_amplitude'))
......@@ -74,9 +74,9 @@ def CepstrumOperator(target, a, k0):
Larger values result in a weaker constraining prior.
"""
a = float(a)
target = DomainTuple.make(target)
if a <= 0:
raise ValueError
target = DomainTuple.make(target)
if len(target) > 1 or target[0].harmonic:
raise TypeError
if isinstance(k0, (float, int)):
......@@ -89,12 +89,12 @@ def CepstrumOperator(target, a, k0):
raise ValueError
qht = QHTOperator(target)
dom = qht.domain[0]
sym = SymmetrizingOperator(target)
# Compute cepstrum field
dim = len(dom.shape)
dom = qht.domain[0]
shape = dom.shape
dim = len(shape)
q_array = dom.get_k_array()
# Fill all non-zero modes
no_zero_modes = (slice(1, None),)*dim
......