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

Merge branch 'byebye_cython' into 'master'

Replace direct smoothing code by pure Python

See merge request !96
parents 7b445aed 260b19e1
Pipeline #12463 failed with stages
in 6 minutes and 28 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
......
......@@ -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
efficiency.
NIFTY offers a toolkit that abstracts discretized representations of
......@@ -71,7 +71,6 @@ Installation
- [Python](http://www.python.org/) (v2.7.x)
- [NumPy](http://www.numpy.org/)
- [Cython](http://cython.org/)
### 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:
......
numpy
cython
mpi4py
matplotlib
plotly
......
......@@ -24,7 +24,7 @@ from diagonal_operator import DiagonalOperator
from endomorphic_operator import EndomorphicOperator
from smoothing_operator import SmoothingOperator
from smoothing_operator import *
from fft_operator import *
......
......@@ -75,7 +75,8 @@ class LinearOperator(Loggable, object):
def __init__(self, default_spaces=None):
self.default_spaces = default_spaces
def _parse_domain(self, domain):
@staticmethod
def _parse_domain(domain):
return utilities.parse_domain(domain)
@abc.abstractproperty
......
......@@ -6,6 +6,7 @@ from nifty.operators.smoothing_operator import SmoothingOperator
from nifty.operators.composed_operator import ComposedOperator
from nifty.operators.diagonal_operator import DiagonalOperator
class ResponseOperator(LinearOperator):
""" NIFTy ResponseOperator (example)
......
......@@ -16,4 +16,4 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
from smoothing_operator import SmoothingOperator
from .smoothing_operator import SmoothingOperator
# -*- coding: utf8 -*-
import numpy as np
from d2o import STRATEGIES
from .smoothing_operator import SmoothingOperator
class DirectSmoothingOperator(SmoothingOperator):
def __init__(self, domain, sigma, log_distances=False,
default_spaces=None):
super(DirectSmoothingOperator, self).__init__(domain,
sigma,
log_distances,
default_spaces)
self.effective_smoothing_width = 3.01
def _precompute(self, x, sigma, dxmax=None):
""" Does precomputations for Gaussian smoothing on a 1D irregular grid.
Parameters
----------
x: 1D floating point array or list containing the individual grid
positions. Points must be given in ascending order.
sigma: The sigma of the Gaussian with which the function living on x
should be smoothed, in the same units as x.
dxmax: (optional) The maximum distance up to which smoothing is
performed, in the same units as x. Default is 3.01*sigma.
Returns
-------
ibegin: integer array of the same size as x
ibegin[i] is the minimum grid index to consider when computing the
smoothed value at grid index i
nval: integer array of the same size as x
nval[i] is the number of indices to consider when computing the
smoothed value at grid index i.
wgt: list with the same number of entries as x
wgt[i] is an array with nval[i] entries containing the
normalized smoothing weights.
"""
if dxmax is None:
dxmax = self.effective_smoothing_width*sigma
x = np.asarray(x)
ibegin = np.searchsorted(x, x-dxmax)
nval = np.searchsorted(x, x+dxmax) - ibegin
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)
return ibegin, nval, wgt
def _apply_kernel_along_array(self, power, startindex, endindex,
distances, smooth_length, smoothing_width,
ibegin, nval, wgt):
if smooth_length == 0.0:
return power[startindex:endindex]
p_smooth = np.zeros(endindex-startindex, dtype=power.dtype)
for i in xrange(startindex, endindex):
imin = max(startindex, ibegin[i])
imax = min(endindex, ibegin[i]+nval[i])
p_smooth[imin:imax] += (power[i] *
wgt[i][imin-ibegin[i]:imax-imin+ibegin[i]])
return p_smooth
def _apply_along_axis(self, axis, arr, startindex, endindex, distances,
smooth_length, smoothing_width):
nd = arr.ndim
if axis < 0:
axis += nd
if (axis >= nd):
raise ValueError(
"axis must be less than arr.ndim; axis=%d, rank=%d."
% (axis, nd))
ibegin, nval, wgt = self._precompute(
distances, smooth_length, smooth_length*smoothing_width)
ind = np.zeros(nd-1, dtype=np.int)
i = np.zeros(nd, dtype=object)
shape = arr.shape
indlist = np.asarray(range(nd))
indlist = np.delete(indlist, axis)
i[axis] = slice(None, None)
outshape = np.asarray(shape).take(indlist)
i.put(indlist, ind)
Ntot = np.product(outshape)
holdshape = outshape
slicedArr = arr[tuple(i.tolist())]
res = self._apply_kernel_along_array(
slicedArr, startindex, endindex, distances,
smooth_length, smoothing_width, ibegin, nval, wgt)
outshape = np.asarray(arr.shape)
outshape[axis] = endindex - startindex
outarr = np.zeros(outshape, dtype=arr.dtype)
outarr[tuple(i.tolist())] = res
k = 1
while k < Ntot:
# increment the index
ind[nd-1] += 1
n = -1
while (ind[n] >= holdshape[n]) and (n > (1-nd)):
ind[n-1] += 1
ind[n] = 0
n -= 1
i.put(indlist, ind)
slicedArr = arr[tuple(i.tolist())]
res = self._apply_kernel_along_array(
slicedArr, startindex, endindex, distances,
smooth_length, smoothing_width, ibegin, nval, wgt)
outarr[tuple(i.tolist())] = res
k += 1
return outarr
def _smooth(self, x, spaces, inverse):
# infer affected axes
# we rely on the knowledge, that `spaces` is a tuple with length 1.
affected_axes = x.domain_axes[spaces[0]]
if len(affected_axes) > 1:
raise ValueError("By this implementation only one-dimensional "
"spaces can be smoothed directly.")
affected_axis = affected_axes[0]
distance_array = x.domain[spaces[0]].get_distance_array(
distribution_strategy='not')
distance_array = distance_array.get_local_data(copy=False)
if self.log_distances:
np.log(distance_array, out=distance_array)
# collect the local data + ghost cells
local_data_Q = False
if x.distribution_strategy == 'not':
local_data_Q = True
elif x.distribution_strategy in STRATEGIES['slicing']:
# infer the local start/end based on the slicing information of
# x's d2o. Only gets non-trivial for axis==0.
if 0 != affected_axis:
local_data_Q = True
else:
start_index = x.val.distributor.local_start
start_distance = distance_array[start_index]
augmented_start_distance = \
(start_distance -
self.effective_smoothing_width*self.sigma)
augmented_start_index = \
np.searchsorted(distance_array, augmented_start_distance)
true_start = start_index - augmented_start_index
end_index = x.val.distributor.local_end
end_distance = distance_array[end_index-1]
augmented_end_distance = \
(end_distance + self.effective_smoothing_width*self.sigma)
augmented_end_index = \
np.searchsorted(distance_array, augmented_end_distance)
true_end = true_start + x.val.distributor.local_length
augmented_slice = slice(augmented_start_index,
augmented_end_index)
augmented_data = x.val.get_data(augmented_slice,
local_keys=True,
copy=False)
augmented_data = augmented_data.get_local_data(copy=False)
augmented_distance_array = distance_array[augmented_slice]
else:
raise ValueError("Direct smoothing not implemented for given"
"distribution strategy.")
if local_data_Q:
# if the needed data resides on the nodes already, the necessary
# are the same; no matter what the distribution strategy was.
augmented_data = x.val.get_local_data(copy=False)
augmented_distance_array = distance_array
true_start = 0
true_end = x.shape[affected_axis]
# perform the convolution along the affected axes
# currently only one axis is supported
data_axis = affected_axes[0]
if inverse:
true_sigma = 1. / self.sigma
else:
true_sigma = self.sigma
local_result = self._apply_along_axis(
data_axis,
augmented_data,
startindex=true_start,
endindex=true_end,
distances=augmented_distance_array,
smooth_length=true_sigma,
smoothing_width=self.effective_smoothing_width)
result = x.copy_empty()
result.val.set_local_data(local_result, copy=False)
return result
# -*- coding: utf-8 -*-
import numpy as np
from nifty.operators.fft_operator import FFTOperator
from .smoothing_operator import SmoothingOperator
class FFTSmoothingOperator(SmoothingOperator):
def _smooth(self, x, spaces, inverse):
Transformator = FFTOperator(x.domain[spaces[0]])
# transform to the (global-)default codomain and perform all remaining
# steps therein
transformed_x = Transformator(x, spaces=spaces)
codomain = transformed_x.domain[spaces[0]]
coaxes = transformed_x.domain_axes[spaces[0]]
# create the kernel using the knowledge of codomain about itself
axes_local_distribution_strategy = \
transformed_x.val.get_axes_local_distribution_strategy(axes=coaxes)
kernel = codomain.get_distance_array(
distribution_strategy=axes_local_distribution_strategy)
if self.log_distances:
kernel.apply_scalar_function(np.log, inplace=True)
kernel.apply_scalar_function(
codomain.get_fft_smoothing_kernel_function(self.sigma),
inplace=True)
# now, apply the kernel to transformed_x
# this is done node-locally utilizing numpys reshaping in order to
# apply the kernel to the correct axes
local_transformed_x = transformed_x.val.get_local_data(copy=False)
local_kernel = kernel.get_local_data(copy=False)
reshaper = [transformed_x.shape[i] if i in coaxes else 1
for i in xrange(len(transformed_x.shape))]
local_kernel = np.reshape(local_kernel, reshaper)
# apply the kernel
if inverse:
local_transformed_x /= local_kernel
else:
local_transformed_x *= local_kernel
transformed_x.val.set_local_data(local_transformed_x, copy=False)
smoothed_x = Transformator.adjoint_times(transformed_x, spaces=spaces)
result = x.copy_empty()
result.set_val(smoothed_x, copy=False)
return result
#cython: nonecheck=False
#cython: boundscheck=True
#cython: wraparound=False
#cython: cdivision=True
import numpy as np
cimport numpy as np
#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
cdef int SIGMA_RANGE = 3
cpdef np.float32_t gaussianKernel_f(np.ndarray[np.float32_t, ndim=1] mpower,
np.ndarray[np.float32_t, ndim=1] mk,
np.float32_t mu,
np.float32_t smooth_length):
cdef np.ndarray[np.float32_t, ndim=1] C = \
np.exp(-(mk - mu) ** 2 / (2. * smooth_length ** 2))
return np.sum(C * mpower) / np.sum(C)
cpdef np.float64_t gaussianKernel(np.ndarray[np.float64_t, ndim=1] mpower,
np.ndarray[np.float64_t, ndim=1] mk,
np.float64_t mu,
np.float64_t smooth_length):
cdef np.ndarray[np.float64_t, ndim=1] C = \
np.exp(-(mk - mu) ** 2 / (2. * smooth_length ** 2))
return np.sum(C * mpower) / np.sum(C)
cpdef np.ndarray[np.float64_t, ndim=1] apply_kernel_along_array(
np.ndarray[np.float64_t, ndim=1] power,
np.int_t startindex,
np.int_t endindex,
np.ndarray[np.float64_t, ndim=1] distances,
np.float64_t smooth_length,
np.float64_t smoothing_width):
if smooth_length == 0.0:
return power[startindex:endindex]
if (smooth_length is None) or (smooth_length < 0):
smooth_length = distances[1]-distances[0]
cdef np.ndarray[np.float64_t, ndim=1] p_smooth = \
np.empty(endindex-startindex, dtype=np.float64)
cdef np.int64_t i, l, u
for i in xrange(startindex, endindex):
current_location = distances[i]
l = np.searchsorted(distances,
current_location - smoothing_width*smooth_length)
u = np.searchsorted(distances,
current_location + smoothing_width*smooth_length)
p_smooth[i-startindex] = gaussianKernel(power[l:u],
distances[l:u],
distances[i],
smooth_length)
return p_smooth
cpdef np.ndarray[np.float32_t, ndim=1] apply_kernel_along_array_f(
np.ndarray[np.float32_t, ndim=1] power,
np.int_t startindex,
np.int_t endindex,
np.ndarray[np.float32_t, ndim=1] distances,
np.float64_t smooth_length,
np.float64_t smoothing_width):
if smooth_length == 0.0:
return power[startindex:endindex]
if (smooth_length is None) or (smooth_length < 0):
smooth_length = distances[1]-distances[0]
cdef np.ndarray[np.float32_t, ndim=1] p_smooth = \
np.empty(endindex-startindex, dtype=np.float64)
cdef np.int64_t i, l, u
for i in xrange(startindex, endindex):
current_location = distances[i]
l = np.searchsorted(distances,
current_location - smoothing_width*smooth_length)
u = np.searchsorted(distances,
current_location + smoothing_width*smooth_length)
p_smooth[i-startindex] = gaussianKernel(power[l:u],
distances[l:u],
distances[i],
smooth_length)
return p_smooth
def getShape(a):
return tuple(a.shape)
cpdef np.ndarray[np.float64_t, ndim=1] apply_along_axis(
np.int_t axis,
np.ndarray arr,
np.int_t startindex,
np.int_t endindex,
np.ndarray[np.float64_t, ndim=1] distances,
np.float64_t smooth_length,
np.float64_t smoothing_width):
cdef np.int_t nd = arr.ndim
if axis < 0:
axis += nd
if (axis >= nd):
raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
% (axis, nd))
cdef np.ndarray[np.int_t, ndim=1] ind = np.zeros(nd-1, dtype=np.int)
cdef np.ndarray i = np.zeros(nd, dtype=object)
cdef tuple shape = getShape(arr)
cdef np.ndarray[np.int_t, ndim=1] indlist = np.asarray(range(nd))
indlist = np.delete(indlist,axis)
i[axis] = slice(None, None)
cdef np.ndarray[np.int_t, ndim=1] outshape = \
np.asarray(shape).take(indlist)
i.put(indlist, ind)
cdef np.ndarray[np.float64_t, ndim=1] slicedArr
cdef np.ndarray[np.float64_t, ndim=1] res
cdef np.int_t Ntot = np.product(outshape)
cdef np.ndarray[np.int_t, ndim=1] holdshape = outshape
slicedArr = arr[tuple(i.tolist())]
res = apply_kernel_along_array(slicedArr,
startindex,
endindex,
distances,
smooth_length,
smoothing_width)
outshape = np.asarray(getShape(arr))
outshape[axis] = endindex - startindex
outarr = np.zeros(outshape, dtype=np.float64)
outarr[tuple(i.tolist())] = res
cdef np.int_t n, k = 1
while k < Ntot:
# increment the index
ind[nd-1] += 1
n = -1
while (ind[n] >= holdshape[n]) and (n > (1-nd)):
ind[n-1] += 1
ind[n] = 0
n -= 1
i.put(indlist, ind)
slicedArr = arr[tuple(i.tolist())]
res = apply_kernel_along_array(slicedArr,
startindex,
endindex,
distances,
smooth_length,
smoothing_width)
outarr[tuple(i.tolist())] = res
k += 1
return outarr
cpdef np.ndarray[np.float32_t, ndim=1] apply_along_axis_f(
np.int_t axis,
np.ndarray arr,
np.int_t startindex,
np.int_t endindex,
np.ndarray[np.float32_t, ndim=1] distances,
np.float64_t smooth_length,
np.float64_t smoothing_width):
cdef np.int_t nd = arr.ndim
if axis < 0:
axis += nd
if (axis >= nd):
raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
% (axis, nd))
cdef np.ndarray[np.int_t, ndim=1] ind = np.zeros(nd-1, dtype=np.int)
cdef np.ndarray i = np.zeros(nd, dtype=object)
cdef tuple shape = getShape(arr)
cdef np.ndarray[np.int_t, ndim=1] indlist = np.asarray(range(nd))
indlist = np.delete(indlist,axis)
i[axis] = slice(None, None)
cdef np.ndarray[np.int_t, ndim=1] outshape = \
np.asarray(shape).take(indlist)
i.put(indlist, ind)
cdef np.ndarray[np.float32_t, ndim=1] slicedArr
cdef np.ndarray[np.float32_t, ndim=1] res
cdef np.int_t Ntot = np.product(outshape)
cdef np.ndarray[np.int_t, ndim=1] holdshape = outshape
slicedArr = arr[tuple(i.tolist())]
res = apply_kernel_along_array_f(slicedArr,
startindex,
endindex,
distances,
smooth_length,
smoothing_width)
outshape = np.asarray(getShape(arr))
outshape[axis] = endindex - startindex
outarr = np.zeros(outshape, dtype=np.float32)
outarr[tuple(i.tolist())] = res
cdef np.int_t n, k = 1
while k < Ntot:
# increment the index
ind[nd-1] += 1
n = -1
while (ind[n] >= holdshape[n]) and (n > (1-nd)):
ind[n-1] += 1
ind[n] = 0
n -= 1
i.put(indlist, ind)
slicedArr = arr[tuple(i.tolist())]
res = apply_kernel_along_array(slicedArr,
startindex,