Commit b82dc0dc authored by Martin Reinecke's avatar Martin Reinecke

Merge remote-tracking branch 'origin/NIFTy_5' into improving_consistency_checks

parents 5655533a 35dc4f9b
Pipeline #47425 passed with stages
in 9 minutes and 40 seconds
......@@ -10,10 +10,11 @@ RUN apt-get update && apt-get install -y \
# Testing dependencies
python3-pytest-cov jupyter \
# Optional NIFTy dependencies
libfftw3-dev python3-mpi4py python3-matplotlib \
libfftw3-dev python3-mpi4py python3-matplotlib python3-pynfft \
# 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)
......@@ -60,7 +62,7 @@ Optional dependencies:
The current version of Nifty5 can be obtained by cloning the repository and
switching to the NIFTy_5 branch:
git clone https://gitlab.mpcdf.mpg.de/ift/NIFTy.git
git clone https://gitlab.mpcdf.mpg.de/ift/nifty.git
### Installation
......@@ -70,7 +72,7 @@ distributions, the "apt" lines will need slight changes.
NIFTy5 and its mandatory dependencies can be installed via:
sudo apt-get install git python3 python3-pip python3-dev
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_5
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_5
Plotting support is added via:
......@@ -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))
import numpy as np
import nifty5 as ift
def convtest(test_signal, delta, func):
domain = test_signal.domain
# Create Convolution Operator
conv_op = ift.FuncConvolutionOperator(domain, func)
# Convolve, Adjoint-Convolve
conv_signal = conv_op(test_signal)
cac_signal = conv_op.adjoint_times(conv_signal)
print(test_signal.integrate(), conv_signal.integrate(),
cac_signal.integrate())
# generate kernel image
conv_delta = conv_op(delta)
# Plot results
plot = ift.Plot()
plot.add(signal, title='Signal')
plot.add(conv_signal, title='Signal Convolved')
plot.add(cac_signal, title='Signal, Conv, Adj-Conv')
plot.add(conv_delta, title='Kernel')
plot.output()
# Healpix test
nside = 64
npix = 12 * nside * nside
domain = ift.HPSpace(nside)
# Define test signal (some point sources)
signal_vals = np.zeros(npix, dtype=np.float64)
for i in range(0, npix, npix//12 + 27):
signal_vals[i] = 500.
signal = ift.from_global_data(domain, signal_vals)
delta_vals = np.zeros(npix, dtype=np.float64)
delta_vals[0] = 1.0
delta = ift.from_global_data(domain, delta_vals)
# Define kernel function
def func(theta):
ct = np.cos(theta)
return 1. * np.logical_and(ct > 0.7, ct <= 0.8)
convtest(signal, delta, func)
domain = ift.RGSpace((100, 100))
# Define test signal (some point sources)
signal_vals = np.zeros(domain.shape, dtype=np.float64)
signal_vals[35, 70] = 5000.
signal_vals[45, 8] = 5000.
signal = ift.from_global_data(domain, signal_vals)
# Define delta signal, generate kernel image
delta_vals = np.zeros(domain.shape, dtype=np.float64)
delta_vals[0, 0] = 1.0
delta = ift.from_global_data(domain, delta_vals)
# Define kernel function
def func(dist):
return 1. * np.logical_and(dist > 0.1, dist <= 0.2)
convtest(signal, delta, func)
This diff is collapsed.
......@@ -37,7 +37,7 @@ NIFTy comes with reimplemented MAP and VI estimators.
Free Theory & Implicit Operators
--------------------------------
A free IFT appears when the signal field :math:`{s}` and the noise :math:`{n}` of the data :math:`{d}` are independent, zero-centered Gaussian processes of kown covariances :math:`{S}` and :math:`{N}`, respectively,
A free IFT appears when the signal field :math:`{s}` and the noise :math:`{n}` of the data :math:`{d}` are independent, zero-centered Gaussian processes of known covariances :math:`{S}` and :math:`{N}`, respectively,
.. math::
......@@ -49,7 +49,7 @@ and the measurement equation is linear in both signal and noise,
d= R\, s + n,
with :math:`{R}` being the measurement response, which maps the continous signal field into the discrete data space.
with :math:`{R}` being the measurement response, which maps the continuous signal field into the discrete data space.
This is called a free theory, as the information Hamiltonian
......@@ -96,8 +96,8 @@ These implicit operators can be combined into new operators, e.g. to :math:`{D^{
The invocation of an inverse operator applied to a vector might trigger the execution of a numerical linear algebra solver.
Thus, when NIFTy calculates :math:`{m = D\, j}`, it actually solves :math:`{D^{-1} m = j}` for :math:`{m}` behind the scenes.
The advantage of implicit operators to explicit matrices is the reduced memory requirements.
The reconstruction of only a Megapixel image would otherwithe require the storage and processing of matrices with sizes of several Terabytes.
The advantage of implicit operators compared to explicit matrices is the reduced memory consumption;
for the reconstruction of just a Megapixel image the latter would already require several Terabytes.
Larger images could not be dealt with due to the quadratic memory requirements of explicit operator representations.
The demo codes `demos/getting_started_1.py` and `demos/Wiener_Filter.ipynb` illustrate this.
......@@ -106,7 +106,7 @@ The demo codes `demos/getting_started_1.py` and `demos/Wiener_Filter.ipynb` illu
Generative Models
-----------------
For more sophisticated measurement situations, involving non-linear measuremnts, unknown covariances, calibration constants and the like, it is recommended to formulate those as generative models for which NIFTy provides powerful inference algorithms.
For more sophisticated measurement situations (involving non-linear measurements, unknown covariances, calibration constants and the like) it is recommended to formulate those as generative models for which NIFTy provides powerful inference algorithms.
In a generative model, all known or unknown quantities are described as the results of generative processes, which start with simple probability distributions, like the uniform, the i.i.d. Gaussian, or the delta distribution.
......@@ -129,7 +129,7 @@ NIFTy takes advantage of this formulation in several ways:
1) All prior degrees of freedom have unit covariance, which improves the condition number of operators that need to be inverted.
2) The amplitude operator can be regarded as part of the response, :math:`{R'=R\,A}`.
In general, more sophisticated responses can be constructed out of the composition of simpler operators.
In general, more sophisticated responses can be obtained by combining simpler operators.
3) The response can be non-linear, e.g. :math:`{R'(s)=R \exp(A\,\xi)}`, see `demos/getting_started_2.py`.
......@@ -168,7 +168,7 @@ Maximum a Posteriori
One popular field estimation method is Maximum a Posteriori (MAP).
It only requires to minimize the information Hamiltonian, e.g by a gradient descent method that stops when
It only requires minimizing the information Hamiltonian, e.g. by a gradient descent method that stops when
.. math::
......
......@@ -8,7 +8,7 @@ distributions, the "apt" lines will need slight changes.
NIFTy5 and its mandatory dependencies can be installed via::
sudo apt-get install git python3 python3-pip python3-dev
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_5
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_5
Plotting support is added via::
......@@ -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,
......
......@@ -50,6 +50,7 @@ from .operators.value_inserter import ValueInserter
from .operators.energy_operators import (
EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
BernoulliEnergy, StandardHamiltonian, AveragedEnergy)
from .operators.convolution_operators import FuncConvolutionOperator
from .probing import probe_with_posterior_samples, probe_diagonal, \
StatCalculator
......@@ -73,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,
......@@ -84,6 +86,8 @@ from .library.wiener_filter_curvature import WienerFilterCurvature
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
......@@ -93,5 +97,10 @@ from .logger import logger
from .linearization import Linearization
from . import internal_config
_scheme = internal_config.parallelization_scheme()
if _scheme == "Samples":
from .minimization.metric_gaussian_kl_mpi import MetricGaussianKL_MPI
# We deliberately don't set __all__ here, because we don't want people to do a
# "from nifty5 import *"; that would swamp the global namespace.
......@@ -235,7 +235,6 @@ class data_object(object):
for op in ["__add__", "__radd__", "__iadd__",
"__sub__", "__rsub__", "__isub__",
"__mul__", "__rmul__", "__imul__",
"__div__", "__rdiv__", "__idiv__",
"__truediv__", "__rtruediv__", "__itruediv__",
"__floordiv__", "__rfloordiv__", "__ifloordiv__",
"__pow__", "__rpow__", "__ipow__",
......
......@@ -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):
......
......@@ -15,11 +15,19 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
try:
from mpi4py import MPI
if MPI.COMM_WORLD.Get_size() == 1:
from .data_objects.numpy_do import *
else:
from .data_objects.distributed_do import *
except ImportError:
from . import internal_config
_scheme = internal_config.parallelization_scheme()
if _scheme in ("Samples", "None"):
from .data_objects.numpy_do import *
else:
try:
from mpi4py import MPI
if MPI.COMM_WORLD.Get_size() == 1:
from .data_objects.numpy_do import *
else:
from .data_objects.distributed_do import *
except ImportError:
from .data_objects.numpy_do import *
......@@ -103,6 +103,36 @@ class LMSpace(StructuredDomain):
def get_fft_smoothing_kernel_function(self, sigma):
return lambda x: self._kernel(x, sigma)
def get_conv_kernel_from_func(self, func):
"""Creates a convolution kernel defined by a function.
Parameters
----------
func: function
This function needs to take exactly one argument, which is
colatitude in radians, and return the kernel amplitude at that
colatitude.
Assumes the function to be radially symmetric,
e.g. only dependant on theta in radians"""
from .gl_space import GLSpace
from ..operators.harmonic_operators import HarmonicTransformOperator
import pyHealpix
# define azimuthally symmetric spaces for kernel transform
gl = GLSpace(self.lmax + 1, 1)
lm0 = gl.get_default_codomain()
theta = pyHealpix.GL_thetas(gl.nlat)
# evaluate the kernel function at the required thetas
kernel_sphere = Field.from_global_data(gl, func(theta))
# normalize the kernel such that the integral over the sphere is 4pi
kernel_sphere = kernel_sphere * (4 * np.pi / kernel_sphere.integrate())
# compute the spherical harmonic coefficients of the kernel
op = HarmonicTransformOperator(lm0, gl)
kernel_lm = op.adjoint_times(kernel_sphere.weight(1)).to_global_data()
# evaluate the k lengths of the harmonic space
k_lengths = self.get_k_length_array().to_global_data().astype(np.int)
return Field.from_global_data(self, kernel_lm[k_lengths])
@property
def lmax(self):
"""int : maximum allowed :math:`l`
......
......@@ -90,9 +90,7 @@ class RGSpace(StructuredDomain):
def scalar_dvol(self):
return self._dvol
def get_k_length_array(self):
if (not self.harmonic):
raise NotImplementedError
def _get_dist_array(self):
ibegin = dobj.ibegin_from_shape(self._shape)
res = np.arange(self.local_shape[0], dtype=np.float64) + ibegin[0]
res = np.minimum(res, self.shape[0]-res)*self.distances[0]
......@@ -106,6 +104,11 @@ class RGSpace(StructuredDomain):
res = np.add.outer(res, tmp)
return Field.from_local_data(self, np.sqrt(res))
def get_k_length_array(self):
if (not self.harmonic):
raise NotImplementedError
return self._get_dist_array()
def get_unique_k_lengths(self):
if (not self.harmonic):
raise NotImplementedError
......@@ -146,6 +149,27 @@ class RGSpace(StructuredDomain):
raise NotImplementedError
return lambda x: self._kernel(x, sigma)
def get_conv_kernel_from_func(self, func):
"""Creates a convolution kernel defined by a function.
Parameters
----------
func: function
This function needs to take exactly one argument, which is
distance from center (in the same units as the RGSpace distances),
and return the kernel amplitude at that distance.
Assumes the function to be radially symmetric,
e.g. only dependant on distance"""
from ..operators.harmonic_operators import HarmonicTransformOperator
if (not self.harmonic):
raise NotImplementedError
op = HarmonicTransformOperator(self, self.get_default_codomain())
dist = op.target[0]._get_dist_array()
kernel = Field.from_local_data(op.target, func(dist.local_data))
kernel = kernel / kernel.integrate()
return op.adjoint_times(kernel.weight(1))
def get_default_codomain(self):
"""Returns a :class:`RGSpace` object representing the (position or
harmonic) partner domain of `self`, depending on `self.harmonic`.
......
......@@ -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):
......
......@@ -626,6 +626,11 @@ class Field(object):
raise ValueError("domain mismatch")
return self
def extract_part(self, dom):
if dom != self._domain:
raise ValueError("domain mismatch")
return self
def unite(self, other):
return self+other
......@@ -636,6 +641,8 @@ class Field(object):
return 0.5*(1.+self.tanh())
def clip(self, min=None, max=None):
min = min.local_data if isinstance(min, Field) else min
max = max.local_data if isinstance(max, Field) else max
return Field(self._domain, dobj.clip(self._val, min, max))
def one_over(self):
......@@ -656,7 +663,6 @@ class Field(object):
for op in ["__add__", "__radd__",
"__sub__", "__rsub__",
"__mul__", "__rmul__",
"__div__", "__rdiv__",
"__truediv__", "__rtruediv__",
"__floordiv__", "__rfloordiv__",
"__pow__", "__rpow__",
......
# 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.
# Internal configuration switches, typically for experimental features.
# Leave unchanged unless you know what you are doing!
def parallelization_scheme():
return "Standard"
# return "Samples"
# return "None"
......@@ -17,9 +17,8 @@
from ..minimization.energy_adapter import EnergyAdapter
from ..multi_field import MultiField
from ..operators.distributors import PowerDistributor
from ..operators.energy_operators import (StandardHamiltonian,
InverseGammaLikelihood)
from ..operators.energy_operators import (InverseGammaLikelihood,
StandardHamiltonian)
from ..operators.scaling_operator import ScalingOperator
from ..operators.simple_linear_operators import ducktape
......@@ -79,27 +78,23 @@ def make_adjust_variances_hamiltonian(a,
ic_samp=ic_samp)
def do_adjust_variances(position,
</