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

merge master

parents 11c94b91 a220fe56
Pipeline #13349 passed with stage
in 5 minutes and 17 seconds
......@@ -28,17 +28,25 @@ __all__ = ['dependency_injector', 'nifty_configuration']
# Setup the dependency injector
dependency_injector = keepers.DependencyInjector(
[('mpi4py.MPI', 'MPI'),
('pyfftw', 'fftw'),
'pyHealpix',
'plotly'])
dependency_injector.register('pyfftw', lambda z: hasattr(z, 'FFTW_MPI'))
def _fft_module_checker(z):
if z == 'fftw_mpi':
return hasattr(dependency_injector.get('fftw'), 'FFTW_MPI')
if z == 'fftw':
return ('fftw' in dependency_injector)
if z == 'numpy':
return True
return False
# Initialize the variables
variable_fft_module = keepers.Variable(
'fft_module',
['fftw', 'numpy'],
lambda z: (('pyfftw' in dependency_injector)
if z == 'fftw' else True))
['fftw_mpi', 'fftw', 'numpy'],
_fft_module_checker)
def dtype_validator(dtype):
......
......@@ -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'):
......
......@@ -63,9 +63,10 @@ class FFTOperator(LinearOperator):
but for full control, the user should explicitly specify a codomain.
module: String (optional)
Software module employed for carrying out the transform operations.
For RGSpace pairs this can be "numpy" or "fftw", where "numpy" is
always available, but "fftw" offers higher performance and
parallelization. For sphere-related domains, only "pyHealpix" is
For RGSpace pairs this can be "scalar" or "mpi", where "scalar" is
always available (using pyfftw if available, else numpy.fft), and "mpi"
requires pyfftw and offers MPI parallelization.
For sphere-related domains, only "pyHealpix" is
available. If omitted, "fftw" is selected for RGSpaces if available,
else "numpy"; on the sphere the default is "pyHealpix".
domain_dtype: data type (optional)
......
......@@ -25,7 +25,7 @@ import nifty.nifty_utilities as utilities
from keepers import Loggable
pyfftw = gdi.get('pyfftw')
fftw = gdi.get('fftw')
class Transform(Loggable, object):
......@@ -200,20 +200,21 @@ class Transform(Loggable, object):
raise NotImplementedError
class FFTW(Transform):
class MPIFFT(Transform):
"""
The pyfftw pendant of a fft object.
The MPI-parallel FFTW pendant of a fft object.
"""
def __init__(self, domain, codomain):
if 'pyfftw' not in gdi:
raise ImportError("The module pyfftw is needed but not available.")
if not hasattr(fftw, 'FFTW_MPI'):
raise ImportError(
"The MPI FFTW module is needed but not available.")
super(FFTW, self).__init__(domain, codomain)
super(MPIFFT, self).__init__(domain, codomain)
# Enable caching for pyfftw.interfaces
pyfftw.interfaces.cache.enable()
# Enable caching
fftw.interfaces.cache.enable()
# The plan_dict stores the FFTWTransformInfo objects which correspond
# to a certain set of (field_val, domain, codomain) sets.
......@@ -409,7 +410,7 @@ class FFTW(Transform):
def transform(self, val, axes, **kwargs):
"""
The pyfftw transform function.
The MPI-parallel FFTW transform function.
Parameters
----------
......@@ -467,8 +468,9 @@ class FFTW(Transform):
class FFTWTransformInfo(object):
def __init__(self, domain, codomain, axes, local_shape,
local_offset_Q, fftw_context, **kwargs):
if pyfftw is None:
raise ImportError("The module pyfftw is needed but not available.")
if not hasattr(fftw, 'FFTW_MPI'):
raise ImportError(
"The MPI FFTW module is needed but not available.")
shape = (local_shape if axes is None else
[y for x, y in enumerate(local_shape) if x in axes])
......@@ -512,9 +514,9 @@ class FFTWLocalTransformInfo(FFTWTransformInfo):
fftw_context,
**kwargs)
if codomain.harmonic:
self._fftw_interface = pyfftw.interfaces.numpy_fft.fftn
self._fftw_interface = fftw.interfaces.numpy_fft.fftn
else:
self._fftw_interface = pyfftw.interfaces.numpy_fft.ifftn
self._fftw_interface = fftw.interfaces.numpy_fft.ifftn
@property
def fftw_interface(self):
......@@ -531,7 +533,7 @@ class FFTWMPITransfromInfo(FFTWTransformInfo):
local_offset_Q,
fftw_context,
**kwargs)
self._plan = pyfftw.create_mpi_plan(
self._plan = fftw.create_mpi_plan(
input_shape=transform_shape,
input_dtype='complex128',
output_dtype='complex128',
......@@ -545,15 +547,26 @@ class FFTWMPITransfromInfo(FFTWTransformInfo):
return self._plan
class NUMPYFFT(Transform):
class SerialFFT(Transform):
"""
The numpy fft pendant of a fft object.
"""
def __init__(self, domain, codomain, use_fftw):
super(SerialFFT, self).__init__(domain, codomain)
if use_fftw and (fftw is None):
raise ImportError(
"The serial FFTW module is needed but not available.")
self._use_fftw = use_fftw
# Enable caching
if self._use_fftw:
fftw.interfaces.cache.enable()
def transform(self, val, axes, **kwargs):
"""
The pyfftw transform function.
The scalar FFT transform function.
Parameters
----------
......@@ -572,6 +585,7 @@ class NUMPYFFT(Transform):
result : np.ndarray or distributed_data_object
Fourier-transformed pendant of the input field.
"""
# Check if the axes provided are valid given the shape
if axes is not None and \
not all(axis in range(len(val.shape)) for axis in axes):
......@@ -625,10 +639,18 @@ class NUMPYFFT(Transform):
local_val = self._apply_mask(temp_val, mask, axes)
# perform the transformation
if self.codomain.harmonic:
result_val = np.fft.fftn(local_val, axes=axes)
if self._use_fftw:
if self.codomain.harmonic:
result_val = fftw.interfaces.numpy_fft.fftn(
local_val, axes=axes)
else:
result_val = fftw.interfaces.numpy_fft.ifftn(
local_val, axes=axes)
else:
result_val = np.fft.ifftn(local_val, axes=axes)
if self.codomain.harmonic:
result_val = np.fft.fftn(local_val, axes=axes)
else:
result_val = np.fft.ifftn(local_val, axes=axes)
# Apply domain centering mask
if reduce(lambda x, y: x + y, self.domain.zerocenter):
......
......@@ -18,7 +18,7 @@
import numpy as np
from transformation import Transformation
from rg_transforms import FFTW, NUMPYFFT
from rg_transforms import MPIFFT, SerialFFT
from nifty import RGSpace, nifty_configuration
......@@ -30,20 +30,18 @@ class RGRGTransformation(Transformation):
super(RGRGTransformation, self).__init__(domain, codomain, module)
if module is None:
if nifty_configuration['fft_module'] == 'fftw':
self._transform = FFTW(self.domain, self.codomain)
elif nifty_configuration['fft_module'] == 'numpy':
self._transform = NUMPYFFT(self.domain, self.codomain)
else:
raise ValueError('Unsupported default FFT module:' +
nifty_configuration['fft_module'])
module = nifty_configuration['fft_module']
if module == 'fftw_mpi':
self._transform = MPIFFT(self.domain, self.codomain)
elif module == 'fftw':
self._transform = SerialFFT(self.domain, self.codomain,
use_fftw=True)
elif module == 'numpy':
self._transform = SerialFFT(self.domain, self.codomain,
use_fftw=False)
else:
if module == 'fftw':
self._transform = FFTW(self.domain, self.codomain)
elif module == 'numpy':
self._transform = NUMPYFFT(self.domain, self.codomain)
else:
raise ValueError('Unsupported FFT module:' + module)
raise ValueError('Unsupported FFT module:' + module)
# ---Mandatory properties and methods---
......@@ -118,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.")
......
......@@ -53,10 +53,13 @@ class DirectSmoothingOperator(SmoothingOperator):
wgt = []
expfac = 1. / (2. * sigma*sigma)
for i in range(x.size):
t = x[ibegin[i]:ibegin[i]+nval[i]]-x[i]
t = np.exp(-t*t*expfac)
t *= 1./np.sum(t)
wgt.append(t)
if nval[i]>0:
t = x[ibegin[i]:ibegin[i]+nval[i]]-x[i]
t = np.exp(-t*t*expfac)
t *= 1./np.sum(t)
wgt.append(t)
else:
wgt.append(np.array([]))
return ibegin, nval, wgt
......@@ -146,7 +149,7 @@ class DirectSmoothingOperator(SmoothingOperator):
#MR FIXME: this causes calls of log(0.) which should probably be avoided
if self.log_distances:
np.log(distance_array, out=distance_array)
np.log(np.maximum(distance_array,1e-15), out=distance_array)
# collect the local data + ghost cells
local_data_Q = False
......
......@@ -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)
# ---Added properties and methods---
......
......@@ -267,14 +267,16 @@ class RGSpace(Space):
cords = np.ogrid[inds]
dists = ((cords[0] - shape[0]//2)*dk[0])**2
dists = (cords[0] - shape[0]//2)*dk[0]
dists *= dists
# apply zerocenterQ shift
if not self.zerocenter[0]:
dists = np.fft.ifftshift(dists)
# only save the individual slice
dists = dists[slice_of_first_dimension]
for ii in range(1, len(shape)):
temp = ((cords[ii] - shape[ii] // 2) * dk[ii])**2
temp = (cords[ii] - shape[ii] // 2) * dk[ii]
temp *= temp
if not self.zerocenter[ii]:
temp = np.fft.ifftshift(temp)
dists = dists + temp
......@@ -282,7 +284,7 @@ class RGSpace(Space):
return dists
def get_fft_smoothing_kernel_function(self, sigma):
return lambda x: np.exp(-0.5 * np.pi**2 * x**2 * sigma**2)
return lambda x: np.exp(-0.5 * np.pi*np.pi * x*x * sigma*sigma)
# ---Added properties and methods---
......
......@@ -16,9 +16,9 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from nose_parameterized import parameterized
from parameterized import parameterized
from nifty import RGSpace, LMSpace, HPSpace, GLSpace, PowerSpace
from nifty.config import dependency_injector as di
from nifty.config import dependency_injector as gdi
def custom_name_func(testcase_func, param_num, param):
......@@ -36,7 +36,7 @@ def expand(*args, **kwargs):
def generate_spaces():
spaces = [RGSpace(4), PowerSpace(RGSpace((4, 4), harmonic=True)),
LMSpace(5), HPSpace(4)]
if 'pyHealpix' in di:
if 'pyHealpix' in gdi:
spaces.append(GLSpace(4))
return spaces
......
from nifty import *
#This tests if it is possible to import all of nifties methods. Experience shows this is not always possible.
pass
\ No newline at end of file
# This tests if it is possible to import all of Nifty's methods.
# Experience shows this is not always possible.
pass
......@@ -20,7 +20,7 @@ import unittest
import numpy as np
from numpy.testing import assert_equal,\
assert_allclose
from nifty.config import dependency_injector as di
from nifty.config import dependency_injector as gdi
from nifty import Field,\
RGSpace,\
LMSpace,\
......@@ -62,11 +62,15 @@ class FFTOperatorTests(unittest.TestCase):
res = foo.get_distance_array('not')