Commit abb05cee authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge NIFTy_5

parents 09d93746 0aa213dd
...@@ -5,7 +5,7 @@ RUN apt-get update && apt-get install -y \ ...@@ -5,7 +5,7 @@ RUN apt-get update && apt-get install -y \
git \ git \
# Packages needed for NIFTy # Packages needed for NIFTy
libfftw3-dev \ libfftw3-dev \
python3 python3-pip python3-dev python3-future python3-scipy cython3 \ python3 python3-pip python3-dev python3-scipy \
# Documentation build dependencies # Documentation build dependencies
python3-sphinx python3-sphinx-rtd-theme \ python3-sphinx python3-sphinx-rtd-theme \
# Testing dependencies # Testing dependencies
......
...@@ -65,7 +65,7 @@ if __name__ == '__main__': ...@@ -65,7 +65,7 @@ if __name__ == '__main__':
'im': .4, # y-intercept mean 'im': .4, # y-intercept mean
'iv': .3 # relatively high y-intercept variance 'iv': .3 # relatively high y-intercept variance
} }
A = ift.AmplitudeOperator(**dct) A = ift.SLAmplitude(**dct)
# Build the operator for a correlated signal # Build the operator for a correlated signal
power_distributor = ift.PowerDistributor(harmonic_space, power_space) power_distributor = ift.PowerDistributor(harmonic_space, power_space)
......
...@@ -212,7 +212,7 @@ specific inference problems. Currently these are: ...@@ -212,7 +212,7 @@ specific inference problems. Currently these are:
.. currentmodule:: nifty5.library .. 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 - :class:`~inverse_gamma_operator.InverseGammaOperator`, which models point sources which are
distributed according to a inverse-gamma distribution. distributed according to a inverse-gamma distribution.
- :class:`~correlated_fields.CorrelatedField`, which models a diffuse log-normal field. It takes an - :class:`~correlated_fields.CorrelatedField`, which models a diffuse log-normal field. It takes an
......
...@@ -19,7 +19,6 @@ from .field import Field ...@@ -19,7 +19,6 @@ from .field import Field
from .multi_field import MultiField from .multi_field import MultiField
from .operators.operator import Operator from .operators.operator import Operator
from .operators.central_zero_padder import CentralZeroPadder
from .operators.diagonal_operator import DiagonalOperator from .operators.diagonal_operator import DiagonalOperator
from .operators.distributors import DOFDistributor, PowerDistributor from .operators.distributors import DOFDistributor, PowerDistributor
from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter
...@@ -73,7 +72,7 @@ from .minimization.kl_energy import KL_Energy ...@@ -73,7 +72,7 @@ from .minimization.kl_energy import KL_Energy
from .sugar import * from .sugar import *
from .plot import Plot 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.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse from .library.los_response import LOSResponse
from .library.dynamic_operator import (dynamic_operator, from .library.dynamic_operator import (dynamic_operator,
......
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # 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 import numpy as np
from numpy import absolute, clip, cos, cosh, empty, empty_like, exp, full, log from numpy import absolute, clip, cos, cosh, empty, empty_like, exp, full, log
......
...@@ -21,7 +21,7 @@ from .structured_domain import StructuredDomain ...@@ -21,7 +21,7 @@ from .structured_domain import StructuredDomain
class GLSpace(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 Its harmonic partner domain is the
:class:`~nifty5.domains.lm_space.LMSpace`. :class:`~nifty5.domains.lm_space.LMSpace`.
......
...@@ -21,7 +21,7 @@ from .structured_domain import StructuredDomain ...@@ -21,7 +21,7 @@ from .structured_domain import StructuredDomain
class HPSpace(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 Its harmonic partner domain is the
:class:`~nifty5.domains.lm_space.LMSpace`. :class:`~nifty5.domains.lm_space.LMSpace`.
......
...@@ -22,7 +22,7 @@ from .structured_domain import StructuredDomain ...@@ -22,7 +22,7 @@ from .structured_domain import StructuredDomain
class LMSpace(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` Its harmonic partner spaces are :class:`~nifty5.domains.hp_space.HPSpace`
and :class:`~nifty5.domains.gl_space.GLSpace`. and :class:`~nifty5.domains.gl_space.GLSpace`.
......
...@@ -16,30 +16,29 @@ ...@@ -16,30 +16,29 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce from functools import reduce
import numpy as np import numpy as np
from .. import dobj
from ..field import Field from ..field import Field
from ..sugar import exp
from .structured_domain import StructuredDomain from .structured_domain import StructuredDomain
class LogRGSpace(StructuredDomain): class LogRGSpace(StructuredDomain):
"""NIFTy subclass for logarithmic Cartesian grids. '''Represents a logarithmic Cartesian grid.
Parameters Parameters
---------- ----------
shape : int or tuple of int shape : int or tuple of int
Number of grid points or numbers of gridpoints along each axis. Number of grid points or numbers of gridpoints along each axis.
bindistances : float or tuple of float bindistances : float or tuple of float
Distance between two grid points along each axis. These are Logarithmic distance between two grid points along each axis.
measured on logarithmic scale and are constant therfore. Equidistant spacing of bins on logarithmic scale is assumed.
t_0 : float or tuple of float t_0 : float or tuple of float
FIXME Coordinate of pixel ndim*(1,).
harmonic : bool, optional harmonic : bool, optional
Whether the space represents a grid in position or harmonic space. Whether the space represents a grid in position or harmonic space.
Default: False. Default: False.
""" '''
_needed_for_hash = ['_shape', '_bindistances', '_t_0', '_harmonic'] _needed_for_hash = ['_shape', '_bindistances', '_t_0', '_harmonic']
def __init__(self, shape, bindistances, t_0, harmonic=False): def __init__(self, shape, bindistances, t_0, harmonic=False):
...@@ -77,6 +76,7 @@ class LogRGSpace(StructuredDomain): ...@@ -77,6 +76,7 @@ class LogRGSpace(StructuredDomain):
@property @property
def t_0(self): def t_0(self):
"""np.ndarray : array of coordinates of pixel ndim*(1,)."""
return np.array(self._t_0) return np.array(self._t_0)
def __repr__(self): def __repr__(self):
...@@ -84,16 +84,49 @@ class LogRGSpace(StructuredDomain): ...@@ -84,16 +84,49 @@ class LogRGSpace(StructuredDomain):
self.shape, self.harmonic)) self.shape, self.harmonic))
def get_default_codomain(self): 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) codomain_bindistances = 1./(self.bindistances*self.shape)
return LogRGSpace(self.shape, codomain_bindistances, self._t_0, True) return LogRGSpace(self.shape, codomain_bindistances, self._t_0, True)
def get_k_length_array(self): 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: if not self.harmonic:
raise NotImplementedError raise NotImplementedError
ks = self.get_k_array() ks = self.get_k_array()
return Field.from_global_data(self, np.linalg.norm(ks, axis=0)) return Field.from_global_data(self, np.linalg.norm(ks, axis=0))
def get_k_array(self): 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) ndim = len(self.shape)
k_array = np.zeros((ndim,) + self.shape) k_array = np.zeros((ndim,) + self.shape)
dist = self.bindistances dist = self.bindistances
......
...@@ -22,7 +22,7 @@ from .structured_domain import StructuredDomain ...@@ -22,7 +22,7 @@ from .structured_domain import StructuredDomain
class PowerSpace(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 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. k-modes of equal length get mapped to one power index.
......
...@@ -24,7 +24,7 @@ from .structured_domain import StructuredDomain ...@@ -24,7 +24,7 @@ from .structured_domain import StructuredDomain
class RGSpace(StructuredDomain): class RGSpace(StructuredDomain):
"""NIFTy subclass for regular Cartesian grids. """Represents a regular Cartesian grid.
Parameters Parameters
---------- ----------
......
...@@ -21,7 +21,7 @@ from .domain import Domain ...@@ -21,7 +21,7 @@ from .domain import Domain
class StructuredDomain(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 An instance of a space contains information about the manifold's
geometry and enhances the functionality of Domain by methods that geometry and enhances the functionality of Domain by methods that
...@@ -50,7 +50,7 @@ class StructuredDomain(Domain): ...@@ -50,7 +50,7 @@ class StructuredDomain(Domain):
@property @property
def total_volume(self): def total_volume(self):
"""float : Total domain volume """float : Total domain volume.
Returns the sum over all the domain's pixel volumes. Returns the sum over all the domain's pixel volumes.
""" """
...@@ -63,7 +63,7 @@ class StructuredDomain(Domain): ...@@ -63,7 +63,7 @@ class StructuredDomain(Domain):
raise NotImplementedError raise NotImplementedError
def get_k_length_array(self): 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. Returns the length of the k vector for every pixel.
This method is only implemented for harmonic domains. This method is only implemented for harmonic domains.
......
...@@ -23,18 +23,15 @@ from .domain_tuple import DomainTuple ...@@ -23,18 +23,15 @@ from .domain_tuple import DomainTuple
class Field(object): class Field(object):
_scalar_dom = DomainTuple.scalar_domain()
"""The discrete representation of a continuous field over multiple spaces. """The discrete representation of a continuous field over multiple spaces.
In NIFTy, Fields are used to store data arrays and carry all the needed Stores data arrays and carries all the needed metainformation (i.e. the
metainformation (i.e. the domain) for operators to be able to work on them. domain) for operators to be able to operate on them.
Parameters Parameters
---------- ----------
domain : DomainTuple domain : DomainTuple
the domain of the new Field the domain of the new Field
val : data_object val : data_object
This object's global shape must match the domain shape This object's global shape must match the domain shape
After construction, the object will no longer be writeable! After construction, the object will no longer be writeable!
...@@ -42,9 +39,11 @@ class Field(object): ...@@ -42,9 +39,11 @@ class Field(object):
Notes Notes
----- -----
If possible, do not invoke the constructor directly, but use one of the 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): def __init__(self, domain, val):
if not isinstance(domain, DomainTuple): if not isinstance(domain, DomainTuple):
raise TypeError("domain must be of type DomainTuple") raise TypeError("domain must be of type DomainTuple")
...@@ -337,7 +336,7 @@ class Field(object): ...@@ -337,7 +336,7 @@ class Field(object):
""" """
if not isinstance(x, Field): if not isinstance(x, Field):
raise TypeError("The multiplier must be an instance of " + raise TypeError("The multiplier must be an instance of " +
"the NIFTy field class") "the Field class")
from .operators.outer_product_operator import OuterProduct from .operators.outer_product_operator import OuterProduct
return OuterProduct(self, x.domain)(x) return OuterProduct(self, x.domain)(x)
...@@ -360,7 +359,7 @@ class Field(object): ...@@ -360,7 +359,7 @@ class Field(object):
""" """
if not isinstance(x, Field): if not isinstance(x, Field):
raise TypeError("The dot-partner must be an instance of " + raise TypeError("The dot-partner must be an instance of " +
"the NIFTy field class") "the Field class")
if x._domain != self._domain: if x._domain != self._domain:
raise ValueError("Domain mismatch") raise ValueError("Domain mismatch")
......
...@@ -15,64 +15,93 @@ ...@@ -15,64 +15,93 @@
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
from ..domain_tuple import DomainTuple from ..domain_tuple import DomainTuple
from ..operators.contraction_operator import ContractionOperator from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import ducktape from ..operators.simple_linear_operators import ducktape
from ..operators.scaling_operator import ScalingOperator
def CorrelatedField(s_space, amplitude_operator, name='xi'): def CorrelatedField(target, amplitude_operator, name='xi'):
''' '''Constructs an operator which turns a white Gaussian excitation field
Function for construction of correlated fields 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 Parameters
---------- ----------
s_space : Domain target : Domain, DomainTuple or tuple of Domain
Field domain Target of the operator. Must contain exactly one space.
amplitude_operator: Operator amplitude_operator: Operator
operator for correlation structure
name : string name : string
MultiField component name :class:`MultiField` key for xi-field.
Returns
-------
Correlated field : Operator
''' '''
h_space = s_space.get_default_codomain() tgt = DomainTuple.make(target)
ht = HarmonicTransformOperator(h_space, s_space) 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] p_space = amplitude_operator.target[0]
power_distributor = PowerDistributor(h_space, p_space) power_distributor = PowerDistributor(h_space, p_space)
A = power_distributor(amplitude_operator) A = power_distributor(amplitude_operator)
vol = h_space.scalar_dvol vol = h_space.scalar_dvol**-0.5
vol = ScalingOperator(vol**(-0.5), h_space) return ht(vol*A*ducktape(h_space, None, name))
return ht(vol(A)*ducktape(h_space, None, name))
def MfCorrelatedField(s_space_spatial, def MfCorrelatedField(target, amplitudes, name='xi'):
s_space_energy, '''Constructs an operator which turns white Gaussian excitation fields
amplitude_operator_spatial, into a correlated field defined on a DomainTuple with two entries and two
amplitude_operator_energy, separate correlation structures.
name="xi"):
''' This operator may be used as a model for multi-frequency reconstructions
Method for construction of correlated multi-frequency fields 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() tgt = DomainTuple.make(target)
h_space_energy = s_space_energy.get_default_codomain() if len(tgt) != 2:
h_space = DomainTuple.make((h_space_spatial, h_space_energy)) raise ValueError
ht1 = HarmonicTransformOperator(h_space, target=s_space_spatial, space=0) if len(amplitudes) != 2:
ht2 = HarmonicTransformOperator(ht1.target, space=1) raise ValueError
ht = ht2(ht1)
p_space_spatial = amplitude_operator_spatial.target[0] hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
p_space_energy = amplitude_operator_energy.target[0] 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) psp = [aa.target[0] for aa in amplitudes]
pd_energy = PowerDistributor(pd_spatial.domain, p_space_energy, 1) pd0 = PowerDistributor(hsp, psp[0], 0)
pd = pd_spatial(pd_energy) pd1 = PowerDistributor(pd0.domain, psp[1], 1)
pd = pd0 @ pd1
dom_distr_spatial = ContractionOperator(pd.domain, 1).adjoint dd0 = ContractionOperator(pd.domain, 1).adjoint
dom_distr_energy = ContractionOperator(pd.domain, 0).adjoint dd1 = ContractionOperator(pd.domain, 0).adjoint
d = [dd0, dd1]
a_spatial = dom_distr_spatial(amplitude_operator_spatial) a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
a_energy = dom_distr_energy(amplitude_operator_energy) a = reduce(lambda x, y: x*y, a)
a = a_spatial*a_energy A = pd @ a
A = pd(a) vol = reduce(lambda x, y: x*y, [sp.scalar_dvol**-0.5 for sp in hsp])
return ht(A*ducktape(h_space, None, name)) return ht(vol*A*ducktape(hsp, None, name))
...@@ -147,7 +147,7 @@ def dynamic_operator(domain, ...@@ -147,7 +147,7 @@ def dynamic_operator(domain,
Parameters Parameters
---------- ----------
domain : RGSpace 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 harmonic_padding : None, int, list of int
Amount of central padding in harmonic space in pixels. If None the Amount of central padding in harmonic space in pixels. If None the
field is not padded at all. field is not padded at all.
...@@ -158,9 +158,9 @@ def dynamic_operator(domain, ...@@ -158,9 +158,9 @@ def dynamic_operator(domain,
key : String key : String
key for dynamics encoding parameter. key for dynamics encoding parameter.
causal : boolean 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 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 Returns
------- -------
...@@ -197,14 +197,14 @@ def dynamic_lightcone_operator(domain, ...@@ -197,14 +197,14 @@ def dynamic_lightcone_operator(domain,
quant, quant,
causal=True, causal=True,
minimum_phase=False): minimum_phase=False):
'''Constructs an operator encoding the Green's function of a linear '''Extends the functionality of :function: dynamic_operator to a Green's
homogeneous dynamic system. The Green's function is constrained to be function which is constrained to be within a light cone.
within a light cone.
Parameters Parameters
---------- ----------
domain : RGSpace 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 harmonic_padding : None, int, list of int
Amount of central padding in harmonic space in pixels. If None the Amount of central padding in harmonic space in pixels. If None the
field is not padded at all. field is not padded at all.
...@@ -221,18 +221,17 @@ def dynamic_lightcone_operator(domain, ...@@ -221,18 +221,17 @@ def dynamic_lightcone_operator(domain,
quant : float quant : float
Quantization of the light cone in pixels. Quantization of the light cone in pixels.
causal : boolean 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 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.