Commit 102c1aad authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'master' into fftw_for_the_masses

# Conflicts:
#	test/test_spaces/test_power_space.py
parents 693d7194 05d1d092
Pipeline #13333 failed with stage
in 5 minutes and 19 seconds
.. currentmodule:: nifty
The ``PropagatorOperator`` class -- ...
.......................................
.. autoclass:: PropagatorOperator
:show-inheritance:
:members:
.. currentmodule:: nifty
The ``ResponseOperator`` class -- A possible response implementation
....................................................................
.. autoclass:: ResponseOperator
:show-inheritance:
:members:
.. currentmodule:: nifty
The ``RGSpace`` class -- Regular Cartesian grids
................................................
.. autoclass:: RGSpace
:show-inheritance:
:members:
.. currentmodule:: nifty
The ``RGRGTransformation`` class -- A transformation routine
............................................................
.. autoclass:: RGRGTransformation
:show-inheritance:
:members:
.. currentmodule:: nifty
The ``SmoothingOperator`` class -- Smoothing fields
...................................................
.. autoclass:: SmoothingOperator
:show-inheritance:
:members:
Spaces
======
The :py:class:`Space` classes of NIFTY represent geometrical spaces approximated by grids in the computer environment. Each subclass of the base class corresponds to a specific grid type and replaces some of the inherited methods with its own methods that are unique to the respective grid. This framework ensures an abstract handling of spaces independent of the underlying geometrical grid and the grid's resolution.
Each instance of a :py:class:`Space` needs to capture all structural and dimensional specifics of the grid and all computationally relevant quantities such as the data type of associated field values. These parameters are stored as properties of an instance of the class at its initialization, and they do not need to be accessed explicitly by the user thereafter. This prevents the writing of grid or resolution dependent code.
Spatial symmetries of a system can be exploited by corresponding coordinate transformations. Often, transformations from one basis to its harmonic counterpart can greatly reduce the computational complexity of algorithms. The harmonic basis is defined by the eigenbasis of the Laplace operator; e.g., for a flat position space it is the Fourier basis. This conjugation of bases is implemented in NIFTY by distinguishing conjugate space classes, which can be obtained by the instance method *get_codomain* (and checked for by *check_codomain*). Moreover, transformations between conjugate spaces are performed automatically if required.
Space classes
-------------
Next to the generic :py:class:`Space` class, NIFTY has implementations of five subclasses, representing specific geometrical spaces and their discretizations.
.. toctree::
:maxdepth: 1
rg_space
hp_space
gl_space
lm_space
power_space
.. currentmodule:: nifty
The ``Space`` class -- The base Space object
--------------------------------------------
.. autoclass:: Space
:show-inheritance:
:members:
Transformations
---------------
NIFTY provides transformations
.. toctree::
:maxdepth: 1
rgrg_transformation
gllm_transformation
hplm_transformation
lmgl_transformation
lmhp_transformation
......@@ -269,7 +269,7 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods---
def power_analyze(self, spaces=None, logarithmic=False, nbin=None,
binbounds=None, decompose_power=True):
binbounds=None, keep_phase_information=False):
""" Computes the square root power spectrum for a subspace of `self`.
Creates a PowerSpace for the space addressed by `spaces` with the given
......@@ -283,7 +283,6 @@ class Field(Loggable, Versionable, object):
spaces : int *optional*
The subspace for which the powerspectrum shall be computed
(default : None).
if spaces==None : Tries to synthesize for the whole domain
logarithmic : boolean *optional*
True if the output PowerSpace should use logarithmic binning.
{default : False}
......@@ -294,11 +293,16 @@ class Field(Loggable, Versionable, object):
binbounds : array-like *optional*
Inner bounds of the bins (default : None).
if binbounds==None : bins are inferred. Overwrites nbins and log
decompose_power : boolean, *optional*
Whether the analysed signal-space Field is intrinsically real or
complex and if the power spectrum shall therefore be computed
for the real and the imaginary part of the Field separately
(default : True).
keep_phase_information : boolean, *optional*
If False, return a real-valued result containing the power spectrum
of the input Field.
If True, return a complex-valued result whose real component
contains the power spectrum computed from the real part of the
input Field, and whose imaginary component contains the power
spectrum computed from the imaginary part of the input Field.
The absolute value of this result should be identical to the output
of power_analyze with keep_phase_information=False.
(default : False).
Raise
-----
......@@ -331,23 +335,47 @@ class Field(Loggable, Versionable, object):
# check if the `spaces` input is valid
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None:
if len(self.domain) == 1:
spaces = (0,)
else:
raise ValueError(
"Field has multiple spaces as domain "
"but `spaces` is None.")
spaces = range(len(self.domain))
if len(spaces) == 0:
raise ValueError(
"No space for analysis specified.")
elif len(spaces) > 1:
raise ValueError(
"Conversion of only one space at a time is allowed.")
space_index = spaces[0]
if keep_phase_information:
parts_val = self._hermitian_decomposition(
domain=self.domain,
val=self.val,
spaces=spaces,
domain_axes=self.domain_axes,
preserve_gaussian_variance=False)
parts = [self.copy_empty().set_val(part_val, copy=False)
for part_val in parts_val]
else:
parts = [self]
parts = [abs(part)**2 for part in parts]
for space_index in spaces:
parts = [self._single_power_analyze(
work_field=part,
space_index=space_index,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds)
for part in parts]
if keep_phase_information:
result_field = parts[0] + 1j*parts[1]
else:
result_field = parts[0]
return result_field
@classmethod
def _single_power_analyze(cls, work_field, space_index, logarithmic, nbin,
binbounds):
if not self.domain[space_index].harmonic:
if not work_field.domain[space_index].harmonic:
raise ValueError(
"The analyzed space must be harmonic.")
......@@ -358,10 +386,10 @@ class Field(Loggable, Versionable, object):
# If it was complex, all the power is put into a real power spectrum.
distribution_strategy = \
self.val.get_axes_local_distribution_strategy(
self.domain_axes[space_index])
work_field.val.get_axes_local_distribution_strategy(
work_field.domain_axes[space_index])
harmonic_domain = self.domain[space_index]
harmonic_domain = work_field.domain[space_index]
power_domain = PowerSpace(harmonic_partner=harmonic_domain,
distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
......@@ -371,35 +399,18 @@ class Field(Loggable, Versionable, object):
pindex = power_domain.pindex
rho = power_domain.rho
if decompose_power:
hermitian_part, anti_hermitian_part = \
harmonic_domain.hermitian_decomposition(
self.val,
axes=self.domain_axes[space_index])
[hermitian_power, anti_hermitian_power] = \
[self._calculate_power_spectrum(
x=part,
pindex=pindex,
rho=rho,
axes=self.domain_axes[space_index])
for part in [hermitian_part, anti_hermitian_part]]
power_spectrum = hermitian_power + 1j * anti_hermitian_power
else:
power_spectrum = self._calculate_power_spectrum(
x=self.val,
pindex=pindex,
rho=rho,
axes=self.domain_axes[space_index])
power_spectrum = cls._calculate_power_spectrum(
field_val=work_field.val,
pindex=pindex,
rho=rho,
axes=work_field.domain_axes[space_index])
# create the result field and put power_spectrum into it
result_domain = list(self.domain)
result_domain = list(work_field.domain)
result_domain[space_index] = power_domain
result_dtype = power_spectrum.dtype
result_field = self.copy_empty(
result_field = work_field.copy_empty(
domain=result_domain,
dtype=result_dtype,
distribution_strategy=power_spectrum.distribution_strategy)
......@@ -407,17 +418,16 @@ class Field(Loggable, Versionable, object):
return result_field
def _calculate_power_spectrum(self, x, pindex, rho, axes=None):
fieldabs = abs(x)
fieldabs **= 2
@classmethod
def _calculate_power_spectrum(cls, field_val, pindex, rho, axes=None):
if axes is not None:
pindex = self._shape_up_pindex(
pindex=pindex,
target_shape=x.shape,
target_strategy=x.distribution_strategy,
axes=axes)
power_spectrum = pindex.bincount(weights=fieldabs,
pindex = cls._shape_up_pindex(
pindex=pindex,
target_shape=field_val.shape,
target_strategy=field_val.distribution_strategy,
axes=axes)
power_spectrum = pindex.bincount(weights=field_val,
axis=axes)
if axes is not None:
new_rho_shape = [1, ] * len(power_spectrum.shape)
......@@ -425,10 +435,10 @@ class Field(Loggable, Versionable, object):
rho = rho.reshape(new_rho_shape)
power_spectrum /= rho
power_spectrum **= 0.5
return power_spectrum
def _shape_up_pindex(self, pindex, target_shape, target_strategy, axes):
@staticmethod
def _shape_up_pindex(pindex, target_shape, target_strategy, axes):
if pindex.distribution_strategy not in \
DISTRIBUTION_STRATEGIES['global']:
raise ValueError("pindex's distribution strategy must be "
......@@ -543,6 +553,8 @@ class Field(Loggable, Versionable, object):
# components
spec = self.val.get_full_data()
spec = np.sqrt(spec)
for power_space_index in spaces:
spec = self._spec_to_rescaler(spec, result_list, power_space_index)
local_rescaler = spec
......@@ -561,10 +573,11 @@ class Field(Loggable, Versionable, object):
if real_signal:
result_val_list = [self._hermitian_decomposition(
result_domain,
result_val,
spaces,
result_list[0].domain_axes)[0]
result_domain,
result_val,
spaces,
result_list[0].domain_axes,
preserve_gaussian_variance=True)[0]
for result_val in result_val_list]
# store the result into the fields
......@@ -579,12 +592,13 @@ class Field(Loggable, Versionable, object):
return result
@staticmethod
def _hermitian_decomposition(domain, val, spaces, domain_axes):
def _hermitian_decomposition(domain, val, spaces, domain_axes,
preserve_gaussian_variance=False):
# hermitianize for the first space
(h, a) = domain[spaces[0]].hermitian_decomposition(
val,
domain_axes[spaces[0]],
preserve_gaussian_variance=True)
val,
domain_axes[spaces[0]],
preserve_gaussian_variance=preserve_gaussian_variance)
# hermitianize all remaining spaces using the iterative formula
for space in xrange(1, len(spaces)):
(hh, ha) = domain[space].hermitian_decomposition(
......@@ -1142,8 +1156,8 @@ class Field(Loggable, Versionable, object):
"""
return_field = self.copy_empty()
new_val = abs(self.get_val(copy=False))
return_field = self.copy_empty(dtype=new_val.dtype)
return_field.set_val(new_val, copy=False)
return return_field
......
......@@ -23,29 +23,29 @@ from .line_search import LineSearch
class LineSearchStrongWolfe(LineSearch):
"""Class for finding a step size that satisfies the strong Wolfe conditions.
Algorithm contains two stages. It begins whit a trial step length and it
Algorithm contains two stages. It begins whit a trial step length and it
keeps increasing the it until it finds an acceptable step length or an
interval. If it does not satisfy the Wolfe conditions it performs the Zoom
algorithm (second stage). By interpolating it decreases the size of the
interval until an acceptable step length is found.
interval. If it does not satisfy the Wolfe conditions it performs the Zoom
algorithm (second stage). By interpolating it decreases the size of the
interval until an acceptable step length is found.
Parameters
----------
c1 : float
c1 : float
Parameter for Armijo condition rule. (Default: 1e-4)
c2 : float
Parameter for curvature condition rule. (Default: 0.9)
max_step_size : float
Maximum step allowed in to be made in the descent direction.
Maximum step allowed in to be made in the descent direction.
(Default: 50)
max_iterations : integer
Maximum number of iterations performed by the line search algorithm.
(Default: 10)
max_zoom_iterations : integer
Maximum number of iterations performed by the zoom algorithm.
Maximum number of iterations performed by the zoom algorithm.
(Default: 10)
Attributes
----------
c1 : float
......@@ -53,19 +53,18 @@ class LineSearchStrongWolfe(LineSearch):
c2 : float
Parameter for curvature condition rule.
max_step_size : float
Maximum step allowed in to be made in the descent direction.
Maximum step allowed in to be made in the descent direction.
max_iterations : integer
Maximum number of iterations performed by the line search algorithm.
max_zoom_iterations : integer
Maximum number of iterations performed by the zoom algorithm.
"""
def __init__(self, c1=1e-4, c2=0.9,
max_step_size=50, max_iterations=10,
max_zoom_iterations=10):
super(LineSearchStrongWolfe, self).__init__()
self.c1 = np.float(c1)
......@@ -76,11 +75,11 @@ class LineSearchStrongWolfe(LineSearch):
def perform_line_search(self, energy, pk, f_k_minus_1=None):
"""Performs the first stage of the algorithm.
It starts with a trial step size and it keeps increasing it until it
satisfy the strong Wolf conditions. It also performs the descent and
It starts with a trial step size and it keeps increasing it until it
satisfy the strong Wolf conditions. It also performs the descent and
returns the optimal step length and the new enrgy.
Parameters
----------
energy : Energy object
......@@ -89,9 +88,9 @@ class LineSearchStrongWolfe(LineSearch):
pk : Field
Unit vector pointing into the search direction.
f_k_minus_1 : float
Value of the fuction (which is being minimized) at the k-1
Value of the fuction (which is being minimized) at the k-1
iteration of the line search procedure. (Default: None)
Returns
-------
alpha_star : float
......@@ -100,9 +99,9 @@ class LineSearchStrongWolfe(LineSearch):
Value of the energy after the performed descent.
energy_star : Energy object
The new Energy object on the new position.
"""
"""
self._set_line_energy(energy, pk, f_k_minus_1=f_k_minus_1)
c1 = self.c1
c2 = self.c2
......@@ -195,11 +194,11 @@ class LineSearchStrongWolfe(LineSearch):
def _zoom(self, alpha_lo, alpha_hi, phi_0, phiprime_0,
phi_lo, phiprime_lo, phi_hi, c1, c2):
"""Performs the second stage of the line search algorithm.
If the first stage was not successful then the Zoom algorithm tries to
find a suitable step length by using bisection, quadratic, cubic
If the first stage was not successful then the Zoom algorithm tries to
find a suitable step length by using bisection, quadratic, cubic
interpolation.
Parameters
----------
alpha_lo : float
......@@ -207,24 +206,24 @@ class LineSearchStrongWolfe(LineSearch):
alph_hi : float
The upper boundary for the step length interval.
phi_0 : float
Value of the energy at the starting point of the line search
Value of the energy at the starting point of the line search
algorithm.
phiprime_0 : Field
Gradient at the starting point of the line search algorithm.
phi_lo : float
Value of the energy if we perform a step of length alpha_lo in
Value of the energy if we perform a step of length alpha_lo in
descent direction.
phiprime_lo : Field
Gradient at the nwe position if we perform a step of length
Gradient at the nwe position if we perform a step of length
alpha_lo in descent direction.
phi_hi : float
Value of the energy if we perform a step of length alpha_hi in
Value of the energy if we perform a step of length alpha_hi in
descent direction.
c1 : float
Parameter for Armijo condition rule.
c2 : float
Parameter for curvature condition rule.
Returns
-------
alpha_star : float
......@@ -233,7 +232,7 @@ class LineSearchStrongWolfe(LineSearch):
Value of the energy after the performed descent.
energy_star : Energy object
The new Energy object on the new position.
"""
max_iterations = self.max_zoom_iterations
# define the cubic and quadratic interpolant checks
......@@ -306,12 +305,13 @@ class LineSearchStrongWolfe(LineSearch):
def _cubicmin(self, a, fa, fpa, b, fb, c, fc):
"""Estimating the minimum with cubic interpolation.
Finds the minimizer for a cubic polynomial that goes through the
points ( a,f(a) ), ( b,f(b) ), and ( c,f(c) ) with derivative at point a of fpa.
points ( a,f(a) ), ( b,f(b) ), and ( c,f(c) ) with derivative at point
a of fpa.
f(x) = A *(x-a)^3 + B*(x-a)^2 + C*(x-a) + D
If no minimizer can be found return None
Parameters
----------
a : float
......@@ -328,12 +328,12 @@ class LineSearchStrongWolfe(LineSearch):
Selected point.
fc : float
Value of polynomial at point c.
Returns
-------
xmin : float
Position of the approximated minimum.
"""
with np.errstate(divide='raise', over='raise', invalid='raise'):
......@@ -341,12 +341,12 @@ class LineSearchStrongWolfe(LineSearch):
C = fpa
db = b - a
dc = c - a
denom = (db * dc) ** 2 * (db - dc)
denom = db * db * dc * dc * (db - dc)
d1 = np.empty((2, 2))
d1[0, 0] = dc ** 2
d1[0, 1] = -db ** 2
d1[1, 0] = -dc ** 3
d1[1, 1] = db ** 3
d1[0, 0] = dc * dc
d1[0, 1] = -(db*db)
d1[1, 0] = -(dc*dc*dc)
d1[1, 1] = db*db*db
[A, B] = np.dot(d1, np.asarray([fb - fa - C * db,
fc - fa - C * dc]).flatten())
A /= denom
......@@ -361,11 +361,11 @@ class LineSearchStrongWolfe(LineSearch):
def _quadmin(self, a, fa, fpa, b, fb):
"""Estimating the minimum with quadratic interpolation.
Finds the minimizer for a quadratic polynomial that goes through
the points ( a,f(a) ), ( b,f(b) ) with derivative at point a of fpa.
f(x) = B*(x-a)^2 + C*(x-a) + D
Parameters
----------
a : float
......@@ -378,11 +378,11 @@ class LineSearchStrongWolfe(LineSearch):
Selected point.
fb : float
Value of polynomial at point b.
Returns
-------
xmin : float
Position of the approximated minimum.
Position of the approximated minimum.
"""
# f(x) = B*(x-a)^2 + C*(x-a) + D
with np.errstate(divide='raise', over='raise', invalid='raise'):
......
......@@ -116,7 +116,7 @@ class RGRGTransformation(Transformation):
np.absolute(np.array(domain.shape) *
np.array(domain.distances) *
np.array(codomain.distances) - 1) <
10**-7):
1e-7):
raise AttributeError("The grid-distances of domain and codomain "
"do not match.")
......
......@@ -93,11 +93,11 @@ class HPSpace(Space):
@property
def shape(self):
return (np.int(12 * self.nside ** 2),)
return (np.int(12 * self.nside * self.nside),)
@property
def dim(self):
return np.int(12 * self.nside ** 2)
return np.int(12 * self.nside * self.nside)
@property
def total_volume(self):
......@@ -108,7 +108,7 @@ class HPSpace(Space):
def weight(self, x, power=1, axes=None, inplace=False):
weight = ((4 * np.pi) / (12 * self.nside**2)) ** np.float(power)
weight = ((4*np.pi) / (12*self.nside*self.nside)) ** np.float(power)
if inplace:
x *= weight
......
......@@ -130,7 +130,7 @@ class LMSpace(Space):
# dim = (((2*(l+1)-1)+1)**2/4 - 2 * (l-m)(l-m+1)/2
# dim = np.int((l+1)**2 - (l-m)*(l-m+1.))
# We fix l == m
return np.int((l+1)**2)
return np.int((l+1)*(l+1))
@property
def total_volume(self):
......@@ -166,7 +166,7 @@ class LMSpace(Space):
def get_fft_smoothing_kernel_function(self, sigma):
# FIXME why x(x+1) ? add reference to paper!
return lambda x: np.exp(-0.5 * x * (x + 1) * sigma**2)
return lambda x: np.exp(-0.5 * x * (x + 1) * sigma*sigma)