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

Commit 1046cb1d authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'NIFTy_5' into adjust_variances

parents 5c4d0aba c26a26ea
Pipeline #47210 passed with stages
in 8 minutes and 12 seconds
......@@ -14,6 +14,7 @@ RUN apt-get update && apt-get install -y \
# more optional NIFTy dependencies
&& pip3 install pyfftw \
&& pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
&& pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git \
&& pip3 install jupyter \
&& rm -rf /var/lib/apt/lists/*
......
......@@ -52,6 +52,8 @@ Optional dependencies:
- [pyFFTW](https://pypi.python.org/pypi/pyFFTW) for faster Fourier transforms
- [pyHealpix](https://gitlab.mpcdf.mpg.de/ift/pyHealpix) (for harmonic
transforms involving domains on the sphere)
- [nifty_gridder](https://gitlab.mpcdf.mpg.de/ift/nifty_gridder) (for radio
interferometry responses)
- [mpi4py](https://mpi4py.scipy.org) (for MPI-parallel execution)
- [matplotlib](https://matplotlib.org/) (for field plotting)
......@@ -97,6 +99,10 @@ Support for spherical harmonic transforms is added via:
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
Support for the radio interferometry gridder is added via:
pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
MPI support is added via:
sudo apt-get install python3-mpi4py
......
from time import time
import matplotlib.pyplot as plt
import numpy as np
import nifty5 as ift
ift.fft.enable_fftw()
np.random.seed(40)
N0s, a0s, b0s, c0s = [], [], [], []
N1s, a1s, b1s, c1s = [], [], [], []
for ii in range(10, 23):
nu = 1024
nv = 1024
N = int(2**ii)
print('N = {}'.format(N))
uv = np.random.rand(N, 2) - 0.5
vis = np.random.randn(N) + 1j*np.random.randn(N)
uvspace = ift.RGSpace((nu, nv))
visspace = ift.UnstructuredDomain(N)
img = np.random.randn(nu*nv)
img = img.reshape((nu, nv))
img = ift.from_global_data(uvspace, img)
t0 = time()
GM = ift.GridderMaker(uvspace, eps=1e-7)
idx = GM.getReordering(uv)
uv = uv[idx]
vis = vis[idx]
vis = ift.from_global_data(visspace, vis)
op = GM.getFull(uv).adjoint
t1 = time()
op(img).to_global_data()
t2 = time()
op.adjoint(vis).to_global_data()
t3 = time()
N0s.append(N)
a0s.append(t1 - t0)
b0s.append(t2 - t1)
c0s.append(t3 - t2)
t0 = time()
op = ift.NFFT(uvspace, uv)
t1 = time()
op(img).to_global_data()
t2 = time()
op.adjoint(vis).to_global_data()
t3 = time()
N1s.append(N)
a1s.append(t1 - t0)
b1s.append(t2 - t1)
c1s.append(t3 - t2)
print('Measure rest operator')
sc = ift.StatCalculator()
op = GM.getRest().adjoint
for _ in range(10):
t0 = time()
res = op(img)
sc.add(time() - t0)
t_fft = sc.mean
print('FFT shape', res.shape)
plt.scatter(N0s, a0s, label='Gridder mr')
plt.scatter(N1s, a1s, marker='^', label='NFFT')
plt.legend()
# no idea why this is necessary, but if it is omitted, the range is wrong
plt.ylim(min(a0s+a1s), max(a0s+a1s))
plt.ylabel('time [s]')
plt.title('Initialization')
plt.loglog()
plt.savefig('bench0.png')
plt.close()
plt.scatter(N0s, b0s, color='k', marker='^', label='Gridder mr times')
plt.scatter(N1s, b1s, color='r', marker='^', label='NFFT times')
plt.scatter(N0s, c0s, color='k', label='Gridder mr adjoint times')
plt.scatter(N1s, c1s, color='r', label='NFFT adjoint times')
plt.axhline(sc.mean, label='FFT')
plt.axhline(sc.mean + np.sqrt(sc.var))
plt.axhline(sc.mean - np.sqrt(sc.var))
plt.legend()
plt.ylabel('time [s]')
plt.title('Apply')
plt.loglog()
plt.savefig('bench1.png')
plt.close()
......@@ -45,13 +45,6 @@ def make_random_mask():
return mask.to_global_data()
def mask_to_nan(mask, field):
# Set masked pixels to nan for plotting
masked_data = field.local_data.copy()
masked_data[mask.local_data == 0] = np.nan
return ift.from_local_data(field.domain, masked_data)
if __name__ == '__main__':
np.random.seed(42)
......@@ -64,7 +57,7 @@ if __name__ == '__main__':
if mode == 0:
# One-dimensional regular grid
position_space = ift.RGSpace([1024])
mask = np.ones(position_space.shape)
mask = np.zeros(position_space.shape)
elif mode == 1:
# Two-dimensional regular grid with checkerboard mask
position_space = ift.RGSpace([128, 128])
......@@ -101,23 +94,22 @@ if __name__ == '__main__':
# Build instrument response consisting of a discretization, mask
# and harmonic transformaion
# Data is defined on a geometry-free space, thus the geometry is removed
GR = ift.GeometryRemover(position_space)
# Masking operator to model that parts of the field have not been observed
mask = ift.Field.from_global_data(position_space, mask)
Mask = ift.DiagonalOperator(mask)
Mask = ift.MaskOperator(mask)
# The response operator consists of
# - an harmonic transform (to get to image space)
# - a harmonic transform (to get to image space)
# - the application of the mask
# - the removal of geometric information
# The removal of geometric information is included in the MaskOperator
# it can also be implemented with a GeometryRemover
# Operators can be composed either with parenthesis
R = GR(Mask(HT))
R = Mask(HT)
# or with @
R = GR @ Mask @ HT
R = Mask @ HT
data_space = GR.target
data_space = R.target
# Set the noise covariance N
noise = 5.
......@@ -144,16 +136,17 @@ if __name__ == '__main__':
filename = "getting_started_1_mode_{}.png".format(mode)
if rg and len(position_space.shape) == 1:
plot.add(
[HT(MOCK_SIGNAL), GR.adjoint(data),
[HT(MOCK_SIGNAL), Mask.adjoint(data),
HT(m)],
label=['Mock signal', 'Data', 'Reconstruction'],
alpha=[1, .3, 1])
plot.add(mask_to_nan(mask, HT(m - MOCK_SIGNAL)), title='Residuals')
plot.add(Mask.adjoint(Mask(HT(m - MOCK_SIGNAL))), title='Residuals')
plot.output(nx=2, ny=1, xsize=10, ysize=4, name=filename)
else:
plot.add(HT(MOCK_SIGNAL), title='Mock Signal')
plot.add(mask_to_nan(mask, (GR(Mask)).adjoint(data)), title='Data')
plot.add(Mask.adjoint(data), title='Data')
plot.add(HT(m), title='Reconstruction')
plot.add(mask_to_nan(mask, HT(m - MOCK_SIGNAL)), title='Residuals')
plot.add(Mask.adjoint(Mask(HT(m) - HT(MOCK_SIGNAL))),
title='Residuals')
plot.output(nx=2, ny=2, xsize=10, ysize=10, name=filename)
print("Saved results as '{}'.".format(filename))
This diff is collapsed.
......@@ -35,6 +35,24 @@ Support for spherical harmonic transforms is added via::
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
Support for the radio interferometry gridder is added via:
pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
MPI support is added via::
sudo apt-get install python3-mpi4py
NIFTy documentation is provided by Sphinx. To build the documentation::
sudo apt-get install python3-sphinx-rtd-theme dvipng
cd <nifty_directory>
sh docs/generate.sh
To view the documentation in firefox::
firefox docs/build/index.html
(Note: Make sure that you reinstall nifty after each change since sphinx
imports nifty from the Python path.)
Discretization and Volume in NIFTy
Discretisation and Volume in NIFTy
==================================
.. note:: Some of this discussion is rather technical and may be skipped in a first read-through.
......@@ -160,15 +160,21 @@ Often, log-likelihoods contain integrals over the quantity of interest :math:`s`
\int_\Omega \text{d}x\, s(x) \approx \sum_i s^i\int_{\Omega_i}\text{d}x\, 1
Here the domain of the integral :math:`\Omega = \dot{\bigcup_q} \; \Omega_i` is the disjoint union over smaller :math:`\Omega_i`, e.g. the pixels of the space, and :math:`s_i` is the discretized field value on the :math:`i`-th pixel.
Here the domain of the integral :math:`\Omega = \dot{\bigcup_q} \; \Omega_i` is the disjoint union over smaller :math:`\Omega_i`, e.g. the pixels of the space, and :math:`s_i` is the discretised field value on the :math:`i`-th pixel.
This introduces the weighting :math:`V_i=\int_{\Omega_i}\text{d}x\, 1`, also called the volume factor, a property of the space.
NIFTy aids you in constructing your own log-likelihood by providing methods like :func:`~field.Field.weight`, which weights all pixels of a field with their corresponding volume.
An integral over a :class:`~field.Field` :code:`s` can be performed by calling :code:`s.weight(1).sum()`, which is equivalent to :code:`s.integrate()`.
Volume factors are also applied automatically in the following places:
- :class:`~operators.harmonic_operators.FFTOperator` as well as all other harmonic operators. Here the zero mode of the transformed field is the integral over the original field, thus the whole field is weighted once.
- some response operators, such as the :class:`~library.los_response.LOSResponse`. In this operator a line integral is descritized, so a 1-dimensional volume factor is applied.
- In :class:`~library.correlated_fields.CorrelatedField` as well :class:`~library.correlated_fields.MfCorrelatedField`, the field is multiplied by the square root of the total volume in configuration space. This ensures that the same field reconstructed over a larger domain has the same variance in position space in the limit of infinite resolution. It also ensures that power spectra in NIFTy behave according to the definition of a power spectrum, namely the power of a k-mode is the expectation of the k-mode square, divided by the volume of the space.
- :class:`~operators.harmonic_operators.FFTOperator` as well as all other harmonic operators.
Here the zero mode of the transformed field is the integral over the original field, thus the whole field is weighted once.
- Some response operators, such as the :class:`~library.los_response.LOSResponse`.
In this operator a line integral is discretised, so a 1-dimensional volume factor is applied.
- In :class:`~library.correlated_fields.CorrelatedField` as well as :class:`~library.correlated_fields.MfCorrelatedField`.
Both describe fields with a smooth, a priori unknown correlation structure specified by a power spectrum.
The field is multiplied by the square root of the total volume of it domain's harmonic counterpart.
This ensures that the same power spectrum can be used regardless of the chosen resolution, provided the total volume of the space remains the same.
It also guarantees that the power spectra in NIFTy behave according to their definition, i.e. the power of a mode :math:`s_k` is the expectation value of that mode squared, divided by the volume of its space :math:`P(k) = \left\langle s_k^2 \right\rangle / V_k`.
Note that in contrast to some older versions of NIFTy, the dot product :code:`s.vdot(t)` of fields does **not** apply a volume factor, but instead just sums over the field components,
......
......@@ -74,7 +74,8 @@ from .minimization.metric_gaussian_kl import MetricGaussianKL
from .sugar import *
from .plot import Plot
from .library.smooth_linear_amplitude import SLAmplitude, CepstrumOperator
from .library.smooth_linear_amplitude import (
SLAmplitude, LinearSLAmplitude, CepstrumOperator)
from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse
from .library.dynamic_operator import (dynamic_operator,
......@@ -86,6 +87,7 @@ from .library.correlated_fields import CorrelatedField, MfCorrelatedField
from .library.adjust_variances import (make_adjust_variances_hamiltonian,
do_adjust_variances)
from .library.nfft import NFFT
from .library.gridder import GridderMaker
from . import extra
......
......@@ -149,7 +149,7 @@ def ensure_default_distributed(arr):
def absmax(arr):
return np.linalg.norm(arr.rehape(-1), ord=np.inf)
return np.linalg.norm(arr.reshape(-1), ord=np.inf)
def norm(arr, ord=2):
......
......@@ -33,7 +33,8 @@ def _assert_allclose(f1, f2, atol, rtol):
_assert_allclose(val, f2[key], atol=atol, rtol=rtol)
def _adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol):
def _adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol,
only_r_linear):
needed_cap = op.TIMES | op.ADJOINT_TIMES
if (op.capability & needed_cap) != needed_cap:
return
......@@ -41,6 +42,8 @@ def _adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol):
f2 = from_random("normal", op.target, dtype=target_dtype)
res1 = f1.vdot(op.adjoint_times(f2))
res2 = op.times(f1).vdot(f2)
if only_r_linear:
res1, res2 = res1.real, res2.real
np.testing.assert_allclose(res1, res2, atol=atol, rtol=rtol)
......@@ -57,8 +60,10 @@ def _inverse_implementation(op, domain_dtype, target_dtype, atol, rtol):
_assert_allclose(res, foo, atol=atol, rtol=rtol)
def _full_implementation(op, domain_dtype, target_dtype, atol, rtol):
_adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol)
def _full_implementation(op, domain_dtype, target_dtype, atol, rtol,
only_r_linear):
_adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol,
only_r_linear)
_inverse_implementation(op, domain_dtype, target_dtype, atol, rtol)
......@@ -72,7 +77,7 @@ def _check_linearity(op, domain_dtype, atol, rtol):
def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
atol=0, rtol=1e-7):
atol=0, rtol=1e-7, only_r_linear=False):
"""
Checks an operator for algebraic consistency of its capabilities.
......@@ -98,15 +103,21 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
Relative tolerance for the check. If atol is specified,
then satisfying any tolerance will let the check pass.
Default: 0.
only_r_linear: bool
set to True if the operator is only R-linear, not C-linear.
This will relax the adjointness test accordingly.
"""
if not isinstance(op, LinearOperator):
raise TypeError('This test tests only linear operators.')
_check_linearity(op, domain_dtype, atol, rtol)
_full_implementation(op, domain_dtype, target_dtype, atol, rtol)
_full_implementation(op.adjoint, target_dtype, domain_dtype, atol, rtol)
_full_implementation(op.inverse, target_dtype, domain_dtype, atol, rtol)
_full_implementation(op, domain_dtype, target_dtype, atol, rtol,
only_r_linear)
_full_implementation(op.adjoint, target_dtype, domain_dtype, atol, rtol,
only_r_linear)
_full_implementation(op.inverse, target_dtype, domain_dtype, atol, rtol,
only_r_linear)
_full_implementation(op.adjoint.inverse, domain_dtype, target_dtype, atol,
rtol)
rtol, only_r_linear)
def _get_acceptable_location(op, loc, lin):
......
# 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) 2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from ..domain_tuple import DomainTuple
from ..domains.rg_space import RGSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..fft import hartley
from ..operators.linear_operator import LinearOperator
from ..sugar import from_global_data, makeDomain
class GridderMaker(object):
def __init__(self, domain, eps=1e-15):
domain = makeDomain(domain)
if (len(domain) != 1 or not isinstance(domain[0], RGSpace) or
not len(domain.shape) == 2):
raise ValueError("need domain with exactly one 2D RGSpace")
nu, nv = domain.shape
if nu % 2 != 0 or nv % 2 != 0:
raise ValueError("dimensions must be even")
rat = 3 if eps < 1e-11 else 2
nu2, nv2 = rat*nu, rat*nv
nspread = int(-np.log(eps)/(np.pi*(rat-1)/(rat-.5)) + .5) + 1
nu2 = max([nu2, 2*nspread])
nv2 = max([nv2, 2*nspread])
r2lamb = rat*rat*nspread/(rat*(rat-.5))
oversampled_domain = RGSpace(
[nu2, nv2], distances=[1, 1], harmonic=False)
self._nspread = nspread
self._r2lamb = r2lamb
self._rest = _RestOperator(domain, oversampled_domain, r2lamb)
def getReordering(self, uv):
from nifty_gridder import peanoindex
nu2, nv2 = self._rest._domain.shape
return peanoindex(uv, nu2, nv2)
def getGridder(self, uv):
return RadioGridder(self._rest.domain, self._nspread, self._r2lamb, uv)
def getRest(self):
return self._rest
def getFull(self, uv):
return self.getRest() @ self.getGridder(uv)
class _RestOperator(LinearOperator):
def __init__(self, domain, oversampled_domain, r2lamb):
self._domain = makeDomain(oversampled_domain)
self._target = domain
nu, nv = domain.shape
nu2, nv2 = oversampled_domain.shape
# compute deconvolution operator
rng = np.arange(nu)
k = np.minimum(rng, nu-rng)
c = np.pi*r2lamb/nu2**2
self._deconv_u = np.roll(np.exp(c*k**2), -nu//2).reshape((-1, 1))
rng = np.arange(nv)
k = np.minimum(rng, nv-rng)
c = np.pi*r2lamb/nv2**2
self._deconv_v = np.roll(
np.exp(c*k**2)/r2lamb, -nv//2).reshape((1, -1))
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
self._check_input(x, mode)
nu, nv = self._target.shape
res = x.to_global_data()
if mode == self.TIMES:
res = hartley(res)
res = np.roll(res, (nu//2, nv//2), axis=(0, 1))
res = res[:nu, :nv]
res *= self._deconv_u
res *= self._deconv_v
else:
res = res*self._deconv_u
res *= self._deconv_v
nu2, nv2 = self._domain.shape
res = np.pad(res, ((0, nu2-nu), (0, nv2-nv)), mode='constant',
constant_values=0)
res = np.roll(res, (-nu//2, -nv//2), axis=(0, 1))
res = hartley(res)
return from_global_data(self._tgt(mode), res)
class RadioGridder(LinearOperator):
def __init__(self, target, nspread, r2lamb, uv):
self._domain = DomainTuple.make(
UnstructuredDomain((uv.shape[0],)))
self._target = DomainTuple.make(target)
self._capability = self.TIMES | self.ADJOINT_TIMES
self._nspread, self._r2lamb = int(nspread), float(r2lamb)
self._uv = uv # FIXME: should we write-protect this?
def apply(self, x, mode):
from nifty_gridder import (to_grid, to_grid_post,
from_grid, from_grid_pre)
self._check_input(x, mode)
nu2, nv2 = self._target.shape
x = x.to_global_data()
if mode == self.TIMES:
res = to_grid(self._uv, x, nu2, nv2, self._nspread, self._r2lamb)
res = to_grid_post(res)
else:
x = from_grid_pre(x)
res = from_grid(self._uv, x, nu2, nv2, self._nspread, self._r2lamb)
return from_global_data(self._tgt(mode), res)
......@@ -169,6 +169,18 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
which returns on its target a power spectrum which consists out of a
smooth and a linear part.
'''
return LinearSLAmplitude(target=target, n_pix=n_pix, a=a, k0=k0, sm=sm,
sv=sv, im=im, iv=iv, keys=keys).exp()
def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
keys=['tau', 'phi']):
'''LinearOperator for parametrizing smooth log-amplitudes (square roots of
power spectra).
Logarithm of SLAmplitude
See documentation of SLAmplitude for more details
'''
if not (isinstance(n_pix, int) and isinstance(target, PowerSpace)):
raise TypeError
......@@ -196,4 +208,4 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
loglog_ampl = 0.5*(smooth + linear)
# Go from loglog-space to linear-linear-space
return et @ loglog_ampl.exp()
return et @ loglog_ampl
......@@ -45,6 +45,13 @@ def FuncConvolutionOperator(domain, func, space=None):
The index of the subdomain on which the operator should act
If None, it is set to 0 if `domain` contains exactly one space.
`domain[space]` must be of type `RGSpace`, `HPSpace`, or `GLSpace`.
Notes
-----
The operator assumes periodic boundaries in the input domain. This means
for a sufficiently broad function a point source close to the boundary will
blur into the opposite side of the image. Zero padding can be applied to
avoid this behaviour.
"""
domain = DomainTuple.make(domain)
space = utilities.infer_space(domain, space)
......
......@@ -20,6 +20,7 @@ import numpy as np
from .. import utilities
from ..domain_tuple import DomainTuple
from ..field import Field
from ..multi_field import MultiField
from ..linearization import Linearization
from ..sugar import makeDomain, makeOp
from .linear_operator import LinearOperator
......@@ -121,7 +122,7 @@ class GaussianEnergy(EnergyOperator):
"""
def __init__(self, mean=None, covariance=None, domain=None):
if mean is not None and not isinstance(mean, Field):
if mean is not None and not isinstance(mean, (Field, MultiField)):
raise TypeError
if covariance is not None and not isinstance(covariance,
LinearOperator):
......@@ -307,7 +308,6 @@ class StandardHamiltonian(EnergyOperator):
Tells an internal :class:`SamplingEnabler` which convergence criterion
to use to draw Gaussian samples.
See also
--------
`Encoding prior knowledge in the structure of the likelihood`,
......
......@@ -361,6 +361,13 @@ class HarmonicTransformOperator(LinearOperator):
The index of the domain on which the operator should act
If None, it is set to 0 if domain contains exactly one subdomain.
domain[space] must be a harmonic domain.
Notes
-----
HarmonicTransformOperator uses a Hartley transformation to transform
between harmonic and non-harmonic RGSpaces. This has the advantage that all
field values are real in either space. If you require a true Fourier
transform you should use FFTOperator instead.
"""
def __init__(self, domain, target=None, space=None):
......
......@@ -384,7 +384,6 @@ class _OpSum(Operator):
v = x._val if lin else x
v1 = v.extract(self._op1.domain)
v2 = v.extract(self._op2.domain)
res = None
if not lin:
return self._op1(v1).unite(self._op2(v2))
wm = x.want_metric
......@@ -393,7 +392,10 @@ class _OpSum(Operator):
op = lin1._jac._myadd(lin2._jac, False)
res = lin1.new(lin1._val.unite(lin2._val), op(x.jac))
if lin1._metric is not None and lin2._metric is not None:
res = res.add_metric(lin1._metric + lin2._metric)
from .sandwich_operator import SandwichOperator
met = lin1._metric._myadd(lin2._metric, False)
met = SandwichOperator.make(x.jac, met)
res = res.add_metric(met)
return res
def _simplify_for_constant_input_nontrivial(self, c_inp):
......
......@@ -66,9 +66,12 @@ class ScalingOperator(EndomorphicOperator):
return x
if fct == 0.:
return full(self.domain, 0.)
if (mode & 10) != 0:
MODES_WITH_ADJOINT = self.ADJOINT_TIMES | self.ADJOINT_INVERSE_TIMES
MODES_WITH_INVERSE = self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES
if (mode & MODES_WITH_ADJOINT) != 0:
fct = np.conj(fct)
if (mode & 12) != 0:
if (mode & MODES_WITH_INVERSE) != 0:
fct = 1./fct
return x*fct
......
......@@ -29,6 +29,7 @@ from .multi_field import MultiField
from .operators.block_diagonal_operator import BlockDiagonalOperator
from .operators.diagonal_operator import DiagonalOperator
from .operators.distributors import PowerDistributor
from .plot import Plot
__all__ = ['PS_field', 'power_analyze', 'create_power_operator',
'create_harmonic_smoothing_operator', 'from_random',
......@@ -37,7 +38,7 @@ __all__ = ['PS_field', 'power_analyze', 'create_power_operator',
'sin', 'cos', 'tan', 'sinh', 'cosh',
'absolute', 'one_over', 'clip', 'sinc',
'conjugate', 'get_signal_variance', 'makeOp', 'domain_union',
'get_default_codomain']
'get_default_codomain', 'single_plot']
def PS_field(pspace, func):
......@@ -434,3 +435,14 @@ def get_default_codomain(domainoid, space=None):
ret = [dom for dom in domainoid]
ret[space] = domainoid[space].get_default_codomain()
return DomainTuple.make(ret)
def single_plot(field, **kwargs):
"""Creates a single plot using `Plot`.
Keyword arguments are passed to both `Plot.add` and `Plot.output`.
"""
p = Plot()
p.add(field, **kwargs)
if 'title' in kwargs:
del(kwargs['title'])
p.output(**kwargs)
......@@ -71,6 +71,20 @@ def testLinearInterpolator():
ift.extra.consistency_check(op)
@pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
def testRealizer(sp):
op = ift.Realizer(sp)
ift.extra.consistency_check(op, np.complex128, np.float64,
only_r_linear=True)
@pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
def testConjugationOperator(sp):
op = ift.ConjugationOperator(sp)
ift.extra.consistency_check(op, np.complex128, np.complex128,
only_r_linear=True)
@pmp('args', [(ift.RGSpace(10, harmonic=True), 4, 0), (ift.RGSpace(
(24, 31), distances=(0.4, 2.34), harmonic=True), 3, 0),
(ift.LMSpace(4), 10, 0)])
......
# 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) 2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
import pytest
from numpy.testing import assert_allclose
import nifty5 as ift
np.random.seed(40)
pmp = pytest.mark.parametrize
@pmp('nu', [12, 128])
@pmp('nv', [4, 12, 128])
@pmp('N', [1, 10, 100])
def test_gridding(nu, nv, N):
uv = np.random.rand(N, 2) - 0.5
vis = np.random.randn(N) + 1j*np.random.randn(N)
# Nifty
GM = ift.GridderMaker(ift.RGSpace((nu, nv)))
# re-order for performance
idx = GM.getReordering(uv)
uv, vis = uv[idx], vis[idx]
vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis)
Op = GM.getFull(uv)
pynu = Op(vis2).to_global_data()
# DFT
x, y = np.meshgrid(
*[-ss/2 + np.arange(ss) for ss in [nu, nv]], indexing='ij')
dft = pynu*0.
for i in range(N):
dft += (vis[i]*np.exp(2j*np.pi*(x*uv[i, 0] + y*uv[i, 1]))).real
assert_allclose(dft, pynu)
@pmp('eps', [1e-2, 1e-6, 1e-15])
@pmp('nu', [12, 128])
@pmp('nv', [4, 12, 128])
@pmp('N', [1, 10, 100])
def test_build(nu, nv, N, eps):
dom = ift.RGSpace([nu, nv])
uv = np.random.rand(N, 2) - 0.5
GM = ift.GridderMaker(dom)
# re-order for performance
idx = GM.getReordering(uv)
<