Commit 5689ba2d authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'docstrings_pa' into 'NIFTy_5'

Docstrings pa

See merge request ift/nifty-dev!169
parents 9c813309 55124bb3
......@@ -65,7 +65,7 @@ if __name__ == '__main__':
'im': .4, # y-intercept mean
'iv': .3 # relatively high y-intercept variance
}
A = ift.AmplitudeOperator(**dct)
A = ift.SLAmplitude(**dct)
# Build the operator for a correlated signal
power_distributor = ift.PowerDistributor(harmonic_space, power_space)
......
......@@ -212,7 +212,7 @@ specific inference problems. Currently these are:
.. currentmodule:: nifty5.library
- :class:`~amplitude_operator.AmplitudeOperator`, which returns a smooth power spectrum.
- :class:`~smooth_linear_amplitude.SLAmplitude`, which returns a smooth power spectrum.
- :class:`~inverse_gamma_operator.InverseGammaOperator`, which models point sources which are
distributed according to a inverse-gamma distribution.
- :class:`~correlated_fields.CorrelatedField`, which models a diffuse log-normal field. It takes an
......
......@@ -73,7 +73,7 @@ from .minimization.kl_energy import KL_Energy
from .sugar import *
from .plot import Plot
from .library.amplitude_operator import AmplitudeOperator
from .library.smooth_linear_amplitude import (SLAmplitude, CepstrumOperator)
from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse
from .library.dynamic_operator import (dynamic_operator,
......
......@@ -15,7 +15,7 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
# Data object module for NIFTy that uses simple numpy ndarrays.
# Data object module that uses simple numpy ndarrays.
import numpy as np
from numpy import absolute, clip, cos, cosh, empty, empty_like, exp, full, log
......
......@@ -21,7 +21,7 @@ from .structured_domain import StructuredDomain
class GLSpace(StructuredDomain):
"""NIFTy subclass for Gauss-Legendre pixelizations of the two-sphere.
"""Represents a 2-sphere with Gauss-Legendre pixelization.
Its harmonic partner domain is the
:class:`~nifty5.domains.lm_space.LMSpace`.
......
......@@ -21,7 +21,7 @@ from .structured_domain import StructuredDomain
class HPSpace(StructuredDomain):
"""NIFTy subclass for HEALPix discretizations of the two-sphere.
"""Represents 2-sphere with HEALPix discretization.
Its harmonic partner domain is the
:class:`~nifty5.domains.lm_space.LMSpace`.
......
......@@ -22,7 +22,7 @@ from .structured_domain import StructuredDomain
class LMSpace(StructuredDomain):
"""NIFTy subclass for sets of spherical harmonic coefficients.
"""Represents a set of spherical harmonic coefficients.
Its harmonic partner spaces are :class:`~nifty5.domains.hp_space.HPSpace`
and :class:`~nifty5.domains.gl_space.GLSpace`.
......
......@@ -16,30 +16,29 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
import numpy as np
from .. import dobj
from ..field import Field
from ..sugar import exp
from .structured_domain import StructuredDomain
class LogRGSpace(StructuredDomain):
"""NIFTy subclass for logarithmic Cartesian grids.
'''Represents a logarithmic Cartesian grid.
Parameters
----------
shape : int or tuple of int
Number of grid points or numbers of gridpoints along each axis.
bindistances : float or tuple of float
Distance between two grid points along each axis. These are
measured on logarithmic scale and are constant therfore.
Logarithmic distance between two grid points along each axis.
Equidistant spacing of bins on logarithmic scale is assumed.
t_0 : float or tuple of float
FIXME
Coordinate of pixel ndim*(1,).
harmonic : bool, optional
Whether the space represents a grid in position or harmonic space.
Default: False.
"""
'''
_needed_for_hash = ['_shape', '_bindistances', '_t_0', '_harmonic']
def __init__(self, shape, bindistances, t_0, harmonic=False):
......@@ -77,6 +76,7 @@ class LogRGSpace(StructuredDomain):
@property
def t_0(self):
"""np.ndarray : array of coordinates of pixel ndim*(1,)."""
return np.array(self._t_0)
def __repr__(self):
......@@ -84,16 +84,49 @@ class LogRGSpace(StructuredDomain):
self.shape, self.harmonic))
def get_default_codomain(self):
"""Returns a :class:`LogRGSpace` object representing the (position or
harmonic) partner domain of `self`, depending on `self.harmonic`. The
`bindistances` are transformed and `t_0` stays the same.
Returns
-------
LogRGSpace
The parter domain
"""
codomain_bindistances = 1./(self.bindistances*self.shape)
return LogRGSpace(self.shape, codomain_bindistances, self._t_0, True)
def get_k_length_array(self):
"""Generates array of distances to origin of the space.
Returns
-------
numpy.ndarray
Distances to origin of the space. If any index of the array is
zero then the distance is np.nan if self.harmonic True.
The dtype is float64, the shape is `self.shape`.
Raises
------
NotImplementedError
If `self.harmonic` is False.
"""
if not self.harmonic:
raise NotImplementedError
ks = self.get_k_array()
return Field.from_global_data(self, np.linalg.norm(ks, axis=0))
def get_k_array(self):
"""Generates coordinates of the space.
Returns
-------
numpy.ndarray
Coordinates of the space. If one index of the array is zero the
corresponding coordinate is -np.inf (np.nan) if self.harmonic is
False (True).
The dtype is float64 and shape: `(len(self.shape),) + self.shape`.
"""
ndim = len(self.shape)
k_array = np.zeros((ndim,) + self.shape)
dist = self.bindistances
......
......@@ -22,7 +22,7 @@ from .structured_domain import StructuredDomain
class PowerSpace(StructuredDomain):
"""NIFTy class for spaces of power spectra.
"""Represents non-equidistantly binned spaces for power spectra.
A power space is the result of a projection of a harmonic domain where
k-modes of equal length get mapped to one power index.
......
......@@ -24,7 +24,7 @@ from .structured_domain import StructuredDomain
class RGSpace(StructuredDomain):
"""NIFTy subclass for regular Cartesian grids.
"""Represents a regular Cartesian grid.
Parameters
----------
......
......@@ -21,7 +21,7 @@ from .domain import Domain
class StructuredDomain(Domain):
"""The abstract base class for all structured NIFTy domains.
"""The abstract base class for all structured domains.
An instance of a space contains information about the manifold's
geometry and enhances the functionality of Domain by methods that
......@@ -50,7 +50,7 @@ class StructuredDomain(Domain):
@property
def total_volume(self):
"""float : Total domain volume
"""float : Total domain volume.
Returns the sum over all the domain's pixel volumes.
"""
......@@ -63,7 +63,7 @@ class StructuredDomain(Domain):
raise NotImplementedError
def get_k_length_array(self):
"""k vector lengths, if applicable,
"""k vector lengths, if applicable.
Returns the length of the k vector for every pixel.
This method is only implemented for harmonic domains.
......
......@@ -23,18 +23,15 @@ from .domain_tuple import DomainTuple
class Field(object):
_scalar_dom = DomainTuple.scalar_domain()
"""The discrete representation of a continuous field over multiple spaces.
In NIFTy, Fields are used to store data arrays and carry all the needed
metainformation (i.e. the domain) for operators to be able to work on them.
Stores data arrays and carries all the needed metainformation (i.e. the
domain) for operators to be able to operate on them.
Parameters
----------
domain : DomainTuple
the domain of the new Field
val : data_object
This object's global shape must match the domain shape
After construction, the object will no longer be writeable!
......@@ -42,9 +39,11 @@ class Field(object):
Notes
-----
If possible, do not invoke the constructor directly, but use one of the
many convenience functions for Field construction!
many convenience functions for instantiation!
"""
_scalar_dom = DomainTuple.scalar_domain()
def __init__(self, domain, val):
if not isinstance(domain, DomainTuple):
raise TypeError("domain must be of type DomainTuple")
......@@ -337,7 +336,7 @@ class Field(object):
"""
if not isinstance(x, Field):
raise TypeError("The multiplier must be an instance of " +
"the NIFTy field class")
"the Field class")
from .operators.outer_product_operator import OuterProduct
return OuterProduct(self, x.domain)(x)
......@@ -360,7 +359,7 @@ class Field(object):
"""
if not isinstance(x, Field):
raise TypeError("The dot-partner must be an instance of " +
"the NIFTy field class")
"the Field class")
if x._domain != self._domain:
raise ValueError("Domain mismatch")
......
......@@ -15,64 +15,93 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
from ..domain_tuple import DomainTuple
from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import ducktape
from ..operators.scaling_operator import ScalingOperator
def CorrelatedField(s_space, amplitude_operator, name='xi'):
'''
Function for construction of correlated fields
def CorrelatedField(target, amplitude_operator, name='xi'):
'''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.
Parameters
----------
s_space : Domain
Field domain
target : Domain, DomainTuple or tuple of Domain
Target of the operator. Must contain exactly one space.
amplitude_operator: Operator
operator for correlation structure
name : string
MultiField component name
:class:`MultiField` key for xi-field.
Returns
-------
Correlated field : Operator
'''
h_space = s_space.get_default_codomain()
ht = HarmonicTransformOperator(h_space, s_space)
tgt = DomainTuple.make(target)
if len(tgt) > 1:
raise ValueError
h_space = tgt[0].get_default_codomain()
ht = HarmonicTransformOperator(h_space, tgt[0])
p_space = amplitude_operator.target[0]
power_distributor = PowerDistributor(h_space, p_space)
A = power_distributor(amplitude_operator)
vol = h_space.scalar_dvol
vol = ScalingOperator(vol**(-0.5), h_space)
return ht(vol(A)*ducktape(h_space, None, name))
vol = h_space.scalar_dvol**-0.5
return ht(vol*A*ducktape(h_space, None, name))
def MfCorrelatedField(s_space_spatial,
s_space_energy,
amplitude_operator_spatial,
amplitude_operator_energy,
name="xi"):
'''
Method for construction of correlated multi-frequency fields
def MfCorrelatedField(target, amplitudes, name='xi'):
'''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.
This operator may be used as a model for multi-frequency reconstructions
with a correlation structure in both spatial and energy direction.
Parameters
----------
target : Domain, DomainTuple or tuple of Domain
Target of the operator. Must contain exactly one space.
amplitudes: iterable of Operator
List of two amplitude operators.
name : string
:class:`MultiField` key for xi-field.
Returns
-------
Correlated field : Operator
'''
h_space_spatial = s_space_spatial.get_default_codomain()
h_space_energy = s_space_energy.get_default_codomain()
h_space = DomainTuple.make((h_space_spatial, h_space_energy))
ht1 = HarmonicTransformOperator(h_space, target=s_space_spatial, space=0)
ht2 = HarmonicTransformOperator(ht1.target, space=1)
ht = ht2(ht1)
tgt = DomainTuple.make(target)
if len(tgt) != 2:
raise ValueError
if len(amplitudes) != 2:
raise ValueError
p_space_spatial = amplitude_operator_spatial.target[0]
p_space_energy = amplitude_operator_energy.target[0]
hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
ht1 = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
ht2 = HarmonicTransformOperator(ht1.target, space=1)
ht = ht2 @ ht1
pd_spatial = PowerDistributor(h_space, p_space_spatial, 0)
pd_energy = PowerDistributor(pd_spatial.domain, p_space_energy, 1)
pd = pd_spatial(pd_energy)
psp = [aa.target[0] for aa in amplitudes]
pd0 = PowerDistributor(hsp, psp[0], 0)
pd1 = PowerDistributor(pd0.domain, psp[1], 1)
pd = pd0 @ pd1
dom_distr_spatial = ContractionOperator(pd.domain, 1).adjoint
dom_distr_energy = ContractionOperator(pd.domain, 0).adjoint
dd0 = ContractionOperator(pd.domain, 1).adjoint
dd1 = ContractionOperator(pd.domain, 0).adjoint
d = [dd0, dd1]
a_spatial = dom_distr_spatial(amplitude_operator_spatial)
a_energy = dom_distr_energy(amplitude_operator_energy)
a = a_spatial*a_energy
A = pd(a)
return ht(A*ducktape(h_space, None, name))
a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
a = reduce(lambda x, y: x*y, a)
A = pd @ a
vol = reduce(lambda x, y: x*y, [sp.scalar_dvol**-0.5 for sp in hsp])
return ht(vol*A*ducktape(hsp, None, name))
......@@ -147,7 +147,7 @@ def dynamic_operator(domain,
Parameters
----------
domain : RGSpace
The space under consideration.
The position 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.
......@@ -158,9 +158,9 @@ def dynamic_operator(domain,
key : String
key for dynamics encoding parameter.
causal : boolean
Whether or not the reconstructed dynamics should be causal in time.
Whether or not the Green's function shall be causal in time.
minimum_phase: boolean
Whether or not the reconstructed dynamics should be minimum phase.
Whether or not the Green's function shall be a minimum phase filter.
Returns
-------
......@@ -197,14 +197,14 @@ def dynamic_lightcone_operator(domain,
quant,
causal=True,
minimum_phase=False):
'''Constructs an operator encoding the Green's function of a linear
homogeneous dynamic system. The Green's function is constrained to be
within a light cone.
'''Extends the functionality of :function: dynamic_operator to a Green's
function which is constrained to be within a light cone.
Parameters
----------
domain : RGSpace
The space under consideration. Must have dim > 1.
The position 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
field is not padded at all.
......@@ -221,18 +221,17 @@ def dynamic_lightcone_operator(domain,
quant : float
Quantization of the light cone in pixels.
causal : boolean
Whether or not the reconstructed dynamics should be causal in time.
Whether or not the Green's function shall be causal in time.
minimum_phase: boolean
Whether or not the reconstructed dynamics should be minimum phase.
Whether or not the Green's function shall be a minimum phase filter.
Returns
-------
Operator
The Operator encoding the dynamic Green's function in harmonic space
when evaluated.
dict
A collection of sub-chains of Operator which can be used
for plotting and evaluation.
The Operator encoding the dynamic Green's function in harmonic space.
Dictionary of Operator
A collection of sub-chains of Operator which can be used for plotting
and evaluation.
Notes
-----
......
......@@ -17,88 +17,113 @@
import numpy as np
from ..domain_tuple import DomainTuple
from ..domains.power_space import PowerSpace
from ..field import Field
from ..operators.exp_transform import ExpTransform
from ..operators.offset_operator import OffsetOperator
from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
from ..sugar import makeOp
def _ceps_kernel(dof_space, k, a, k0):
return a**2/(1 + (k/k0)**2)**2
def _ceps_kernel(k, a, k0):
return (a/(1+np.sum((k.T/k0)**2, axis=-1).T))**2
def _create_cepstrum_amplitude_field(domain, cepstrum):
dim = len(domain.shape)
shape = domain.shape
q_array = domain.get_k_array()
def CepstrumOperator(target, a, k0):
'''Turns a white Gaussian random field into a smooth field on a LogRGSpace.
# Fill all non-zero modes
no_zero_modes = (slice(1, None), )*dim
ks = q_array[(slice(None), ) + no_zero_modes]
cepstrum_field = np.zeros(shape)
cepstrum_field[no_zero_modes] = cepstrum(ks)
Composed out of three operators:
# Fill zero-mode subspaces
for i in range(dim):
fst_dims = (slice(None), )*i
sl = fst_dims + (slice(1, None), )
sl2 = fst_dims + (0, )
cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
return Field.from_global_data(domain, cepstrum_field)
sym @ qht @ diag(sqrt_ceps),
where sym is a :class:`SymmetrizingOperator`, qht is a :class:`QHTOperator`
and ceps is the so-called cepstrum:
def CepstrumOperator(domain, a, k0):
'''
.. math::
C(k) = \\left(\\frac{a}{1+(k/k0)^2}\\right)^2
'''
from ..operators.qht_operator import QHTOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
# FIXME a>0 k0>0
qht = QHTOperator(target=domain)
dof_space = qht.domain[0]
sym = SymmetrizingOperator(domain)
kern = lambda k: _ceps_kernel(dof_space, k, a, k0)
cepstrum = _create_cepstrum_amplitude_field(dof_space, kern)
return sym @ qht @ makeOp(cepstrum.sqrt())
\\mathrm{sqrt\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
These operators are combined in this fashion in order to generate:
def SlopeOperator(domain, sm, sv, im, iv):
'''
Parameters
----------
- A field which is smooth, i.e. second derivatives are punished (note
that the sqrt-cepstrum is essentially proportional to 1/k**2).
sm, sv : slope_mean = expected exponent of power law (e.g. -4),
slope_variance (default=1)
- A field which is symmetric around the pixel in the middle of the space.
This is result of the :class:`SymmetrizingOperator` and needed in order
to decouple the degrees of freedom at the beginning and the end of the
amplitude whenever :class:`CepstrumOperator` is used as in
:class:`SLAmplitude`.
im, iv : y-intercept_mean, y-intercept_std of power_slope
The prior on the zero mode (or zero subspaces for more than one dimensions)
is the integral of the prior over all other modes along the corresponding
axis.
Parameters
----------
target : LogRGSpace
Target domain of the operator, needs to be non-harmonic.
a : float
Cutoff of smoothness prior (positive only). Controls the
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
each axis. If float then the strength is the same along each axis.
Larger values result in a weaker constraining prior.
'''
from ..operators.slope_operator import SlopeOperator
from ..operators.offset_operator import OffsetOperator
# sv, iv>0
a = float(a)
target = DomainTuple.make(target)
if a <= 0:
raise ValueError
if len(target) > 1 or target[0].harmonic:
raise TypeError
if isinstance(k0, (float, int)):
k0 = np.array([k0]*len(target.shape))
else:
k0 = np.array(k0)
if len(k0) != len(target.shape):
raise ValueError
if np.any(np.array(k0) <= 0):
raise ValueError
qht = QHTOperator(target)
dom = qht.domain[0]
sym = SymmetrizingOperator(target)
# Compute cepstrum field
dim = len(dom.shape)
shape = dom.shape
q_array = dom.get_k_array()
# Fill all non-zero modes
no_zero_modes = (slice(1, None),)*dim
ks = q_array[(slice(None),) + no_zero_modes]
cepstrum_field = np.zeros(shape)
cepstrum_field[no_zero_modes] = _ceps_kernel(ks, a, k0)
# Fill zero-mode subspaces
for i in range(dim):
fst_dims = (slice(None),)*i
sl = fst_dims + (slice(1, None),)
sl2 = fst_dims + (0,)
cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
cepstrum = Field.from_global_data(dom, cepstrum_field)
phi_mean = np.array([sm, im + sm*domain.t_0[0]])
phi_sig = np.array([sv, iv])
slope = SlopeOperator(domain)
phi_mean = Field.from_global_data(slope.domain, phi_mean)
phi_sig = Field.from_global_data(slope.domain, phi_sig)
return slope(OffsetOperator(phi_mean)(makeOp(phi_sig)))
return sym @ qht @ makeOp(cepstrum.sqrt())
def AmplitudeOperator(
target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
'''Operator for parametrizing smooth power spectra.
def SLAmplitude(target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
'''Operator for parametrizing smooth amplitudes (square roots of power