...
 
Commits (6)
# 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.MultiCorrelatedField(sky_target, [amp_p, amp_e],
va=5.,
vq=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)
......@@ -91,11 +91,13 @@ if __name__ == '__main__':
'sm': -4., # Expected exponent of power law
'sv': 1., # Prior standard deviation of -||-
'im': 0., # Expected y-intercept of power law
'iv': 0., # Prior standard deviation of -||-
'iv': 2., # Prior standard deviation of -||-
# zero-mode prior
'zm': 2.5e-4, # Prior mean
'zv': 1.0e-4 # Prior variance
#'zm': 2.5e-4, # Prior mean
#'zv': 1.0e-4 # Prior variance
'zm': None, # Prior mean
'zv': None # Prior variance
}
amp_e_dict = {
......@@ -111,11 +113,13 @@ if __name__ == '__main__':
'sm': -3., # Expected exponent of power law
'sv': 1., # Prior standard deviation of -||-
'im': 0., # Expected y-intercept of power law
'iv': 0., # Prior standard deviation of -||-
'iv': 1., # Prior standard deviation of -||-
# zero-mode prior
'zm': 3.0e-3, # Prior mean
'zv': 2.0e-3 # Prior variance
#'zm': 3.0e-3, # Prior mean
#'zv': 2.0e-3 # Prior variance
'zm': None, # Prior mean
'zv': None # Prior variance
}
amp_p = ift.SLAmplitude(**amp_p_dict)
......@@ -198,15 +202,20 @@ if __name__ == '__main__':
amplitude_specs = [(p_e, 'energy')]
for p, n in amplitude_specs:
sample_specs = [
ift.power_analyze(p.force(pos + s)) for s in KL.samples
]
sc = ift.StatCalculator()
for sample in KL.samples:
sc.add(ift.power_analyze(p.force(pos + sample)))
for spec in sample_specs:
sc.add(spec)
plot = ift.Plot()
plot.add(
sample_specs +
[sc.mean, ift.power_analyze(p.force(mock_position))],
linewidth=[1., 3.],
label=['Posterior Mean', 'Ground Truth'],
title=f"Posterior {n} spectrum")
linewidth=[1.] * len(sample_specs) + [3., 3.],
label=[None] * len(sample_specs) +
['Posterior Mean', 'Ground Truth'],
title=f"Sampled Posterior {n} spectrum")
plot.output(xsize=12,
ysize=10,
name=filename.format(f'{iteration:02d}-power_{n}'))
......@@ -214,22 +223,10 @@ if __name__ == '__main__':
print("Saved results as '{}'.".format(filename.format('*')))
# Print InverseGamma prior fits
if mode == 1:
keys = ['zm_e', 'zm_p', 'xi_pspec_scaling']
hyperparams = [(amp_e_dict['zm'], amp_e_dict['zv']),
(amp_p_dict['zm'], amp_p_dict['zv']),
(diff_gs_params['gs_alpha'], diff_gs_params['gs_q'])]
elif mode == 0:
keys = ['zm_e', 'xi_pspec_scaling']
hyperparams = [(amp_e_dict['zm'], amp_e_dict['zv']),
(diff_gs_params['gs_alpha'], diff_gs_params['gs_q'])]
print("InverseGamma fits:")
for key, hpar in zip(keys, hyperparams):
if key in ['zm_e', 'zm_p']:
op = InvGammaFromMeanVar(sky.domain[key], *hpar)
else:
op = ift.InverseGammaOperator(sky.domain[key], *hpar)
print(f"{key}:")
print(f"\tmock: {float(op(mock_position[key]).local_data):+1.2e}")
print(f"\t rec: {float(op(pos[key]).local_data):+1.2e}\n")
print("Variance Prior fit:")
key = 'xi_pspec_scaling'
op = ift.InverseGammaOperator(sky.domain[key], diff_gs_params['gs_alpha'],
diff_gs_params['gs_q'])
print(f"{key}:")
print(f"\tmock: {float(op(mock_position[key]).local_data):+1.2e}")
print(f"\t rec: {float(op(pos[key]).local_data):+1.2e}\n")
......@@ -76,8 +76,8 @@ if __name__ == '__main__':
'iv': .3, # y-intercept variance
# Zero mode prior (InverseGammaOperator)
'zm': 2., # Prior mean
'zv': 0.5 # Prior variance
'zm': 1.0, # Prior mean
'zv': 0.01 # Prior variance
}
A = ift.SLAmplitude(**dct)
......
# coding: utf-8
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
df = pd.read_csv('stds.csv')
N1 = np.unique(df['N1'])
N2 = np.unique(df['N2'])
# assume N1 and N2 go through the same values for ease
assert np.allclose(N1, N2)
N = N1
logN = np.log2(N)
nN = len(N)
# read data
# assume no combinations were skipped
stds = np.empty((nN, nN))
df_stds = df['std'].to_numpy()
for i in range(nN):
stds[i] = df_stds[df['N1'] == N1[i]]
# plot all data
logNmin = logN.min()
logNmax = logN.max()
logNdiff = (logNmax - logNmin) / nN
extents = [logNmin - logNdiff/2, logNmax + logNdiff/2]
extents = extents + extents
im = plt.imshow(stds.transpose(), origin='lower', extent=extents, aspect=1.)
plt.title("Pixel Standard Deviation")
plt.xlabel("log2(N1)")
plt.ylabel("log2(N2)")
#plt.ticklabel_format(style='sci', scilimits=(0, 0))
plt.colorbar(im, orientation='horizontal')
plt.savefig("stds.png")
plt.close()
# plot data slices
for cname, xname in [('N1', 'N2'), ('N2', 'N1')]:
for n in N:
idx = df[cname] == n
plt.plot(df[xname][idx], df['std'][idx], label=f"{cname} = {n}")
plt.xscale('log')
plt.xlabel(xname)
plt.ylabel('std')
plt.legend(loc='upper right')
plt.title(f"Pixel Standard Deviation")
plt.savefig(f"stds-{cname}-const.png")
plt.close()
N1,N2,std
16, 16, 0.15656182044290326
16, 32, 0.14258513509868775
16, 64, 0.13403633010259242
16, 128, 0.12688340189273692
16, 256, 0.12209517082737151
16, 512, 0.11924802673814003
16, 1024, 0.11501015017748917
32, 16, 0.15015620262999732
32, 32, 0.13703995151953716
32, 64, 0.1278865371763398
32, 128, 0.12100677520481282
32, 256, 0.1180975109995972
32, 512, 0.1129459053464104
32, 1024, 0.11123990590779961
64, 16, 0.1467013112071079
64, 32, 0.1331952125728305
64, 64, 0.12451510216327419
64, 128, 0.11884546632988117
64, 256, 0.11359660631274206
64, 512, 0.11102942631069113
64, 1024, 0.1088691294817882
128, 16, 0.14372298653287052
128, 32, 0.1312031976090056
128, 64, 0.12301993124270796
128, 128, 0.11593919637894273
128, 256, 0.11254965750398203
128, 512, 0.10939315156963356
128, 1024, 0.10684179540255825
256, 16, 0.14208744072193152
256, 32, 0.13032890412009882
256, 64, 0.12097424710156424
256, 128, 0.11512535626826093
256, 256, 0.11104087581197146
256, 512, 0.10818579324619022
256, 1024, 0.1061634828421742
512, 16, 0.14140313629166448
512, 32, 0.12807787834570547
512, 64, 0.12034793301610947
512, 128, 0.11463176639095547
512, 256, 0.11076795667077056
512, 512, 0.1074963119411124
512, 1024, 0.10496816518193403
1024, 16, 0.14050226645140723
1024, 32, 0.12774949308444447
1024, 64, 0.11973393622424576
1024, 128, 0.11434272196958291
1024, 256, 0.10986964286804798
1024, 512, 0.10681079847980875
1024, 1024, 0.10440265992667878
......@@ -176,7 +176,7 @@ def MultiCorrelatedField(target, amplitudes, gs_alpha, gs_q, name='xi'):
# -> constant amplitude
# Still got to normalize, tough:
# Since the field contains only ones, 2-norm = sqrt(1-norm)
amplitude_integral = Field.full(spr.domain, 1.).integrate().sqrt()
amplitude_integral = np.sqrt(Field.full(spr.domain, 1.).integrate())
fct /= amplitude_integral
fct = ScalingOperator(fct, pd.domain)
......
......@@ -114,72 +114,6 @@ def CepstrumOperator(target, a, k0):
return sym @ qht @ makeOp(cepstrum.sqrt())
def OneDCepstrumOperator(target, a, k0):
"""Turns a white Gaussian random field into a smooth field on a LogRGSpace.
The smooth field is composed out of three operators:
sym @ qht @ diag(sqrt_ceps),
where sym is a :class:`SymmetrizingOperator`, qht is a :class:`QHTOperator`
and ceps is the so-called cepstrum:
.. math::
\\mathrm{sqrt\\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
These operators are combined in this fashion in order to generate:
- A field which is smooth, i.e. second derivatives are punished (note
that the sqrt-cepstrum is essentially proportional to 1/k**2).
- A field which is symmetric around the pixel in the middle of the space.
This is result of the :class:`SymmetrizingOperator` and needed in order
to decouple the degrees of freedom at the beginning and the end of the
amplitude whenever :class:`CepstrumOperator` is used as in
:class:`SLAmplitude`.
Parameters
----------
target : LogRGSpace
Target domain of the operator, needs to be non-harmonic and 1d.
a : float
Cutoff of smoothness prior (positive only). Controls the
regularization of the inverse laplace operator to be finite at zero.
Larger values for the cutoff results in a weaker constraining prior.
k0 : float
Strength of smoothness prior in quefrency space (positive only) along
each axis. If float then the strength is the same along each axis.
Larger values result in a weaker constraining prior.
"""
a = float(a)
if a <= 0:
raise ValueError
target = DomainTuple.make(target)
if len(target.shape) != 1 or target[0].harmonic:
raise TypeError
if not isinstance(k0, (float, int)):
raise ValueError
if k0 <= 0:
raise ValueError
qht = QHTOperator(target)
sym = SymmetrizingOperator(target)
# Compute cepstrum field
dom = qht.domain[0]
shape = dom.shape
q_array = dom.get_k_array()
# Fill all non-zero modes
no_zero_modes = (slice(1, None),)
ks = q_array[(slice(None),) + no_zero_modes]
cepstrum_field = np.zeros(shape)
cepstrum_field[no_zero_modes] = _ceps_kernel(ks, a, k0)
# zero mode stays zero
cepstrum = Field.from_global_data(dom, cepstrum_field)
return sym @ qht @ makeOp(cepstrum.sqrt())
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, zm, zv,
keys=['tau', 'phi', 'zm']):
'''Operator for parametrizing smooth amplitudes (square roots of power
......@@ -231,10 +165,12 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, zm, zv,
LogRGSpace (see :class:`ExpTransform`).
iv : float
Prior standard deviation of y-intercept of power law.
zm : float
zm : float or None
Prior mean value of the zero mode.
zv : float
If None is given, zero mode is set to zero.
zv : float or None
Prior variance of the zero mode.
If None is given, zero mode is set to zero.
Returns
-------
......@@ -263,9 +199,16 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, zm, zv,
sm, sv, im, iv = float(sm), float(sv), float(im), float(iv)
if sv <= 0 or iv < 0:
raise ValueError
zm, zv = float(zm), float(zv)
if zm <= 0 or zv <= 0:
raise ValueError
if zm is not None and zv is not None:
zero_mode = True
zm, zv = float(zm), float(zv)
if zm <= 0 or zv <= 0:
raise ValueError
if len(keys) < 3:
raise ValueError("Need a key for the zero mode domain")
else:
zero_mode = False
et = ExpTransform(target, n_pix)
dom = et.domain[0]
......@@ -283,18 +226,21 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, zm, zv,
zm_mask = np.ones(dom.shape)
zm_mask[(0,) * len(dom.shape)] = 0.
zm_mask = Field.from_global_data(dom, zm_mask)
zm_mask_inv = Field.full(dom, 1.) - zm_mask
zm_mask = makeOp(zm_mask)
smooth = zm_mask @ CepstrumOperator(dom, **dct).ducktape(keys[0])
# Zero Mode
vdot = VdotOperator(zm_mask_inv)
zb = zm**2 / zv + 1.
ig = InverseGammaOperator(vdot.target, alpha=zb+1., q=zm*zb)
zero_mode = vdot.adjoint @ ig.ducktape(keys[2])
zm_mask_op = makeOp(zm_mask)
smooth = zm_mask_op @ CepstrumOperator(dom, **dct).ducktape(keys[0])
# Combine components
loglog_ampl = 0.5*(smooth + linear) + zero_mode
loglog_ampl = 0.5*(smooth + linear)
# Add optional zero mode prior
if zero_mode:
zm_mask_inv = Field.full(dom, 1.) - zm_mask
vdot = VdotOperator(zm_mask_inv)
b = zm**2 / zv + 1.
log_ig = InverseGammaOperator(vdot.target, alpha=b+1., q=zm*b).log()
zero_mode = vdot.adjoint @ log_ig.ducktape(keys[2])
loglog_ampl = loglog_ampl + zero_mode
# Go from loglog-space to linear-linear-space
return et @ loglog_ampl
# ToDo Development
## In `MultiCorrelatedField`:
- before outer product devide by integral
(Reimar: /= np.linalg.norm(a))
- instead of zero-mode prior multiply the global Amplitude
with sqrt(InvGamma) prior
## LinearSLAmplitude:
- zero mode as log(InverseGammaOperator)
- use uninformative prior
## Demo:
- plot spectrum samples
## Tests:
- Statistics of fields drawn from sub-spectrum are the same as
......