Commit 6b02ddc4 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'NIFTy_4' into citations

parents 90cb9662 446692d3
Pipeline #26240 passed with stage
in 5 minutes and 31 seconds
# custom
# never store the git version file
git_version.py
# custom
setup.cfg
.idea
.DS_Store
......
......@@ -12,12 +12,12 @@ before_script:
- sh ci/install_basics.sh
- pip install --process-dependency-links -r ci/requirements.txt
- pip3 install --process-dependency-links -r ci/requirements.txt
- pip install --user .
- pip3 install --user .
test_min:
stage: test
script:
- pip install --user .
- pip3 install --user .
- nosetests -q
- nosetests3 -q
- OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 4 nosetests -q 2>/dev/null
......@@ -25,3 +25,13 @@ test_min:
- nosetests -q --with-coverage --cover-package=nifty4 --cover-branches --cover-erase
- >
coverage report | grep TOTAL | awk '{ print "TOTAL: "$6; }'
pages:
script:
- sh docs/generate.sh
- mv docs/build/ public/
artifacts:
paths:
- public
only:
- NIFTy_4
#FROM ubuntu:artful
FROM debian:testing-slim
# dependencies via apt
RUN apt-get update
ADD ci/install_basics.sh /tmp/install_basics.sh
RUN sh /tmp/install_basics.sh
# python dependencies
ADD ci/requirements.txt /tmp/requirements.txt
RUN pip install --process-dependency-links -r /tmp/requirements.txt
# copy sources and install nifty
COPY . /tmp/NIFTy
RUN pip install /tmp/NIFTy
# Cleanup
RUN rm -r /tmp/*
......@@ -4,7 +4,7 @@ NIFTy - Numerical Information Field Theory
[![coverage report](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/NIFTy_4/coverage.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/NIFTy_4)
**NIFTy** project homepage:
[https://www.mpa-garching.mpg.de/ift/nifty/](https://www.mpa-garching.mpg.de/ift/nifty/)
[http://ift.pages.mpcdf.de/NIFTy](http://ift.pages.mpcdf.de/NIFTy)
Summary
-------
......@@ -31,47 +31,22 @@ point sets, *n*-dimensional regular grids, spherical spaces, their
harmonic counterparts, and product spaces constructed as combinations of
those.
### Class & Feature Overview
The NIFTy library features three main classes: **Space**s that represent
certain grids, **Field**s that are defined on spaces, and **LinearOperator**s
that apply to fields.
- [Spaces](https://www.mpa-garching.mpg.de/ift/nifty/space.html)
- `RGSpace` - *n*-dimensional regular Euclidean grid
- `LMSpace` - spherical harmonics
- `GLSpace` - Gauss-Legendre grid on the 2-sphere
- `HPSpace` - [HEALPix](https://sourceforge.net/projects/healpix/)
grid on the 2-sphere
- [Fields](https://www.mpa-garching.mpg.de/ift/nifty/field.html)
- `Field` - generic class for (discretized) fields
<!-- -->
Field.conjugate Field.dim Field.norm
Field.vdot Field.weight
- [Operators](https://www.mpa-garching.mpg.de/ift/nifty/operator.html)
- `DiagonalOperator` - purely diagonal matrices in a specified
basis
- `FFTOperator` - conversion between spaces and their harmonic
counterparts
- (and more)
- (and more)
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 +54,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:
......@@ -107,7 +110,7 @@ following command in the repository root:
### First Steps
For a quick start, you can browse through the [informal
introduction](https://www.mpa-garching.mpg.de/ift/nifty/start.html) or
introduction](http://ift.pages.mpcdf.de/NIFTy/code.html) or
dive into NIFTy by running one of the demonstrations, e.g.:
python demos/wiener_filter_via_curvature.py
......@@ -119,7 +122,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]
Or one of the following BibTeX entries:
......@@ -174,7 +177,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
parameterized
coverage
git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
sphinx==1.6.7
sphinx_rtd_theme
numpydoc
......@@ -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)"
]
},
{
......@@ -229,14 +229,13 @@
"s_space = ift.RGSpace(N_pixels)\n",
"h_space = s_space.get_default_codomain()\n",
"HT = ift.HarmonicTransformOperator(h_space, target=s_space)\n",
"p_space = ift.PowerSpace(h_space)\n",
"\n",
"# Operators\n",
"Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)\n",
"R = HT #*ift.create_harmonic_smoothing_operator((h_space,), 0, 0.02)\n",
"\n",
"# Fields and data\n",
"sh = ift.power_synthesize(ift.PS_field(p_space, pow_spec),real_signal=True)\n",
"sh = Sh.draw_sample()\n",
"noiseless_data=R(sh)\n",
"noise_amplitude = np.sqrt(0.2)\n",
"N = ift.ScalingOperator(noise_amplitude**2, s_space)\n",
......@@ -281,7 +280,7 @@
}
},
"source": [
"### Create Power Spectra of Signal and Reconstruction"
"### Signal Reconstruction"
]
},
{
......@@ -294,43 +293,15 @@
},
"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 = HT(sh).val.real\n",
"m_data = HT(m).val.real\n",
"s_data = HT(sh).to_global_data()\n",
"m_data = HT(m).to_global_data()\n",
"d_data = d.to_global_data()\n",
"\n",
"d_data = d.val.real"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Signal Reconstruction"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"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 +318,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",
......@@ -377,15 +348,17 @@
},
"outputs": [],
"source": [
"s_power_data = ift.power_analyze(sh).to_global_data()\n",
"m_power_data = ift.power_analyze(m).to_global_data()\n",
"plt.figure(figsize=(15,10))\n",
"plt.loglog()\n",
"plt.xlim(1, int(N_pixels/2))\n",
"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",
......@@ -420,7 +393,7 @@
"# R is defined below\n",
"\n",
"# Fields\n",
"sh = ift.power_synthesize(ift.PS_field(p_space,pow_spec),real_signal=True)\n",
"sh = Sh.draw_sample()\n",
"s = HT(sh)\n",
"n = ift.Field.from_random(domain=s_space, random_type='normal',\n",
" std=noise_amplitude, mean=0)"
......@@ -450,11 +423,14 @@
"l = int(N_pixels * 0.2)\n",
"h = int(N_pixels * 0.2 * 2)\n",
"\n",
"mask = ift.Field(s_space, val=1)\n",
"mask.val[ l : h] = 0\n",
"mask = np.full(s_space.shape, 1.)\n",
"mask[l:h] = 0\n",
"mask = ift.Field.from_global_data(s_space, mask)\n",
"\n",
"R = ift.DiagonalOperator(mask)*HT\n",
"n.val[l:h] = 0\n",
"n = n.to_global_data()\n",
"n[l:h] = 0\n",
"n = ift.Field.from_global_data(s_space, n)\n",
"\n",
"d = R(sh) + n"
]
......@@ -494,7 +470,7 @@
},
"outputs": [],
"source": [
"m_mean, m_var = ift.probe_with_posterior_samples(curv, m, HT, 200)\n"
"m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 200)"
]
},
{
......@@ -518,40 +494,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.to_global_data()\n",
"m_data = HT(m).to_global_data()\n",
"m_var_data = m_var.to_global_data()\n",
"uncertainty = np.sqrt(m_var_data)\n",
"d_data = d.to_global_data()\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,
......@@ -563,11 +516,11 @@
"outputs": [],
"source": [
"fig = plt.figure(figsize=(15,10))\n",
"plt.plot(s_data, 'g', label=\"Signal\", alpha=1, linewidth=4)\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()"
]
......@@ -590,13 +543,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])"
]
......@@ -613,14 +564,13 @@
"source": [
"h_space = s_space.get_default_codomain()\n",
"HT = ift.HarmonicTransformOperator(h_space,s_space)\n",
"p_space = ift.PowerSpace(h_space)\n",
"\n",
"# Operators\n",
"Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)\n",
"N = ift.ScalingOperator(sigma2,s_space)\n",
"\n",
"# Fields and data\n",
"sh = ift.power_synthesize(ift.PS_field(p_space,pow_spec),real_signal=True)\n",
"sh = Sh.draw_sample()\n",
"n = ift.Field.from_random(domain=s_space, random_type='normal',\n",
" std=np.sqrt(sigma2), mean=0)\n",
"\n",
......@@ -629,11 +579,14 @@
"l = int(N_pixels * 0.33)\n",
"h = int(N_pixels * 0.33 * 2)\n",
"\n",
"mask = ift.Field(s_space, val=1)\n",
"mask.val[l:h,l:h] = 0\n",
"mask = np.full(s_space.shape, 1.)\n",
"mask[l:h,l:h] = 0.\n",
"mask = ift.Field.from_global_data(s_space, mask)\n",
"\n",
"R = ift.DiagonalOperator(mask)*HT\n",
"n.val[l:h, l:h] = 0\n",
"n = n.to_global_data()\n",
"n[l:h, l:h] = 0\n",
"n = ift.Field.from_global_data(s_space, n)\n",
"curv = Curvature(R=R, N=N, Sh=Sh)\n",
"D = curv.inverse\n",
"\n",
......@@ -644,18 +597,13 @@
"m = D(j)\n",
"\n",
"# Uncertainty\n",
"m_mean, m_var = ift.probe_with_posterior_samples(curv, m, HT, 20)\n",
"m_mean, m_var = ift.probe_with_posterior_samples(curv, 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",
"\n",
"s_data = HT(sh).to_global_data()\n",
"m_data = HT(m).to_global_data()\n",
"m_var_data = m_var.to_global_data()\n",
"d_data = d.to_global_data()\n",
"uncertainty = np.sqrt(np.abs(m_var_data))"
]
},
......@@ -702,10 +650,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",
"sample = HT(curv.draw_sample()+m).to_global_data()\n",
"post_mean = (m_mean + HT(m)).to_global_data()\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, post_mean, sample, s_data - m_data, uncertainty]\n",
"caption = [\"Signal\", \"Reconstruction\", \"Posterior mean\", \"Sample\", \"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",
......
import numpy as np
import nifty4 as ift
from nifty4.library.nonlinearities import Linear, Tanh, Exponential
import numpy as np
np.random.seed(42)
# np.seterr(all="raise",under="ignore")
def adjust_zero_mode(m0, t0):
mtmp = m0.to_global_data().copy()
zero_position = len(m0.shape)*(0,)
zero_mode = mtmp[zero_position]
mtmp[zero_position] = zero_mode / abs(zero_mode)
ttmp = t0.to_global_data().copy()
ttmp[0] += 2 * np.log(abs(zero_mode))
return (ift.Field.from_global_data(m0.domain, mtmp),
ift.Field.from_global_data(t0.domain, ttmp))
if __name__ == "__main__":
noise_level = 1.
p_spec = (lambda k: (.5 / (k + 1) ** 3))
nonlinearity = Linear()
# Set up position space
s_space = ift.RGSpace([128, 128])
# s_space = ift.HPSpace(32)
s_space = ift.RGSpace((128, 128))
h_space = s_space.get_default_codomain()
# 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, target=s_space)
# Set up power space
# Setting up power space
p_space = ift.PowerSpace(h_space,
binbounds=ift.PowerSpace.useful_binbounds(
h_space, logarithmic=True))
# Choose the prior correlation structure and defining correlation operator
p_spec = (lambda k: (.5 / (k + 1) ** 3))
S = ift.create_power_operator(h_space, power_spectrum=p_spec)
# Choosing the prior correlation structure and defining
# correlation operator
p = ift.PS_field(p_space, p_spec)
log_p = ift.log(p)
S = ift.create_power_operator(h_space, power_spectrum=lambda k: 1e-5)
# Draw a sample sh from the prior distribution in harmonic space
sp = ift.PS_field(p_space, p_spec)
sh = ift.power_synthesize(sp, real_signal=True)
# Drawing a sample sh from the prior distribution in harmonic space
sh = S.draw_sample()
# Choose the measurement instrument
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
Instrument = ift.DiagonalOperator(ift.Field(s_space, 1.))
# Instrument._diagonal.val[200:400, 200:400] = 0
# Instrument._diagonal.val[64:512-64, 64:512-64] = 0
# Add a harmonic transformation to the instrument
R = Instrument*fft
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)
# Create mock data
d = R(sh) + 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),
binbounds=p_space.binbounds))
d_data = d.val
ift.plot(d, name="data.png")
mask = np.ones(s_space.shape)
# mask[6000:8000] = 0.
mask[30:70, 30:70] = 0.
mask = ift.Field.from_global_data(s_space, mask)
MaskOperator = ift.DiagonalOperator(mask)
R = ift.GeometryRemover(s_space)
R = R*MaskOperator
# R = R*HT
# R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0,
# response_sigma)
MeasurementOperator = R
d_space = MeasurementOperator.target
Distributor = ift.PowerDistributor(target=h_space, power_space=p_space)
power = Distributor(ift.exp(0.5*log_p))
# Creating the mock data
true_sky = nonlinearity(HT(power*sh))
noiseless_data = MeasurementOperator(true_sky)
noise_amplitude = noiseless_data.val.std()*noise_level
N = ift.ScalingOperator(noise_amplitude**2, d_space)
n = N.draw_sample()
# Creating the mock data
d = noiseless_data + n
m0 = ift.Field.full(h_space, 1e-7)