Commit 54965760 authored by Martin Reinecke's avatar Martin Reinecke

merge NIFTy_4

parent dc61ecfd
......@@ -64,14 +64,16 @@ Installation
### Requirements
- [Python](https://www.python.org/) (v2.7.x or 3.5.x)
- [NumPy](https://www.numpy.org/)
- [pyFFTW](https://pypi.python.org/pypi/pyFFTW)
- [Python](https://www.python.org/) (v2.7.x or 3.5.x)
- [NumPy](https://www.numpy.org/)
- [pyFFTW](https://pypi.python.org/pypi/pyFFTW)
Optional dependencies:
- [pyHealpix](https://gitlab.mpcdf.mpg.de/ift/pyHealpix)
- [mpi4py](https://mpi4py.scipy.org)
- [matplotlib](https://matplotlib.org/)
- [pyHealpix](https://gitlab.mpcdf.mpg.de/ift/pyHealpix) (for harmonic
transforms involving domains on the sphere)
- [mpi4py](https://mpi4py.scipy.org) (for MPI-parallel execution)
- [matplotlib](https://matplotlib.org/) (for field plotting)
- [SciPy](https://www.scipy.org/) (for additional minimization algorithms)
### Sources
......@@ -79,24 +81,52 @@ The current version of Nifty4 can be obtained by cloning the repository and
switching to the NIFTy_4 branch:
git clone https://gitlab.mpcdf.mpg.de/ift/NIFTy.git
git checkout NIFTy_4
### Installation via pip
### Installation
It is possible to simply install NIFTy with its mandatory dependencies via the command
In the following, we assume a Debian-based distribution. For other
distributions, the "apt" lines will need slight changes.
NIFTy4 and its mandatory dependencies can be installed via:
sudo apt-get install git libfftw3-dev python python-pip python-dev
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_4
The optional dependencies can be installed via
(Note: If you encounter problems related to `pyFFTW`, make sure that you are
using a pip-installed `pyFFTW` package. Unfortunately, some distributions are
shipping an incorrectly configured `pyFFTW` package, which does not cooperate
with the installed `FFTW3` libraries.)
Plotting support is added via:
pip install --user matplotlib
Support for spherical harmonic transforms is added via:
pip install --user mpi4py
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
MPI support is added via:
sudo apt-get install openmpi-bin libopenmpi-dev
pip install --user mpi4py
Scipy-based minimizers are enabled via:
pip install --user scipy
### Installation for Python 3
If you want to run NIFTy with Python 3, you need to make the following changes
to the instructions above:
- in all `apt-get` commands, replace `python-*` by `python3-*`
- in all `pip` commands, replace `pip` by `pip3`
### Running the tests
In oder to run the tests one needs two additional packages:
pip install --user nose parameterized
pip install --user nose parameterized coverage
Afterwards the tests (including a coverage report) can be run using the
following command in the repository root:
......@@ -119,7 +149,7 @@ Please acknowledge the use of NIFTy in your publication(s) by using a
phrase such as the following:
> *"Some of the results in this publication have been derived using the
> NIFTy package [Selig et al., 2013]."*
> NIFTy package [Steininger et al., 2017]."[1]
### References
......@@ -132,7 +162,5 @@ The NIFTy package is licensed under the terms of the
* * * * *
[1] Selig et al., "NIFTy - Numerical Information Field Theory - a
versatile Python library for signal inference", [A&A, vol. 554, id.
A26](https://dx.doi.org/10.1051/0004-6361/201321236), 2013;
[arXiv:1301.4499](https://www.arxiv.org/abs/1301.4499)
[1] Steininger et al., "NIFTy 3 - Numerical Information Field Theory - A Python framework for multicomponent signal inference on HPC clusters", 2017, submitted to PLOS One;
[arXiv:1708.01073](https://arxiv.org/abs/1708.01073)
apt-get install -y build-essential git autoconf libtool pkg-config libfftw3-dev openmpi-bin libopenmpi-dev \
apt-get install -y git libfftw3-dev openmpi-bin libopenmpi-dev \
python python-pip python-dev python-nose python-numpy python-matplotlib python-future python-mpi4py python-scipy \
python3 python3-pip python3-dev python3-nose python3-numpy python3-matplotlib python3-future python3-mpi4py python3-scipy
......@@ -80,8 +80,8 @@
"source": [
"## Wiener Filter: Example\n",
"\n",
"- One-dimensional signal with power spectrum: $$P(k) = P_0\\,\\left(1+\\left(\\frac{k}{k_0}\\right)^2\\right)^{-\\gamma /2},$$\n",
"with $P_0 = 0.2, k_0 = 5, \\gamma = 4$. Recall: $P(k)$ defines an isotropic and homogeneous $S$.\n",
"- We assume statistical homogeneity and isotropy. Therefore the signal covariance $S$ is diagonal in harmonic space, and is described by a one-dimensional power spectrum, assumed here as $$P(k) = P_0\\,\\left(1+\\left(\\frac{k}{k_0}\\right)^2\\right)^{-\\gamma /2},$$\n",
"with $P_0 = 0.2, k_0 = 5, \\gamma = 4$.\n",
"- $N = 0.2 \\cdot \\mathbb{1}$.\n",
"- Number of data points $N_{pix} = 512$.\n",
"- reconstruction in harmonic space.\n",
......@@ -172,7 +172,7 @@
" inverter = ift.ConjugateGradient(controller=IC)\n",
" # WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy\n",
" # helper methods.\n",
" return ift.library.WienerFilterCurvature(R,N,Sh,inverter)\n"
" return ift.library.WienerFilterCurvature(R,N,Sh,inverter)"
]
},
{
......@@ -296,14 +296,14 @@
"source": [
"s_power = ift.power_analyze(sh)\n",
"m_power = ift.power_analyze(m)\n",
"s_power_data = s_power.val.real\n",
"m_power_data = m_power.val.real\n",
"s_power_data = s_power.val\n",
"m_power_data = m_power.val\n",
"\n",
"# Get signal data and reconstruction data\n",
"s_data = HT(sh).val.real\n",
"m_data = HT(m).val.real\n",
"s_data = HT(sh).val\n",
"m_data = HT(m).val\n",
"\n",
"d_data = d.val.real"
"d_data = d.val"
]
},
{
......@@ -328,9 +328,9 @@
"outputs": [],
"source": [
"plt.figure(figsize=(15,10))\n",
"plt.plot(s_data, 'g', label=\"Signal\")\n",
"plt.plot(d_data, 'k+', label=\"Data\")\n",
"plt.plot(m_data, 'r', label=\"Reconstruction\")\n",
"plt.plot(s_data, 'r', label=\"Signal\", linewidth=3)\n",
"plt.plot(d_data, 'k.', label=\"Data\")\n",
"plt.plot(m_data, 'k', label=\"Reconstruction\",linewidth=3)\n",
"plt.title(\"Reconstruction\")\n",
"plt.legend()\n",
"plt.show()"
......@@ -347,9 +347,9 @@
"outputs": [],
"source": [
"plt.figure(figsize=(15,10))\n",
"plt.plot(s_data - s_data, 'g', label=\"Signal\")\n",
"plt.plot(d_data - s_data, 'k+', label=\"Data\")\n",
"plt.plot(m_data - s_data, 'r', label=\"Reconstruction\")\n",
"plt.plot(s_data - s_data, 'r', label=\"Signal\", linewidth=3)\n",
"plt.plot(d_data - s_data, 'k.', label=\"Data\")\n",
"plt.plot(m_data - s_data, 'k', label=\"Reconstruction\",linewidth=3)\n",
"plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5)\n",
"plt.title(\"Residuals\")\n",
"plt.legend()\n",
......@@ -383,9 +383,9 @@
"ymin = min(m_power_data)\n",
"plt.ylim(ymin, 1)\n",
"xs = np.arange(1,int(N_pixels/2),.1)\n",
"plt.plot(xs, pow_spec(xs), label=\"True Power Spectrum\", linewidth=.7, color='k')\n",
"plt.plot(s_power_data, 'g', label=\"Signal\")\n",
"plt.plot(m_power_data, 'r', label=\"Reconstruction\")\n",
"plt.plot(xs, pow_spec(xs), label=\"True Power Spectrum\", color='k',alpha=0.5)\n",
"plt.plot(s_power_data, 'r', label=\"Signal\")\n",
"plt.plot(m_power_data, 'k', label=\"Reconstruction\")\n",
"plt.axhline(noise_amplitude**2 / N_pixels, color=\"k\", linestyle='--', label=\"Noise level\", alpha=.5)\n",
"plt.axhspan(noise_amplitude**2 / N_pixels, ymin, facecolor='0.9', alpha=.5)\n",
"plt.title(\"Power Spectrum\")\n",
......@@ -494,12 +494,7 @@
},
"outputs": [],
"source": [
"sc = ift.probing.utils.StatCalculator()\n",
"for i in range(200):\n",
" print i\n",
" sc.add(HT(curv.generate_posterior_sample()))\n",
"\n",
"m_var = sc.var"
"m_mean, m_var = ift.probe_with_posterior_samples(curv, m, HT, 200)"
]
},
{
......@@ -523,40 +518,17 @@
},
"outputs": [],
"source": [
"s_power = ift.power_analyze(sh)\n",
"m_power = ift.power_analyze(m)\n",
"s_power_data = s_power.val.real\n",
"m_power_data = m_power.val.real\n",
"\n",
"# Get signal data and reconstruction data\n",
"s_data = s.val.real\n",
"m_data = HT(m).val.real\n",
"m_var_data = m_var.val.real\n",
"uncertainty = np.sqrt(np.abs(m_var_data))\n",
"d_data = d.val.real\n",
"s_data = s.val\n",
"m_data = HT(m).val\n",
"m_var_data = m_var.val\n",
"uncertainty = np.sqrt(m_var_data)\n",
"d_data = d.val\n",
"\n",
"# Set lost data to NaN for proper plotting\n",
"d_data[d_data == 0] = np.nan"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"plt.figure(figsize=(15,10))\n",
"plt.plot(s_data, 'g', label=\"Signal\", linewidth=1)\n",
"plt.plot(d_data, 'k+', label=\"Data\", alpha=1)\n",
"plt.axvspan(l, h, facecolor='0.8', alpha=.5)\n",
"plt.title(\"Incomplete Data\")\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": null,
......@@ -568,11 +540,11 @@
"outputs": [],
"source": [
"fig = plt.figure(figsize=(15,10))\n",
"plt.plot(s_data, 'g', label=\"Signal\", alpha=1, linewidth=1)\n",
"plt.plot(d_data, 'k+', label=\"Data\", alpha=.5)\n",
"plt.plot(m_data, 'r', label=\"Reconstruction\")\n",
"plt.axvspan(l, h, facecolor='0.8', alpha=.5)\n",
"plt.fill_between(range(N_pixels), m_data - uncertainty, m_data + uncertainty, facecolor='0')\n",
"plt.axvspan(l, h, facecolor='0.8',alpha=0.5)\n",
"plt.fill_between(range(N_pixels), m_data - uncertainty, m_data + uncertainty, facecolor='0.5', alpha=0.5)\n",
"plt.plot(s_data, 'r', label=\"Signal\", alpha=1, linewidth=3)\n",
"plt.plot(d_data, 'k.', label=\"Data\")\n",
"plt.plot(m_data, 'k', label=\"Reconstruction\", linewidth=3)\n",
"plt.title(\"Reconstruction of incomplete data\")\n",
"plt.legend()"
]
......@@ -595,13 +567,11 @@
"outputs": [],
"source": [
"N_pixels = 256 # Number of pixels\n",
"sigma2 = 2. # Noise variance\n",
"\n",
"sigma2 = 2. # Noise variance\n",
"\n",
"def pow_spec(k):\n",
" P0, k0, gamma = [.2, 2, 4]\n",
" return P0 * (1. + (k/k0)**2)**(- gamma / 2)\n",
"\n",
" return P0 * (1. + (k/k0)**2)**(-gamma/2)\n",
"\n",
"s_space = ift.RGSpace([N_pixels, N_pixels])"
]
......@@ -649,28 +619,17 @@
"m = D(j)\n",
"\n",
"# Uncertainty\n",
"sc = ift.probing.utils.StatCalculator()\n",
"\n",
"IC = ift.GradientNormController(iteration_limit=50000,\n",
" tol_abs_gradnorm=0.1)\n",
"inverter = ift.ConjugateGradient(controller=IC)\n",
"curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter)\n",
"\n",
"for i in range(20):\n",
" print i\n",
" sc.add(HT(curv.generate_posterior_sample()))\n",
"\n",
"m_var = sc.var\n",
"m_mean, m_var = ift.probe_with_posterior_samples(curv, m, HT, 20)\n",
"\n",
"# Get data\n",
"s_power = ift.power_analyze(sh)\n",
"m_power = ift.power_analyze(m)\n",
"s_power_data = s_power.val.real\n",
"m_power_data = m_power.val.real\n",
"s_data = HT(sh).val.real\n",
"m_data = HT(m).val.real\n",
"m_var_data = m_var.val.real\n",
"d_data = d.val.real\n",
"s_power_data = s_power.val\n",
"m_power_data = m_power.val\n",
"s_data = HT(sh).val\n",
"m_data = HT(m).val\n",
"m_var_data = m_var.val\n",
"d_data = d.val\n",
"\n",
"uncertainty = np.sqrt(np.abs(m_var_data))"
]
......@@ -718,10 +677,12 @@
"mi = np.min(s_data)\n",
"ma = np.max(s_data)\n",
"\n",
"fig, axes = plt.subplots(2, 2, figsize=(15, 15))\n",
"fig, axes = plt.subplots(3, 2, figsize=(15, 22.5))\n",
"samp1 = HT(curv.draw_sample()+m).val\n",
"samp2 = HT(curv.draw_sample()+m).val\n",
"\n",
"data = [s_data, m_data, s_data - m_data, uncertainty]\n",
"caption = [\"Signal\", \"Reconstruction\", \"Residuals\", \"Uncertainty Map\"]\n",
"data = [s_data, m_data, samp1, samp2, s_data - m_data, uncertainty]\n",
"caption = [\"Signal\", \"Reconstruction\", \"Sample 1\", \"Sample 2\", \"Residuals\", \"Uncertainty Map\"]\n",
"\n",
"for ax in axes.flat:\n",
" im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi, vmax=ma)\n",
......
......@@ -8,11 +8,11 @@ np.random.seed(42)
if __name__ == "__main__":
# Set up position space
s_space = ift.RGSpace([128, 128])
# s_space = ift.HPSpace(32)
#s_space = ift.HPSpace(32)
# Define harmonic transformation and associated harmonic space
h_space = s_space.get_default_codomain()
fft = ift.HarmonicTransformOperator(h_space, s_space)
HT = ift.HarmonicTransformOperator(h_space, s_space)
# Set up power space
p_space = ift.PowerSpace(h_space,
......@@ -34,34 +34,40 @@ if __name__ == "__main__":
# Instrument._diagonal.val[64:512-64, 64:512-64] = 0
# Add a harmonic transformation to the instrument
R = Instrument*fft
R = Instrument*HT
noise = 1.
N = ift.DiagonalOperator(ift.Field.full(s_space, noise).weight(1))
n = ift.Field.from_random(domain=s_space, random_type='normal',
std=np.sqrt(noise), mean=0)
noiseless_data = R(sh)
signal_to_noise = 1.
noise_amplitude = noiseless_data.val.std()/signal_to_noise
N = ift.ScalingOperator(noise_amplitude**2, s_space)
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=noise_amplitude,
mean=0)
# Create mock data
d = R(sh) + n
d = noiseless_data + n
# The information source
j = R.adjoint_times(N.inverse_times(d))
realized_power = ift.log(ift.power_analyze(sh,
binbounds=p_space.binbounds))
data_power = ift.log(ift.power_analyze(fft.adjoint_times(d),
data_power = ift.log(ift.power_analyze(HT.adjoint_times(d),
binbounds=p_space.binbounds))
d_data = d.val
ift.plot(d, name="data.png")
plotdict = {"colormap": "Planck-like"}
ift.plot(d, name="data.png", **plotdict)
ift.plot(HT(sh), name="mock_signal.png", **plotdict)
IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
tol_abs_gradnorm=0.1)
tol_abs_gradnorm=1e-15)
minimizer = ift.RelaxedNewton(IC1)
ICI = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=0.1)
ICI = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-10)
map_inverter = ift.ConjugateGradient(controller=ICI)
ICI2 = ift.GradientNormController(iteration_limit=200,
tol_abs_gradnorm=1e-5)
tol_abs_gradnorm=1e-10)
power_inverter = ift.ConjugateGradient(controller=ICI2)
# Set starting position
......@@ -91,4 +97,4 @@ if __name__ == "__main__":
# Plot current estimate
ift.dobj.mprint(i)
if i % 50 == 0:
ift.plot(fft(m0), name='map.png')
ift.plot(HT(m0), name='map.png', **plotdict)
import nifty4 as ift
import numpy as np
if __name__ == "__main__":
np.random.seed(42)
# Setting up parameters
correlation_length = 1. # Typical distance over which the field is correlated
field_variance = 2. # Variance of field in position space
response_sigma = 0.02 # Smoothing length of response (in same unit as L)
signal_to_noise = 100 # The signal to noise ratio
np.random.seed(43) # Fixing the random seed
def power_spectrum(k): # Defining the power spectrum
a = 4 * correlation_length * field_variance**2
return a / (1 + k * correlation_length) ** 4
# Setting up the geometry
L = 2. # Total side-length of the domain
N_pixels = 128 # Grid resolution (pixels per axis)
# signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
signal_space = ift.HPSpace(16)
harmonic_space = signal_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space)
power_space = ift.PowerSpace(harmonic_space)
# Creating the mock signal
S = ift.create_power_operator(harmonic_space,
power_spectrum=power_spectrum)
mock_power = ift.PS_field(power_space, power_spectrum)
mock_signal = ift.power_synthesize(mock_power, real_signal=True)
# Setting up an exemplary response
mask = ift.Field.ones(signal_space)
N10 = int(N_pixels/10)
# mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ift.GeometryRemover(signal_space)
R = R*ift.DiagonalOperator(mask)
R = R*ht
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
data_domain = R.target[0]
# Setting up the noise covariance and drawing a random noise realization
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.val.std()/signal_to_noise
N = ift.DiagonalOperator(
ift.Field.full(data_domain, noise_amplitude**2))
noise = ift.Field.from_random(
domain=data_domain, random_type='normal',
std=noise_amplitude, mean=0)
data = noiseless_data + noise
# Wiener filter
m0 = ift.Field.zeros(harmonic_space)
ctrl = ift.GradientNormController(tol_abs_gradnorm=0.0001)
ctrl2 = ift.GradientNormController(tol_abs_gradnorm=0.1, name="outer")
inverter = ift.ConjugateGradient(controller=ctrl)
energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R,
N, S, inverter=inverter)
#minimizer = ift.VL_BFGS(controller=ctrl2, max_history_length=20)
minimizer = ift.RelaxedNewton(controller=ctrl2)
#minimizer = ift.SteepestDescent(controller=ctrl2)
me = minimizer(energy)
m = ht(me[0].position)
# Plotting
plotdict = {"colormap": "Planck-like"}
ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict)
logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape)
ift.plot(ift.Field(signal_space, val=ift.dobj.from_global_data(logdata)),
name="log_of_data.png", **plotdict)
ift.plot(m, name='m.png', **plotdict)
# Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober):
pass
proby = Proby(signal_space, probe_count=1)
proby(lambda z: ht(me2[0].curvature.inverse_times(ht.adjoint_times(z))))
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.02)
variance = sm(proby.diagonal.weight(-1))
ift.plot(variance, name='variance.png', **plotdict)
......@@ -66,8 +66,7 @@ if __name__ == "__main__":
true_sky = nonlinearity(HT(power*sh))
noiseless_data = MeasurementOperator(true_sky)
noise_amplitude = noiseless_data.val.std()*noise_level
N = ift.DiagonalOperator(
ift.Field.full(d_space, noise_amplitude**2))
N = ift.ScalingOperator(noise_amplitude**2, d_space)
n = ift.Field.from_random(
domain=d_space, random_type='normal',
std=noise_amplitude, mean=0)
......
......@@ -83,7 +83,6 @@ if __name__ == "__main__":
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.val.std()/signal_to_noise
# Setting up the noise covariance and drawing a random noise realization
#ndiag = ift.Field.full(data_domain, noise_amplitude**2)
N = ift.ScalingOperator(noise_amplitude**2, data_domain)
noise = ift.Field.from_random(
domain=data_domain, random_type='normal',
......@@ -100,8 +99,7 @@ if __name__ == "__main__":
m_k = wiener_curvature.inverse_times(j)
m = ht(m_k)
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
plotdict = {"colormap": "Planck-like"}
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
ift.plot(ift.Field(plot_space,val=ht(mock_signal).val), name='mock_signal.png',
**plotdict)
......@@ -110,3 +108,4 @@ if __name__ == "__main__":
# sampling the uncertainty map
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 10)
ift.plot(ift.Field(plot_space, val=ift.sqrt(variance).val), name="uncertainty.png", **plotdict)
ift.plot(ift.Field(plot_space, val=(mean+m).val), name="posterior_mean.png", **plotdict)
......@@ -44,8 +44,7 @@ if __name__ == "__main__":
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.val.std()/signal_to_noise
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, noise_amplitude**2)
N = ift.DiagonalOperator(ndiag)
N = ift.ScalingOperator(noise_amplitude**2, data_domain)
noise = ift.Field.from_random(
domain=data_domain, random_type='normal',
std=noise_amplitude, mean=0)
......@@ -60,20 +59,13 @@ if __name__ == "__main__":
m_k = wiener_curvature.inverse_times(j)
m = ht(m_k)
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
plotdict = {"colormap": "Planck-like"}
ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict)
ift.plot(ift.Field(signal_space, val=data.val),
name="data.png", **plotdict)
ift.plot(m, name="map.png", **plotdict)
# sampling the uncertainty map
sample_variance = ift.Field.zeros(signal_space)
sample_mean = ift.Field.zeros(signal_space)
n_samples = 10
for i in range(n_samples):
sample = ht(wiener_curvature.generate_posterior_sample()) + m
sample_variance += sample**2
sample_mean += sample
variance = sample_variance/n_samples - (sample_mean/n_samples)**2
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 5)
ift.plot(ift.sqrt(variance), name="uncertainty.png", **plotdict)
ift.plot(mean+m, name="posterior_mean.png", **plotdict)
import nifty4 as ift
import numpy as np
np.random.seed(42)
class DiagonalProber(ift.DiagonalProberMixin, ift.Prober):
pass
class MultiProber(ift.DiagonalProberMixin, ift.TraceProberMixin, ift.Prober):
pass
np.random.seed(42)
x = ift.RGSpace((8, 8))
f = ift.Field.from_random(domain=x, random_type='normal')
diagOp = ift.DiagonalOperator(f)
diagProber = DiagonalProber(domain=x)
diagProber(diagOp)
ift.dobj.mprint((f - diagProber.diagonal).norm())
multiProber = MultiProber(domain=x)
multiProber(diagOp)
ift.dobj.mprint((f - multiProber.diagonal).norm())
ift.dobj.mprint(f.sum() - multiProber.trace)
diag = ift.probe_diagonal(diagOp, 1000)
ift.dobj.mprint((f - diag).norm())
......@@ -28,7 +28,7 @@ if __name__ == "__main__":
# Set up the geometry
s_space = ift.RGSpace([N_pixels, N_pixels], distances=pixel_width)
h_space = s_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(h_space, s_space)
HT = ift.HarmonicTransformOperator(h_space, s_space)
p_space = ift.PowerSpace(h_space)
# Create mock data
......@@ -38,12 +38,13 @@ if __name__ == "__main__":
sp = ift.PS_field(p_space, pow_spec)
sh = ift.power_synthesize(sp, real_signal=True)
R = ht*ift.create_harmonic_smoothing_operator((h_space,),0,response_sigma)
R = HT*ift.create_harmonic_smoothing_operator((h_space,), 0,
response_sigma)
noiseless_data = R(sh)
signal_to_noise = 1.
noise_amplitude = noiseless_data.val.std()/signal_to_noise
N = ift.DiagonalOperator(ift.Field.full(s_space, noise_amplitude**2))
N = ift.ScalingOperator(noise_amplitude**2, s_space)
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=noise_amplitude,
......@@ -62,9 +63,11 @@ if __name__ == "__main__":
D = ift.InversionEnabler(D, inverter)
m = D(j)
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
ift.plot(ht(sh), name="mock_signal.png", **plotdict)
ift.plot(ift.Field(s_space, val=d.val),
name="data.png", **plotdict)
ift.plot(ht(m), name="map.png", **plotdict)
# Plotting
d_field = ift.Field(s_space, val=d.val)
zmax = max(HT(sh).max(), d_field.max(), HT(m).max())
zmin = min(HT(sh).min(), d_field.min(), HT(m).min())
plotdict = {"colormap": "Planck-like", "zmax": zmax, "zmin": zmin}
ift.plot(HT(sh), name="mock_signal.png", **plotdict)
ift.plot(d_field, name="data.png", **plotdict)
ift.plot(HT(m), name="reconstruction.png", **plotdict)
......@@ -3,10 +3,24 @@ import nifty4 as ift
import numericalunits as nu
def init_nu():
if ift.dobj.ntask > 1:
from mpi4py import MPI
comm = MPI.COMM_WORLD
data = np.random.randint(10000000, size=1)
data = comm.bcast(data, root=0)
nu.reset_units(data[0])
if __name__ == "__main__":
# In MPI mode, the random seed for numericalunits must be set by hand