Commit b81eb80e authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'field_docu' into docstring_operators

parents 8bb5b571 223f71e3
Pipeline #12476 passed with stage
in 9 minutes and 4 seconds
......@@ -13,7 +13,7 @@ before_script:
- apt-get update
- >
apt-get install -y build-essential python python-pip python-dev git
autoconf gsl-bin libgsl-dev wget python-numpy cython
autoconf gsl-bin libgsl-dev wget python-numpy
- pip install --upgrade -r ci/requirements_base.txt
- chmod +x ci/*.sh
# This Makefile implements common tasks needed by developers
# A list of implemented rules can be obtained by the command "make help"
help :
echo " Implemented targets:"
echo " build build pypmc for python2 and python3"
echo " buildX build pypmc for pythonX only where X is one of {2,3}"
echo " build-sdist build pypmc from the dist directory (python 2 and 3)"
echo " build-sdistX build pypmc from the dist directory (pythonX, X in {2,3})"
echo " check use nosetests to test pypmc with python 2.7 and 3"
echo " checkX use nosetests to test pypmc with python 2.7 or 3,"
echo " where X is one of {2,3}"
echo " check-fast use nosetests to run only quick tests of pypmc"
echo " using nosetests-2.7 and nosetests3"
echo " check-sdist use nosetests-2.7 and nosetests3 to test the distribution"
echo " generated by 'make sdist'"
echo " check-sdistX use nosetests-2.7 or nosetests3 to test the distribution"
echo " generated by 'make sdist', where X is one of {2,3}"
echo " clean delete compiled and temporary files"
echo " coverage produce and show a code coverage report"
echo " Note: Cython modules cannot be analyzed"
echo " distcheck runs 'check', check-sdist', 'run-examples' and"
echo " opens a browser with the built documentation"
echo " doc build the html documentation using sphinx"
echo " doc-pdf build the pdf documentation using sphinx"
echo " help show this message"
echo " run-examples run all examples using python 2 and 3"
echo " sdist make a source distribution"
echo " show-todos show todo marks in the source code"
.PHONY : clean
#remove build doc
rm -rf ./doc/_build
#remove .pyc files created by python 2.7
rm -f ./*.pyc
find -P . -name '*.pyc' -delete
#remove .pyc files crated by python 3
rm -rf ./__pycache__
find -P . -name __pycache__ -delete
#remove build folder in root directory
rm -rf ./build
#remove cythonized C source and object files
find -P . -name '*.c' -delete
#remove variational binaries only if command line argument specified
find -P . -name '*.so' -delete
#remove backup files
find -P . -name '*~' -delete
#remove files created by coverage
rm -f .coverage
rm -rf coverage
# remove egg info
rm -rf pypmc.egg-info
# remove downloaded seutptools
rm -f
# remove dist/
rm -rf dist
.PHONY : build
build : build2
.PHONY : build2
build2 :
python2 build_ext --inplace
check : check2
.PHONY : check2
check2 : build2
@ # run tests
nosetests-2.7 --processes=-1 --process-timeout=60
# run tests in parallel
mpirun -n 2 nosetests-2.7
.PHONY : check-fast
check-fast : build
nosetests-2.7 -a '!slow' --processes=-1 --process-timeout=60
nosetests3 -a '!slow' --processes=-1 --process-timeout=60
.PHONY : .build-system-default
.build-system-default :
python build_ext --inplace
.PHONY : doc
doc : .build-system-default
cd doc && make html
.PHONY : doc-pdf
doc-pdf : .build-system-default
cd doc; make latexpdf
.PHONY : run-examples
run-examples : build
cd examples ; \
for file in $$(ls) ; do \
echo running $${file} with python2 && \
python2 $${file} && \
echo running $${file} with python3 && \
python3 $${file} && \
# execute with mpirun if mpi4py appears in the file \
if grep -Fq 'mpi4py' $${file} ; then \
echo "$${file}" is mpi parallelized && \
echo running $${file} in parallel with python2 && \
mpirun -n 2 python2 $${file} && \
echo running $${file} in parallel with python3 && \
mpirun -n 2 python3 $${file} ; \
fi \
; \
.PHONY : sdist
sdist :
python sdist
.PHONY : build-sdist
build-sdist : build-sdist2 build-sdist3
./dist/pypmc*/NUL : sdist
cd dist && tar xaf *.tar.gz && cd *
.PHONY : build-sdist2
build-sdist2 : ./dist/pypmc*/NUL
cd dist/pypmc* && python2 build
.PHONY : build-sdist3
build-sdist3 : ./dist/pypmc*/NUL
cd dist/pypmc* && python3 build
.PHONY : check-sdist
check-sdist : check-sdist2 check-sdist3
.PHONY : check-sdist2
check-sdist2 : build-sdist2
cd dist/*/build/lib*2.7 && \
nosetests-2.7 --processes=-1 --process-timeout=60 && \
mpirun -n 2 nosetests-2.7
.PHONY : check-sdist3
check-sdist3 : build-sdist3
cd dist/*/build/lib*3.* && \
nosetests3 --processes=-1 --process-timeout=60 && \
mpirun -n 2 nosetests3
.PHONY : distcheck
distcheck : check check-sdist doc
@ # execute "run-examples" after all other recipes makes are done
make run-examples
xdg-open link_to_documentation
.PHONY : show-todos
grep_cmd = ack-grep -i --no-html --no-cc [^"au""sphinx.ext."]todo
show-todos :
@ # suppress errors here
@ # note that no todo found is considered as error
$(grep_cmd) doc ; \
$(grep_cmd) pypmc ; \
$(grep_cmd) examples ; echo \
.PHONY : coverage
coverage : .build-system-default
rm -rf coverage
nosetests --with-coverage --cover-package=nifty --cover-html --cover-html-dir=coverage
xdg-open coverage/index.html
Metadata-Version: 1.0
Name: ift_nifty
Version: 1.0.6
Summary: Numerical Information Field Theory
Author: Theo Steininger
License: GPLv3
Description: UNKNOWN
Platform: UNKNOWN
......@@ -15,7 +15,7 @@ Summary
a versatile library designed to enable the development of signal
inference algorithms that operate regardless of the underlying spatial
grid and its resolution. Its object-oriented framework is written in
Python, although it accesses libraries written in Cython, C++, and C for
Python, although it accesses libraries written in C++ and C for
NIFTY offers a toolkit that abstracts discretized representations of
......@@ -71,7 +71,6 @@ Installation
- [Python]( (v2.7.x)
- [NumPy](
- [Cython](
### Download
......@@ -95,7 +94,7 @@ Starting with a fresh Ubuntu installation move to a folder like
- Using pip install numpy etc...:
sudo pip install numpy cython
sudo pip install numpy
- Install pyHealpix:
......@@ -147,10 +146,9 @@ MacPorts, missing ones need to be installed manually. It may also be
mentioned that one should only use one package manager, as multiple ones
may cause trouble.
- Install basic packages numpy and cython:
- Install numpy:
sudo port install py27-numpy
sudo port install py27-cython
- Install pyHealpix:
......@@ -11,16 +11,36 @@ rank = comm.rank
if __name__ == "__main__":
distribution_strategy = 'not'
#Setting up physical constants
#total length of Interval or Volume the field lives on, e.g. in meters
L = 2.
#typical distance over which the field is correlated (in same unit as L)
correlation_length = 0.1
#variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
field_variance = 2.
#smoothing length that response (in same unit as L)
response_sigma = 0.1
#defining resolution (pixels per dimension)
N_pixels = 512
#Setting up derived constants
k_0 = 1./correlation_length
#note that field_variance**2 = a*k_0/4. for this analytic form of power
a = field_variance**2/k_0*4.
pow_spec = (lambda k: a / (1 + k/k_0) ** 4)
pixel_width = L/N_pixels
# Setting up the geometry
s_space = RGSpace([512, 512])
s_space = RGSpace([N_pixels, N_pixels], distances = pixel_width)
fft = FFTOperator(s_space)
h_space =[0]
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Creating the mock data
pow_spec = (lambda k: 42 / (k + 1) ** 3)
S = create_power_operator(h_space, power_spectrum=pow_spec,
......@@ -30,7 +50,7 @@ if __name__ == "__main__":
sh = sp.power_synthesize(real_signal=True)
ss = fft.inverse_times(sh)
R = SmoothingOperator(s_space, sigma=0.1)
R = SmoothingOperator(s_space, sigma=response_sigma)
signal_to_noise = 1
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
......@@ -91,7 +91,7 @@ if __name__ == "__main__":
x1 = HPSpace(64)
k1 = HPLMTransformation.get_codomain(x1)
p1 = PowerSpace(harmonic_domain=k1, log=False)
p1 = PowerSpace(harmonic_partner=k1, logarithmic=False)
# creating Power Operator with given spectrum
spec = (lambda k: a_s / (1 + (k / k_0) ** 2) ** 2)
......@@ -107,7 +107,7 @@ if __name__ == "__main__":
# drawing a random field
sk = p_field.power_synthesize(decompose_power=True, mean=0.)
sk = p_field.power_synthesize(real_power=True, mean=0.)
s = Fft_op.adjoint_times(sk)
signal_to_noise = 1
......@@ -36,6 +36,75 @@ from nifty.random import Random
class Field(Loggable, Versionable, object):
""" Combines D2O-objects with meta-information NIFTY needs for computations.
In NIFTY, Fields are used to store dataarrays and carry all the needed
metainformation (i.e. the domain) for operators to be able to work on them.
In addition Field has functions to analyse the power spectrum of it's
content or create a random field of it.
domain : DomainObject
One of the space types NIFTY supports. RGSpace, GLSpace, HPSpace,
LMSpace or PowerSpace. It might also be a FieldArray, which is
an unstructured domain.
val : scalar, numpy.ndarray, distributed_data_object, Field
The values the array should contain after init. A scalar input will
fill the whole array with this scalar. If an array is provided the
array's dimensions must match the domain's.
dtype : type
A numpy.type. Most common are int, float and complex.
distribution_strategy: optional[{'fftw', 'equal', 'not', 'freeform'}]
Specifies which distributor will be created and used.
'fftw' uses the distribution strategy of pyfftw,
'equal' tries to distribute the data as uniform as possible
'not' does not distribute the data at all
'freeform' distribute the data according to the given local data/shape
copy: boolean
val : distributed_data_object
domain : DomainObject
See Parameters.
domain_axes : tuple of tuples
Enumerates the axes of the Field
dtype : type
Contains the datatype stored in the Field.
distribution_strategy : string
Name of the used distribution_strategy.
Raised if
*the given domain contains something that is not a DomainObject
*val is an array that has a different dimension than the domain
>>> a = Field(RGSpace([4,5]),val=2)
>>> a.val
array([[2, 2, 2, 2, 2],
[2, 2, 2, 2, 2],
[2, 2, 2, 2, 2],
[2, 2, 2, 2, 2]])
>>> a.dtype
See Also
# ---Initialization methods---
def __init__(self, domain=None, val=None, dtype=None,
......@@ -57,6 +126,21 @@ class Field(Loggable, Versionable, object):
self.set_val(new_val=val, copy=copy)
def _parse_domain(self, domain, val=None):
""" Returns a tuple of DomainObjects for nomenclature unification.
domain : all supported NIFTY spaces
The domain over which the Field lives.
val : a NIFTY Field instance
Can be used to make Field infere it's domain by adopting val's
out : tuple
The output object. A tuple with one or multiple DomainObjects.
if domain is None:
if isinstance(val, Field):
domain = val.domain
......@@ -75,6 +159,29 @@ class Field(Loggable, Versionable, object):
return domain
def _get_axes_tuple(self, things_with_shape, start=0):
""" Enumerates all axes of the domain.
This function is used in the greater context of the 'spaces' keyword.
things_with_shape : indexable list of objects with .shape property
Normal input is a domain/ tuple of domains.
start : int
Sets the integer number for the first axis
out : tuple
Incremental numeration of all axes.
The 'spaces' keyword is used in operators in order to carry out
operations only on a certain subspace if the domain of the Field is
a product space.
i = start
axes_list = []
for thing in things_with_shape:
......@@ -86,6 +193,20 @@ class Field(Loggable, Versionable, object):
return tuple(axes_list)
def _infer_dtype(self, dtype, val):
""" Inferes the datatype of the Field
dtype : type
Can be None
val : list of arrays
If the dtype is None, Fields tries to infere the datatype from the
values given to it at initialization.
out : np.dtype
if dtype is None:
dtype = val.dtype
......@@ -121,6 +242,35 @@ class Field(Loggable, Versionable, object):
def from_random(cls, random_type, domain=None, dtype=None,
distribution_strategy=None, **kwargs):
""" Draws a random field with the given parameters.
cls : class
random_type : String
'pm1', 'normal', 'uniform' are the supported arguments for this
domain : DomainObject
The domain of the output random field
dtype : type
The datatype of the output random field
distribution_strategy : all supported distribution strategies
The distribution strategy of the output random field
out : Field
The output object.
See Also
_parse_random_arguments, power_synthesise
# create a initially empty field
f = cls(domain=domain, dtype=dtype,
......@@ -143,7 +293,6 @@ class Field(Loggable, Versionable, object):
def _parse_random_arguments(random_type, f, **kwargs):
if random_type == "pm1":
random_arguments = {}
......@@ -167,8 +316,58 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods---
def power_analyze(self, spaces=None, logarithmic=False, nbin=None, binbounds=None,
def power_analyze(self, spaces=None, log=False, nbin=None, binbounds=None,
""" Computes the powerspectrum of the Field
Creates a PowerSpace with the given attributes and computes the
power spectrum as a field over this PowerSpace.
It's important to note that this can only be done if the subspace to
be analyzed is in harmonic space.
spaces : int, *optional*
The subspace which you want to have the powerspectrum of.
{default : None}
if spaces==None : Tries to synthesize for the whole domain
log : boolean, *optional*
True if the output PowerSpace should have log binning.
{default : False}
nbin : int, None, *optional*
The number of bins the resulting PowerSpace shall have.
{default : None}
if nbin==None : maximum number of bins is used
binbounds : array-like, None, *optional*
Inner bounds of the bins, if specifield
{default : None}
if binbounds==None : bins are inferred. Overwrites nbins and log
real_signal : boolean, *optional*
Whether the analysed signal-space Field is real or complex.
For a real field a complex power spectrum comes out.
For a compex field all power is put in a real power spectrum.
{default : True}
Raised if
*len(spaces) is either 0 or >1
*len(domain) is not 1 with spaces=None
*the analyzed space is not harmonic
out : Field
The output object. It's domain is a PowerSpace and it contains
the power spectrum of 'self's field.
See Also
power_synthesize, PowerSpace
# check if all spaces in `self.domain` are either harmonic or
# power_space instances
for sp in self.domain:
......@@ -210,18 +409,18 @@ class Field(Loggable, Versionable, object):
harmonic_partner = self.domain[space_index]
power_domain = PowerSpace(harmonic_partner=harmonic_partner,
harmonic_domain = self.domain[space_index]
power_domain = PowerSpace(harmonic_domain=harmonic_domain,
logarithmic=logarithmic, nbin=nbin, binbounds=binbounds)
log=log, nbin=nbin, binbounds=binbounds)
# extract pindex and rho from power_domain
pindex = power_domain.pindex
rho = power_domain.rho
if decompose_power:
if real_signal:
hermitian_part, anti_hermitian_part = \
......@@ -245,7 +444,7 @@ class Field(Loggable, Versionable, object):
result_domain = list(self.domain)
result_domain[space_index] = power_domain
if decompose_power:
if real_signal:
result_dtype = np.complex
result_dtype = np.float
......@@ -303,9 +502,55 @@ class Field(Loggable, Versionable, object):
return result_obj
def power_synthesize(self, spaces=None, real_power=True,
real_signal=False, mean=None, std=None):
def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
mean=None, std=None):
"""Yields a random field in harmonic space with this power spectrum.
This method draws a Gaussian random field in the harmic partner domain.
The drawn field has this field as its power spectrum.
For this the domain must be a PowerSpace.
spaces : {tuple, int, None} *optional*
Specifies the subspace in which the power will be synthesized in
case of a product space.
{default : None}
if spaces==None : Tries to synthesize for the whole domain
real_power : boolean *optional*
Determines whether the power spectrum is real or complex
{default : True}
real_signal : boolean *optional*
True will result in a purely real signal-space field.
This means that the created field is symmetric wrt. the origin
after complex conjugation.
{default : True}
mean : {float, None} *optional*
The mean of the noise field the powerspectrum will be multiplied on.
{default : None}
if mean==None : mean will be set to 0
std : float *optional*
The standard deviation of the noise field the powerspectrum will be
multiplied on.
{default : None}
if std==None : std will be set to 1