Commit 52d9f783 authored by Martin Reinecke's avatar Martin Reinecke

merge NIFTy_5

parents acdb0aaf 95837aad
......@@ -5,7 +5,7 @@ RUN apt-get update && apt-get install -y \
git \
# Packages needed for NIFTy
libfftw3-dev \
python3 python3-pip python3-dev python3-future python3-scipy cython3 \
python3 python3-pip python3-dev python3-scipy \
# Documentation build dependencies
python3-sphinx python3-sphinx-rtd-theme \
# Testing dependencies
......
......@@ -39,9 +39,9 @@ Installation
- [Python 3](https://www.python.org/) (3.5.x or later)
- [SciPy](https://www.scipy.org/)
- [pyFFTW](https://pypi.python.org/pypi/pyFFTW)
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)
- [mpi4py](https://mpi4py.scipy.org) (for MPI-parallel execution)
......@@ -61,18 +61,29 @@ distributions, the "apt" lines will need slight changes.
NIFTy5 and its mandatory dependencies can be installed via:
sudo apt-get install git libfftw3-dev python3 python3-pip python3-dev
sudo apt-get install git python3 python3-pip python3-dev
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_5
Plotting support is added via:
pip3 install --user matplotlib
FFTW support is added via:
sudo apt-get install libfftw3-dev
pip3 install --user pyfftw
To actually use FFTW in your Nifty calculations, you need to call
nifty5.fft.enable_fftw()
at the beginning of your code.
(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:
pip3 install --user matplotlib
Support for spherical harmonic transforms is added via:
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
......@@ -86,7 +97,7 @@ MPI support is added via:
To run the tests, additional packages are required:
sudo apt-get install python3-coverage python3-parameterized python3-pytest python3-pytest-cov
sudo apt-get install python3-coverage python3-pytest python3-pytest-cov
Afterwards the tests (including a coverage report) can be run using the
following command in the repository root:
......
......@@ -264,13 +264,13 @@ This functionality is provided by NIFTy's
:class:`~inversion_enabler.InversionEnabler` class, which is itself a linear
operator.
.. currentmodule:: nifty5.operators.linear_operator
.. currentmodule:: nifty5.operators.operator
Direct multiplication and adjoint inverse multiplication transform a field
defined on the operator's :attr:`~LinearOperator.domain` to one defined on the
operator's :attr:`~LinearOperator.target`, whereas adjoint multiplication and
inverse multiplication transform from :attr:`~LinearOperator.target` to
:attr:`~LinearOperator.domain`.
defined on the operator's :attr:`~Operator.domain` to one defined on the
operator's :attr:`~Operator.target`, whereas adjoint multiplication and inverse
multiplication transform from :attr:`~Operator.target` to
:attr:`~Operator.domain`.
.. currentmodule:: nifty5.operators
......@@ -379,7 +379,7 @@ Minimization algorithms
All minimization algorithms in NIFTy inherit from the abstract
:class:`~minimizer.Minimizer` class, which presents a minimalistic interface
consisting only of a :meth:`~minimizer.Minimizer.__call__()` method taking an
consisting only of a :meth:`~minimizer.Minimizer.__call__` method taking an
:class:`~energy.Energy` object and optionally a preconditioning operator, and
returning the energy at the discovered minimum and a status code.
......@@ -399,17 +399,16 @@ Many minimizers for nonlinear problems can be characterized as
This family of algorithms is encapsulated in NIFTy's
:class:`~descent_minimizers.DescentMinimizer` class, which currently has three
concrete implementations: :class:`~descent_minimizers.SteepestDescent`,
:class:`~descent_minimizers.RelaxedNewton`,
:class:`~descent_minimizers.NewtonCG`, :class:`~descent_minimizers.L_BFGS` and
:class:`~descent_minimizers.VL_BFGS`. Of these algorithms, only
:class:`~descent_minimizers.RelaxedNewton` requires the energy object to provide
:class:`~descent_minimizers.NewtonCG` requires the energy object to provide
a :attr:`~energy.Energy.metric` property, the others only need energy values and
gradients.
The flexibility of NIFTy's design allows using externally provided minimizers.
With only small effort, adapters for two SciPy minimizers were written; they are
available under the names :class:`~scipy_minimizer.ScipyCG` and
:class:`L_BFGS_B`.
:class:`~scipy_minimizer.L_BFGS_B`.
Application to operator inversion
......@@ -432,4 +431,4 @@ performs a minimization of a
:class:`~minimization.quadratic_energy.QuadraticEnergy` by means of the
:class:`~minimization.conjugate_gradient.ConjugateGradient` algorithm. An
example is provided in
:func:`~ļibrary.wiener_filter_curvature.WienerFilterCurvature`.
:func:`~library.wiener_filter_curvature.WienerFilterCurvature`.
......@@ -13,6 +13,7 @@ napoleon_use_ivar = True
napoleon_use_admonition_for_notes = True
napoleon_use_admonition_for_examples = True
napoleon_use_admonition_for_references = True
napoleon_include_special_with_doc = True
project = u'NIFTy5'
copyright = u'2013-2019, Max-Planck-Society'
......@@ -27,3 +28,5 @@ add_module_names = False
html_theme = "sphinx_rtd_theme"
html_logo = 'nifty_logo_black.png'
exclude_patterns = ['mod/modules.rst']
This diff is collapsed.
NIFTy -- Numerical Information Field Theory
===========================================
**NIFTy** [1]_, "\ **N**\umerical **I**\nformation **F**\ield **T**\heor\ **y**\ ", is a versatile library designed to enable the development of signal inference algorithms that are independent of the underlying spatial grid and its resolution.
**NIFTy** [1]_, [2]_, "\ **N**\umerical **I**\nformation **F**\ield **T**\heor\ **y**\ ", is a versatile library designed to enable the development of signal inference algorithms that are independent of the underlying spatial grid and its resolution.
Its object-oriented framework is written in Python, although it accesses libraries written in C++ and C for efficiency.
NIFTy offers a toolkit that abstracts discretized representations of continuous spaces, fields in these spaces, and operators acting on fields into classes.
......@@ -13,22 +13,18 @@ The set of spaces on which NIFTy operates comprises point sets, *n*-dimensional
References
----------
.. [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>`_
.. [1] Selig et al., "NIFTY - Numerical Information Field Theory. A versatile PYTHON library for signal inference ", 2013, Astronmy and Astrophysics 554, 26; `[DOI] <https://ui.adsabs.harvard.edu/link_gateway/2013A&A...554A..26S/doi:10.1051/0004-6361/201321236>`_, `[arXiv:1301.4499] <https://arxiv.org/abs/1301.4499>`_
.. [2] Steininger et al., "NIFTy 3 - Numerical Information Field Theory - A Python framework for multicomponent signal inference on HPC clusters", 2017, accepted by Annalen der Physik; `[arXiv:1708.01073] <https://arxiv.org/abs/1708.01073>`_
Contents
........
.. toctree::
:maxdepth: 2
ift
Gallery <http://wwwmpa.mpa-garching.mpg.de/ift/nifty/gallery/>
installation
code
citations
Indices and tables
..................
* :any:`Module Index <mod/modules>`
* :ref:`search`
Package Documentation <mod/nifty5>
......@@ -7,18 +7,29 @@ distributions, the "apt" lines will need slight changes.
NIFTy5 and its mandatory dependencies can be installed via::
sudo apt-get install git libfftw3-dev python3 python3-pip python3-dev
sudo apt-get install git python3 python3-pip python3-dev
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_5
Plotting support is added via::
pip3 install --user matplotlib
FFTW support is added via:
sudo apt-get install libfftw3-dev
pip3 install --user pyfftw
To actually use FFTW in your Nifty calculations, you need to call
nifty5.fft.enable_fftw()
at the beginning of your code.
(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::
pip3 install --user matplotlib
Support for spherical harmonic transforms is added via::
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
......
......@@ -19,7 +19,6 @@ from .field import Field
from .multi_field import MultiField
from .operators.operator import Operator
from .operators.central_zero_padder import CentralZeroPadder
from .operators.diagonal_operator import DiagonalOperator
from .operators.distributors import DOFDistributor, PowerDistributor
from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter
......@@ -47,9 +46,10 @@ from .operators.outer_product_operator import OuterProduct
from .operators.simple_linear_operators import (
VdotOperator, ConjugationOperator, Realizer,
FieldAdapter, ducktape, GeometryRemover, NullOperator)
from .operators.value_inserter import ValueInserter
from .operators.energy_operators import (
EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
BernoulliEnergy, StandardHamiltonian, SampledKullbachLeiblerDivergence)
BernoulliEnergy, StandardHamiltonian, AveragedEnergy)
from .probing import probe_with_posterior_samples, probe_diagonal, \
StatCalculator
......
......@@ -16,10 +16,10 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
from ..utilities import NiftyMetaBase
from ..utilities import NiftyMeta
class Domain(NiftyMetaBase()):
class Domain(metaclass=NiftyMeta):
"""The abstract class repesenting a (structured or unstructured) domain.
"""
def __repr__(self):
......
......@@ -19,23 +19,63 @@ from .utilities import iscomplextype
import numpy as np
_use_fftw = True
_use_fftw = False
_fftw_prepped = False
_fft_extra_args = {}
if _use_fftw:
import pyfftw
from pyfftw.interfaces.numpy_fft import fftn, rfftn, ifftn
pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(1000.)
# Optional extra arguments for the FFT calls
# if exact reproducibility is needed,
# set "planner_effort" to "FFTW_ESTIMATE"
import os
nthreads = int(os.getenv("OMP_NUM_THREADS", "1"))
_fft_extra_args = dict(planner_effort='FFTW_ESTIMATE', threads=nthreads)
else:
from numpy.fft import fftn, rfftn, ifftn
_fft_extra_args = {}
def enable_fftw():
global _use_fftw
_use_fftw = True
def disable_fftw():
global _use_fftw
_use_fftw = False
def _init_pyfftw():
global _fft_extra_args, _fftw_prepped
if not _fftw_prepped:
import pyfftw
from pyfftw.interfaces.numpy_fft import fftn, rfftn, ifftn
pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(1000.)
# Optional extra arguments for the FFT calls
# if exact reproducibility is needed,
# set "planner_effort" to "FFTW_ESTIMATE"
import os
nthreads = int(os.getenv("OMP_NUM_THREADS", "1"))
_fft_extra_args = dict(planner_effort='FFTW_ESTIMATE',
threads=nthreads)
_fftw_prepped = True
def fftn(a, axes=None):
if _use_fftw:
from pyfftw.interfaces.numpy_fft import fftn
_init_pyfftw()
return fftn(a, axes=axes, **_fft_extra_args)
else:
return np.fft.fftn(a, axes=axes)
def rfftn(a, axes=None):
if _use_fftw:
from pyfftw.interfaces.numpy_fft import rfftn
_init_pyfftw()
return rfftn(a, axes=axes, **_fft_extra_args)
else:
return np.fft.rfftn(a, axes=axes)
def ifftn(a, axes=None):
if _use_fftw:
from pyfftw.interfaces.numpy_fft import ifftn
_init_pyfftw()
return ifftn(a, axes=axes, **_fft_extra_args)
else:
return np.fft.ifftn(a, axes=axes)
def hartley(a, axes=None):
......@@ -46,7 +86,7 @@ def hartley(a, axes=None):
if iscomplextype(a.dtype):
raise TypeError("Hartley transform requires real-valued arrays.")
tmp = rfftn(a, axes=axes, **_fft_extra_args)
tmp = rfftn(a, axes=axes)
def _fill_array(tmp, res, axes):
if axes is None:
......@@ -89,7 +129,7 @@ def my_fftn_r2c(a, axes=None):
if iscomplextype(a.dtype):
raise TypeError("Transform requires real-valued input arrays.")
tmp = rfftn(a, axes=axes, **_fft_extra_args)
tmp = rfftn(a, axes=axes)
def _fill_complex_array(tmp, res, axes):
if axes is None:
......@@ -123,4 +163,4 @@ def my_fftn_r2c(a, axes=None):
def my_fftn(a, axes=None):
return fftn(a, axes=axes, **_fft_extra_args)
return fftn(a, axes=axes)
......@@ -25,7 +25,7 @@ from .domain_tuple import DomainTuple
class Field(object):
"""The discrete representation of a continuous field over multiple spaces.
Stores data arrays and carries all the needed metainformation (i.e. the
Stores data arrays and carries all the needed meta-information (i.e. the
domain) for operators to be able to operate on them.
Parameters
......@@ -442,7 +442,7 @@ class Field(object):
----------
spaces : None, int or tuple of int
The summation is only carried out over the sub-domains in this
tuple. If None, it is carried out over all sub-domains. Default: None.
tuple. If None, it is carried out over all sub-domains.
Returns
-------
......@@ -463,7 +463,6 @@ class Field(object):
spaces : None, int or tuple of int
The summation is only carried out over the sub-domains in this
tuple. If None, it is carried out over all sub-domains.
Default: None.
Returns
-------
......@@ -547,7 +546,7 @@ class Field(object):
----------
spaces : None, int or tuple of int
The operation is only carried out over the sub-domains in this
tuple. If None, it is carried out over all sub-domains. Default: None.
tuple. If None, it is carried out over all sub-domains.
Returns
-------
......
......@@ -41,7 +41,7 @@ def make_adjust_variances(a,
Operator which gives the amplitude when evaluated at a position
xi : Operator
Operator which gives the excitation when evaluated at a position
postion : Field, MultiField
position : Field, MultiField
Position of the whole problem
samples : Field, MultiField
Residual samples of the whole problem
......
......@@ -25,15 +25,15 @@ from ..operators.simple_linear_operators import ducktape
def CorrelatedField(target, amplitude_operator, name='xi'):
'''Constructs an operator which turns a white Gaussian excitation field
"""Constructs an operator which turns a white Gaussian excitation field
into a correlated field.
This function returns an operator which implements:
ht @ (vol * A * xi),
where `ht` is a harmonic transform operator, `A` is the sqare root of the
prior covariance an `xi` is the excitation field.
where `ht` is a harmonic transform operator, `A` is the square root of the
prior covariance and `xi` is the excitation field.
Parameters
----------
......@@ -41,12 +41,12 @@ def CorrelatedField(target, amplitude_operator, name='xi'):
Target of the operator. Must contain exactly one space.
amplitude_operator: Operator
name : string
:class:`MultiField` key for xi-field.
:class:`MultiField` key for the xi-field.
Returns
-------
Correlated field : Operator
'''
"""
tgt = DomainTuple.make(target)
if len(tgt) > 1:
raise ValueError
......@@ -60,7 +60,7 @@ def CorrelatedField(target, amplitude_operator, name='xi'):
def MfCorrelatedField(target, amplitudes, name='xi'):
'''Constructs an operator which turns white Gaussian excitation fields
"""Constructs an operator which turns white Gaussian excitation fields
into a correlated field defined on a DomainTuple with two entries and two
separate correlation structures.
......@@ -79,7 +79,7 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
Returns
-------
Correlated field : Operator
'''
"""
tgt = DomainTuple.make(target)
if len(tgt) != 2:
raise ValueError
......
......@@ -34,7 +34,7 @@ def _float_or_listoffloat(inp):
return [float(x) for x in inp] if isinstance(inp, list) else float(inp)
def _make_dynamic_operator(domain,
def _make_dynamic_operator(target,
harmonic_padding,
sm_s0,
sm_x0,
......@@ -44,8 +44,10 @@ def _make_dynamic_operator(domain,
minimum_phase,
sigc=None,
quant=None):
if not isinstance(domain, RGSpace):
if not isinstance(target, RGSpace):
raise TypeError("RGSpace required")
if not target.harmonic:
raise TypeError("Target space must be harmonic")
if not (isinstance(harmonic_padding, int) or harmonic_padding is None
or all(isinstance(ii, int) for ii in harmonic_padding)):
raise TypeError
......@@ -62,7 +64,7 @@ def _make_dynamic_operator(domain,
if cone and (sigc is None or quant is None):
raise RuntimeError
dom = DomainTuple.make(domain)
dom = DomainTuple.make(target.get_default_codomain())
ops = {}
FFT = FFTOperator(dom)
Real = Realizer(dom)
......@@ -134,20 +136,32 @@ def _make_dynamic_operator(domain,
return m, ops
def dynamic_operator(domain,
def dynamic_operator(*,
target,
harmonic_padding,
sm_s0,
sm_x0,
key,
causal=True,
minimum_phase=False):
'''Constructs an operator encoding the Green's function of a linear
"""Constructs an operator encoding the Green's function of a linear
homogeneous dynamic system.
When evaluated, this operator returns the Green's function representation
in harmonic space. This result can be used as a convolution kernel to
construct solutions of the homogeneous stochastic differential equation
encoded in this operator. Note that if causal is True, the Green's function
is convolved with a step function in time, where the temporal axis is the
first axis of the space. In this case the resulting function only extends
up to half the length of the first axis of the space to avoid boundary
effects during convolution. If minimum_phase is true then the spectrum of
the Green's function is used to construct a corresponding minimum phase
filter.
Parameters
----------
domain : RGSpace
The position space in which the Green's function shall be constructed.
target : RGSpace
The harmonic space in which the Green's function shall be constructed.
harmonic_padding : None, int, list of int
Amount of central padding in harmonic space in pixels. If None the
field is not padded at all.
......@@ -159,13 +173,15 @@ def dynamic_operator(domain,
key for dynamics encoding parameter.
causal : boolean
Whether or not the Green's function shall be causal in time.
Default is True.
minimum_phase: boolean
Whether or not the Green's function shall be a minimum phase filter.
Default is False.
Returns
-------
Operator
The Operator encoding the dynamic Green's function in harmonic space.
The Operator encoding the dynamic Green's function in target space.
Dictionary of Operator
A collection of sub-chains of Operator which can be used for plotting
and evaluation.
......@@ -173,9 +189,9 @@ def dynamic_operator(domain,
Notes
-----
The first axis of the domain is interpreted the time axis.
'''
"""
dct = {
'domain': domain,
'target': target,
'harmonic_padding': harmonic_padding,
'sm_s0': sm_s0,
'sm_x0': sm_x0,
......@@ -187,7 +203,8 @@ def dynamic_operator(domain,
return _make_dynamic_operator(**dct)
def dynamic_lightcone_operator(domain,
def dynamic_lightcone_operator(*,
target,
harmonic_padding,
sm_s0,
sm_x0,
......@@ -199,11 +216,16 @@ def dynamic_lightcone_operator(domain,
minimum_phase=False):
'''Extends the functionality of :function: dynamic_operator to a Green's
function which is constrained to be within a light cone.
The resulting Green's function is constrained to be within a light cone.
This is achieved via convolution of the function with a light cone in
space-time. Thereby the first axis of the space is set to be the teporal
axis.
Parameters
----------
domain : RGSpace
The position space in which the Green's function shall be constructed.
target : RGSpace
The harmonic space in which the Green's function shall be constructed.
It needs to have at least two dimensions.
harmonic_padding : None, int, list of int
Amount of central padding in harmonic space in pixels. If None the
......@@ -222,8 +244,10 @@ def dynamic_lightcone_operator(domain,
Quantization of the light cone in pixels.
causal : boolean
Whether or not the Green's function shall be causal in time.
Default is True.
minimum_phase: boolean
Whether or not the Green's function shall be a minimum phase filter.
Default is False.
Returns
-------
......@@ -238,10 +262,10 @@ def dynamic_lightcone_operator(domain,
The first axis of the domain is interpreted the time axis.
'''
if len(domain.shape) < 2:
if len(target.shape) < 2:
raise ValueError("Space must be at least 2 dimensional!")
dct = {
'domain': domain,
'target': target,
'harmonic_padding': harmonic_padding,
'sm_s0': sm_s0,
'sm_x0': sm_x0,
......
......@@ -26,33 +26,33 @@ from ..sugar import makeOp
class InverseGammaOperator(Operator):
"""Operator which transforms a Gaussian into an inverse gamma distribution.
The pdf of the inverse gamma distribution is defined as follows:
.. math ::
\\frac{q^\\alpha}{\\Gamma(\\alpha)}x^{-\\alpha -1}\\exp \\left(-\\frac{q}{x}\\right)
That means that for large x the pdf falls off like :math:`x^(-\\alpha -1)`.
The mean of the pdf is at :math:`q / (\\alpha - 1)` if :math:`\\alpha > 1`.
The mode is :math:`q / (\\alpha + 1)`.
This transformation is implemented as a linear interpolation which maps a
Gaussian onto a inverse gamma distribution.
Parameters
----------
domain : Domain, tuple of Domain or DomainTuple
The domain on which the field shall be defined. This is at the same
time the domain and the target of the operator.
alpha : float
The alpha-parameter of the inverse-gamma distribution.
q : float
The q-parameter of the inverse-gamma distribution.
delta : float
distance between sampling points for linear interpolation.
"""
def __init__(self, domain, alpha, q, delta=0.001):
"""Operator which transforms a Gaussian into an inverse gamma distribution.
The pdf of the inverse gamma distribution is defined as follows:
.. math::
\frac {\beta ^{\alpha }}{\Gamma (\alpha )}}x^{-\alpha -1}\exp \left(-{\frac {\beta }{x}}\right)
That means that for large x the pdf falls off like x^(-alpha -1).
The mean of the pdf is at q / (alpha - 1) if alpha > 1.
The mode is q / (alpha + 1).
This transformation is implemented as a linear interpolation which
maps a Gaussian onto a inverse gamma distribution.
Parameters
----------
domain : Domain, tuple of Domain or DomainTuple
The domain on which the field shall be defined. This is at the same
time the domain and the target of the operator.
alpha : float
The alpha-parameter of the inverse-gamma distribution.
q : float
The q-parameter of the inverse-gamma distribution.
delta : float
distance between sampling points for linear interpolation.
"""
self._domain = self._target = DomainTuple.make(domain)
self._alpha, self._q, self._delta = float(alpha), float(q), float(delta)
self._xmin, self._xmax = -8.2, 8.2
......
......@@ -47,9 +47,9 @@ def _make_coords(domain, absolute=False):
return k_array
class LightConeDerivative(LinearOperator):
class _LightConeDerivative(LinearOperator):
def __init__(self, domain, target, derivatives):
super(LightConeDerivative, self).__init__()
super(_LightConeDerivative, self).__init__()
self._domain = domain
self._target = target
self._derivatives = derivatives
......@@ -67,7 +67,7 @@ class LightConeDerivative(LinearOperator):
return Field.from_global_data(self._tgt(mode), res)
def cone_arrays(c, domain, sigx, want_gradient):
def _cone_arrays(c, domain, sigx, want_gradient):
x = _make_coords(domain)
a = np.zeros(domain.shape, dtype=np.complex)
if want_gradient:
......@@ -96,6 +96,9 @@ def cone_arrays(c, domain, sigx, want_gradient):
class LightConeOperator(Operator):
'''
FIXME
'''
def __init__(self, domain, target, sigx):
self._domain = domain
self._target = target
......@@ -104,9 +107,9 @@ class LightConeOperator(Operator):
def apply(self, x):
islin = isinstance(x, Linearization)
val = x.val.to_global_data() if islin else x.to_global_data()
a, derivs = cone_arrays(val, self.target, self._sigx, islin)
a, derivs = _cone_arrays(val, self.target, self._sigx, islin)
res = Field.from_global_data(self.target, a)
if not islin:
return res
jac = LightConeDerivative(x.jac.target, self.target, derivs)(x.jac)
jac = _LightConeDerivative(x.jac.target, self.target, derivs)(x.jac)
return Linearization(res, jac, want_metric=x.want_metric)
......@@ -104,7 +104,7 @@ def apply_erf(wgt, dist, lo, mid, hi, sig, erf):
class LOSResponse(LinearOperator):
"""Line-of-sight response operator
This operator transforms from a single RGSpace to an unstructured domain
This operator transforms from a single RGSpace to an UnstructuredDomain
with as many entries as there were lines of sight passed to the
constructor. Adjoint application is also provided.
......
......@@ -33,7 +33,7 @@ def _ceps_kernel(k, a, k0):
def CepstrumOperator(target, a, k0):
'''Turns a white Gaussian random field into a smooth field on a LogRGSpace.
"""Turns a white Gaussian random field into a smooth field on a LogRGSpace.
Composed out of three operators:
......@@ -43,7 +43,7 @@ def CepstrumOperator(target, a, k0):
and ceps is the so-called cepstrum:
.. math::
\\mathrm{sqrt\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
\\mathrm{sqrt\\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
These operators are combined in this fashion in order to generate:
......@@ -69,10 +69,10 @@ def CepstrumOperator(target, a, k0):
regularization of the inverse laplace operator to be finite at zero.
Larger values for the cutoff results in a weaker constraining prior.
k0 : float, list of float
Strength of smothness prior in quefrency space (positive only) along
Strength of smoothness prior in quefrency space (positive only) along
each axis. If float then the strength is the same along each axis.
Larger values result in a weaker constraining prior.
'''
"""
a = float(a)
target = DomainTuple.make(target)
if a <= 0:
......@@ -112,7 +112,7 @@ def CepstrumOperator(target, a, k0):
return sym @ qht @ makeOp(cepstrum.sqrt())
def SLAmplitude(target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
'''Operator for parametrizing smooth amplitudes (square roots of power
spectra).
......
......@@ -187,6 +187,18 @@ class Linearization(object):
return self.__mul__(other)
def outer(self, other):
"""Computes the outer product of this Linearization with a Field or