Skip to content

GitLab

Commit 24147f38 authored by Torsten Ensslin's avatar Torsten Ensslin
Browse files

Merge branch 'NIFTy_5' into 'docstrings_torsten' 

# Conflicts:
#   nifty5/operators/energy_operators.py
parents 9a8c07ff 9c813309
Expand all Show whitespace changes
Inline Side-by-side
.gitignore
View file @ 24147f38
......@@ -11,6 +11,7 @@ setup.cfg
.svn/
*.csv
.pytest_cache/
*.png
# from https://github.com/github/gitignore/blob/master/Python.gitignore
......
demos/getting_started_3.py
View file @ 24147f38
......@@ -46,32 +46,32 @@ if __name__ == '__main__':
mode = 1
position_space = ift.RGSpace([128, 128])
harmonic_space = position_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
power_space = ift.PowerSpace(harmonic_space)
# Set up an amplitude operator for the field
# The parameters mean:
# 64 spectral bins
#
dct = {
'target': power_space,
'n_pix': 64, # 64 spectral bins
# Spectral smoothness (affects Gaussian process part)
# 3 = relatively high variance of spectral curbvature
# 0.4 = quefrency mode below which cepstrum flattens
#
'a': 3, # relatively high variance of spectral curbvature
'k0': .4, # quefrency mode below which cepstrum flattens
# Power-law part of spectrum:
# -5 = preferred power-law slope
# 0.5 = low variance of power-law slope
# 0.4 = y-intercept mean
# 0.3 = relatively high y-intercept variance
A = ift.AmplitudeOperator(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3)
'sm': -5, # preferred power-law slope
'sv': .5, # low variance of power-law slope
'im': .4, # y-intercept mean
'iv': .3 # relatively high y-intercept variance
}
A = ift.AmplitudeOperator(**dct)
# Build the operator for a correlated signal
harmonic_space = position_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
power_space = A.target[0]
power_distributor = ift.PowerDistributor(harmonic_space, power_space)
vol = ift.ScalingOperator(harmonic_space.scalar_dvol**(-0.5),
harmonic_space)
correlated_field = ht(
vol(power_distributor(A))*ift.ducktape(harmonic_space, None, 'xi'))
vol = harmonic_space.scalar_dvol**-0.5
xi = ift.ducktape(harmonic_space, None, 'xi')
correlated_field = ht(vol*power_distributor(A)*xi)
# Alternatively, one can use:
# correlated_field = ift.CorrelatedField(position_space, A)
......
docs/generate.sh
View file @ 24147f38
sphinx-apidoc -l -e -d 2 -o docs/source/mod nifty5
sphinx-apidoc -e -o docs/source/mod nifty5
sphinx-build -b html docs/source/ docs/build/
docs/source/code.rst
View file @ 24147f38
This diff is collapsed.
docs/source/conf.py
View file @ 24147f38
......@@ -10,6 +10,9 @@ master_doc = 'index'
napoleon_google_docstring = False
napoleon_numpy_docstring = True
napoleon_use_ivar = True
napoleon_use_admonition_for_notes = True
napoleon_use_admonition_for_examples = True
napoleon_use_admonition_for_references = True
project = u'NIFTy5'
copyright = u'2013-2019, Max-Planck-Society'
......
nifty5/__init__.py
View file @ 24147f38
......@@ -76,7 +76,8 @@ from .plot import Plot
from .library.amplitude_operator import AmplitudeOperator
from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse
from .library.dynamic_operator import dynamic_operator, dynamic_lightcone_operator
from .library.dynamic_operator import (dynamic_operator,
dynamic_lightcone_operator)
from .library.light_cone_operator import LightConeOperator
from .library.wiener_filter_curvature import WienerFilterCurvature
......
nifty5/domain_tuple.py
View file @ 24147f38
......@@ -103,7 +103,7 @@ class DomainTuple(object):
"""tuple of int: number of pixels along each axis
The shape of the array-like object required to store information
living on the DomainTuple.
defined on the DomainTuple.
"""
return self._shape
......@@ -112,7 +112,7 @@ class DomainTuple(object):
"""tuple of int: number of pixels along each axis on the local task
The shape of the array-like object required to store information
living on part of the domain which is stored on the local MPI task.
defined on part of the domain which is stored on the local MPI task.
"""
from .dobj import local_shape
return local_shape(self._shape)
......
nifty5/domains/domain.py
View file @ 24147f38
......@@ -80,7 +80,7 @@ class Domain(NiftyMetaBase()):
"""tuple of int: number of pixels along each axis
The shape of the array-like object required to store information
living on the domain.
defined on the domain.
"""
raise NotImplementedError
......@@ -89,7 +89,7 @@ class Domain(NiftyMetaBase()):
"""tuple of int: number of pixels along each axis on the local task
The shape of the array-like object required to store information
living on part of the domain which is stored on the local MPI task.
defined on part of the domain which is stored on the local MPI task.
"""
from ..dobj import local_shape
return local_shape(self.shape)
......
nifty5/domains/lm_space.py
View file @ 24147f38
......@@ -32,13 +32,13 @@ class LMSpace(StructuredDomain):
lmax : int
The maximum :math:`l` value of any spherical harmonic coefficient
:math:`a_{lm}` that is represented by this object.
Must be :math:`\ge 0`.
Must be :math:`\\ge 0`.
mmax : int, optional
The maximum :math:`m` value of any spherical harmonic coefficient
:math:`a_{lm}` that is represented by this object.
If not supplied, it is set to `lmax`.
Must be :math:`\ge 0` and :math:`\le` `lmax`.
Must be :math:`\\ge 0` and :math:`\\le` `lmax`.
"""
_needed_for_hash = ["_lmax", "_mmax"]
......
nifty5/domains/log_rg_space.py
View file @ 24147f38
......@@ -38,7 +38,7 @@ class LogRGSpace(StructuredDomain):
FIXME
harmonic : bool, optional
Whether the space represents a grid in position or harmonic space.
(default: False).
Default: False.
"""
_needed_for_hash = ['_shape', '_bindistances', '_t_0', '_harmonic']
......
nifty5/domains/power_space.py
View file @ 24147f38
......@@ -31,13 +31,13 @@ class PowerSpace(StructuredDomain):
----------
harmonic_partner : StructuredDomain
The harmonic domain of which this is the power space.
binbounds : None, or tuple of float (default: None)
if None:
There will be as many bins as there are distinct k-vector lengths
in the harmonic partner space.
The `binbounds` property of the PowerSpace will also be None.
binbounds : None, or tuple of float
By default (binbounds=None):
There are as many bins as there are distinct k-vector lengths in
the harmonic partner space.
The `binbounds` property of the PowerSpace will be None.
else:
the bin bounds requested for this PowerSpace. The array
The bin bounds requested for this PowerSpace. The array
must be sorted and strictly ascending. The first entry is the right
boundary of the first bin, and the last entry is the left boundary
of the last bin, i.e. thee will be `len(binbounds)+1` bins in
......
nifty5/domains/rg_space.py
View file @ 24147f38
......@@ -31,19 +31,21 @@ class RGSpace(StructuredDomain):
shape : int or tuple of int
Number of grid points or numbers of gridpoints along each axis.
distances : None or float or tuple of float, optional
Distance between two grid points along each axis
(default: None).
Distance between two grid points along each axis.
If distances is None:
- if harmonic==True, all distances will be set to 1
- if harmonic==False, the distance along each axis will be
By default (distances=None):
- If harmonic==True, all distances will be set to 1
- If harmonic==False, the distance along each axis will be
set to the inverse of the number of points along that axis.
harmonic : bool, optional
Whether the space represents a grid in position or harmonic space.
(default: False).
Default: False.
Notes
-----
Topologically, a n-dimensional RGSpace is a n-Torus, i.e. it has periodic
boundary conditions.
"""
_needed_for_hash = ["_distances", "_shape", "_harmonic"]
......
nifty5/field.py
View file @ 24147f38
......@@ -25,7 +25,7 @@ from .domain_tuple import DomainTuple
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
metainformation (i.e. the domain) for operators to be able to work on them.
......@@ -159,13 +159,13 @@ class Field(object):
Returns
-------
Field
Field living on `new_domain`, but with the same data as `self`.
Field defined on `new_domain`, but with the same data as `self`.
"""
return Field(DomainTuple.make(new_domain), self._val)
@staticmethod
def from_random(random_type, domain, dtype=np.float64, **kwargs):
""" Draws a random field with the given parameters.
"""Draws a random field with the given parameters.
Parameters
----------
......@@ -284,7 +284,7 @@ class Field(object):
return res
def weight(self, power=1, spaces=None):
""" Weights the pixels of `self` with their invidual pixel-volume.
"""Weights the pixels of `self` with their invidual pixel-volume.
Parameters
----------
......@@ -325,7 +325,7 @@ class Field(object):
return Field.from_local_data(self._domain, aout)
def outer(self, x):
""" Computes the outer product of 'self' with x.
"""Computes the outer product of 'self' with x.
Parameters
----------
......@@ -342,21 +342,21 @@ class Field(object):
return OuterProduct(self, x.domain)(x)
def vdot(self, x=None, spaces=None):
""" Computes the dot product of 'self' with x.
"""Computes the dot product of 'self' with x.
Parameters
----------
x : Field
x must be defined on the same domain as `self`.
spaces : None, int or tuple of int (default: None)
spaces : None, int or tuple of int
The dot product is only carried out over the sub-domains in this
tuple. If None, it is carried out over all sub-domains.
Default: None.
Returns
-------
float, complex, either scalar (for full dot products)
or Field (for partial dot products)
float, complex, either scalar (for full dot products) or Field (for partial dot products).
"""
if not isinstance(x, Field):
raise TypeError("The dot-partner must be an instance of " +
......@@ -375,12 +375,12 @@ class Field(object):
return (self.conjugate()*x).sum(spaces=spaces)
def norm(self, ord=2):
""" Computes the L2-norm of the field values.
"""Computes the L2-norm of the field values.
Parameters
----------
ord : int, default=2
accepted values: 1, 2, ..., np.inf
ord : int
Accepted values: 1, 2, ..., np.inf. Default: 2.
Returns
-------
......@@ -390,7 +390,7 @@ class Field(object):
return dobj.norm(self._val, ord)
def conjugate(self):
""" Returns the complex conjugate of the field.
"""Returns the complex conjugate of the field.
Returns
-------
......@@ -441,9 +441,9 @@ class Field(object):
Parameters
----------
spaces : None, int or tuple of int (default: None)
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.
tuple. If None, it is carried out over all sub-domains. Default: None.
Returns
-------
......@@ -461,9 +461,10 @@ class Field(object):
Parameters
----------
spaces : None, int or tuple of int (default: None)
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
-------
......@@ -484,9 +485,10 @@ class Field(object):
Parameters
----------
spaces : None, int or tuple of int (default: None)
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.
Returns
-------
......@@ -544,9 +546,9 @@ class Field(object):
Parameters
----------
spaces : None, int or tuple of int (default: None)
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.
tuple. If None, it is carried out over all sub-domains. Default: None.
Returns
-------
......@@ -566,9 +568,10 @@ class Field(object):
Parameters
----------
spaces : None, int or tuple of int (default: None)
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.
Returns
-------
......@@ -594,9 +597,10 @@ class Field(object):
Parameters
----------
spaces : None, int or tuple of int (default: None)
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.
Returns
-------
......
nifty5/library/adjust_variances.py
View file @ 24147f38
......@@ -30,7 +30,7 @@ def make_adjust_variances(a,
samples=[],
scaling=None,
ic_samp=None):
""" Creates a Hamiltonian for constant likelihood optimizations.
"""Creates a Hamiltonian for constant likelihood optimizations.
Constructs a Hamiltonian to solve constant likelihood optimizations of the
form phi = a * xi under the constraint that phi remains constant.
......
nifty5/library/amplitude_operator.py
View file @ 24147f38
......@@ -19,76 +19,51 @@ import numpy as np
from ..domains.power_space import PowerSpace
from ..field import Field
from ..sugar import makeOp, sqrt
from ..sugar import makeOp
def _ceps_kernel(dof_space, k, a, k0):
return a**2/(1 + (k/k0)**2)**2
def create_cepstrum_amplitude_field(domain, cepstrum):
"""Creates a ...
Writes the sum of all modes into the zero-mode.
Parameters
----------
domain: ???
???
cepstrum: Callable
???
"""
def _create_cepstrum_amplitude_field(domain, cepstrum):
dim = len(domain.shape)
shape = domain.shape
q_array = domain.get_k_array()
# Fill cepstrum field (all non-zero modes)
no_zero_modes = (slice(1, None),)*dim
ks = q_array[(slice(None),) + no_zero_modes]
# 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)
# Fill cepstrum field (zero-mode subspaces)
# Fill zero-mode subspaces
for i in range(dim):
# Prepare indices
fst_dims = (slice(None),)*i
sl = fst_dims + (slice(1, None),)
sl2 = fst_dims + (0,)
# Do summation
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)
def _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
def CepstrumOperator(domain, a, k0):
'''
Parameters
----------
ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_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
qht = QHTOperator(target=logk_space)
# FIXME a>0 k0>0
qht = QHTOperator(target=domain)
dof_space = qht.domain[0]
sym = SymmetrizingOperator(logk_space)
kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
res = sym(qht(makeOp(sqrt(cepstrum))))
if not zero_mode:
shp = res.target.shape
mask = np.ones(shp)
mask[(0,)*len(shp)] = 0.
mask = makeOp(Field.from_global_data(res.target, mask))
res = mask(res)
return res
def _SlopePowerSpectrum(logk_space, sm, sv, im, iv):
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())
def SlopeOperator(domain, sm, sv, im, iv):
'''
Parameters
----------
......@@ -97,54 +72,90 @@ def _SlopePowerSpectrum(logk_space, sm, sv, im, iv):
slope_variance (default=1)
im, iv : y-intercept_mean, y-intercept_std of power_slope
'''
'''
from ..operators.slope_operator import SlopeOperator
from ..operators.offset_operator import OffsetOperator
phi_mean = np.array([sm, im + sm*logk_space.t_0[0]])
# sv, iv>0
phi_mean = np.array([sm, im + sm*domain.t_0[0]])
phi_sig = np.array([sv, iv])
slope = SlopeOperator(logk_space)
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)))
def AmplitudeOperator(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi'], zero_mode=True):
''' Operator for parametrizing smooth power spectra.
def AmplitudeOperator(
target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
'''Operator for parametrizing smooth power spectra.
The general guideline for setting up generative models in IFT is to
transform the problem into the eigenbase of the prior and formulate the
generative model in this base. This is done here for the case of a power
spectrum which is smooth and has a linear component (both on
double-logarithmic scale).
This function assembles an :class:`Operator` which maps two a-priori white
Gaussian random fields to a smooth power spectrum which is composed out of
a linear and a smooth component.
On double-logarithmic scale, i.e. both x and y-axis on logarithmic scale,
the output of the generated operator is:
AmplitudeOperator = 0.5*(smooth_component + linear_component)
Computes a smooth power spectrum.
Output is defined on a PowerSpace.
This is then exponentiated and exponentially binned (via a :class:`ExpTransform`).
The prior on the linear component is parametrized by four real numbers,
being expected value and prior variance on the slope and the y-intercept
of the linear function.
The prior on the smooth component is parametrized by two real numbers: the
strength and the cutoff of the smoothness prior (see :class:`CepstrumOperator`).
Parameters
----------
Npixdof : int
#pix in dof_space
ceps_a : float
Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a, k0 = ceps_k0
ceps_k0 : float
Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a, k0 = ceps_k0
n_pix : int
Number of pixels of the space in which the .
target : PowerSpace
Target of the Operator.
a : float
Strength of smoothness prior (see :class:`CepstrumOperator`).
k0 : float
Cutoff of smothness prior in quefrency space (see :class:`CepstrumOperator`).
sm : float
slope_mean = expected exponent of power law (e.g. -4)
Expected exponent of power law (see :class:`SlopeOperator`).
sv : float
slope_variance (default=1)
Prior variance of exponent of power law (see :class:`SlopeOperator`).
im : float
y-intercept_mean
Expected y-intercept of power law (see :class:`SlopeOperator`).
iv : float
y-intercept_variance of power_slope
Prior variance of y-intercept of power law (see :class:`SlopeOperator`).
Returns
-------
Operator
Operator which is defined on the space of white excitations fields and
which returns on its target a power spectrum which consists out of a
smooth and a linear part.
'''
from ..operators.exp_transform import ExpTransform
from ..operators.scaling_operator import ScalingOperator
h_space = s_space.get_default_codomain()
et = ExpTransform(PowerSpace(h_space), Npixdof)
logk_space = et.domain[0]
if not (isinstance(n_pix, int) and isinstance(target, PowerSpace)):
raise TypeError
a, k0 = float(a), float(k0)
sm, sv, im, iv = float(sm), float(sv), float(im), float(iv)
et = ExpTransform(target, n_pix)
dom = et.domain[0]
dct = {'a': a, 'k0': k0}
smooth = CepstrumOperator(dom, **dct).ducktape(keys[0])
smooth = _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode)
smooth = smooth.ducktape(keys[0])
linear = _SlopePowerSpectrum(logk_space, sm, sv, im, iv)
linear = linear.ducktape(keys[1])
dct = {'sm': sm, 'sv': sv, 'im': im, 'iv': iv}
linear = SlopeOperator(dom, **dct).ducktape(keys[1])
fac = ScalingOperator(0.5, smooth.target)
return et((fac(smooth + linear)).exp())
return et((0.5*(smooth + linear)).exp())
nifty5/library/dynamic_operator.py
View file @ 24147f38
......@@ -18,8 +18,8 @@
import numpy as np
from ..domain_tuple import DomainTuple
from ..domains.unstructured_domain import UnstructuredDomain
from ..domains.rg_space import RGSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from ..operators.diagonal_operator import DiagonalOperator
from ..operators.field_zero_padder import FieldZeroPadder
......@@ -30,12 +30,39 @@ from ..sugar import makeOp
from .light_cone_operator import LightConeOperator, _field_from_function
def _make_dynamic_operator(domain, harmonic_padding, sm_s0, sm_x0, keys=['f', 'c'],
causal=True, cone=True, minimum_phase=False, sigc=3.,
quant=5.):
dom = DomainTuple.make(domain)
if not isinstance(dom[0], RGSpace):
def _float_or_listoffloat(inp):
return [float(x) for x in inp] if isinstance(inp, list) else float(inp)
def _make_dynamic_operator(domain,
harmonic_padding,
sm_s0,
sm_x0,
cone,
keys,
causal,