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

first try

parent 574443a8
...@@ -13,7 +13,7 @@ before_script: ...@@ -13,7 +13,7 @@ before_script:
- apt-get update - apt-get update
- > - >
apt-get install -y build-essential python python-pip python-dev git 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 - pip install --upgrade -r ci/requirements_base.txt
- chmod +x ci/*.sh - chmod +x ci/*.sh
......
...@@ -15,7 +15,7 @@ Summary ...@@ -15,7 +15,7 @@ Summary
a versatile library designed to enable the development of signal a versatile library designed to enable the development of signal
inference algorithms that operate regardless of the underlying spatial inference algorithms that operate regardless of the underlying spatial
grid and its resolution. Its object-oriented framework is written in 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. efficiency.
NIFTY offers a toolkit that abstracts discretized representations of NIFTY offers a toolkit that abstracts discretized representations of
...@@ -71,7 +71,6 @@ Installation ...@@ -71,7 +71,6 @@ Installation
- [Python](http://www.python.org/) (v2.7.x) - [Python](http://www.python.org/) (v2.7.x)
- [NumPy](http://www.numpy.org/) - [NumPy](http://www.numpy.org/)
- [Cython](http://cython.org/)
### Download ### Download
...@@ -95,7 +94,7 @@ Starting with a fresh Ubuntu installation move to a folder like ...@@ -95,7 +94,7 @@ Starting with a fresh Ubuntu installation move to a folder like
- Using pip install numpy etc...: - Using pip install numpy etc...:
sudo pip install numpy cython sudo pip install numpy
- Install pyHealpix: - Install pyHealpix:
...@@ -147,10 +146,9 @@ MacPorts, missing ones need to be installed manually. It may also be ...@@ -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 mentioned that one should only use one package manager, as multiple ones
may cause trouble. may cause trouble.
- Install basic packages numpy and cython: - Install numpy:
sudo port install py27-numpy sudo port install py27-numpy
sudo port install py27-cython
- Install pyHealpix: - Install pyHealpix:
......
numpy numpy
cython
mpi4py mpi4py
matplotlib matplotlib
plotly plotly
......
#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,
endindex,
distances,
smooth_length,
smoothing_width)
outarr[tuple(i.tolist())] = res
k += 1
return outarr
...@@ -21,13 +21,136 @@ import numpy as np ...@@ -21,13 +21,136 @@ import numpy as np
import nifty.nifty_utilities as utilities import nifty.nifty_utilities as utilities
from nifty.operators.endomorphic_operator import EndomorphicOperator from nifty.operators.endomorphic_operator import EndomorphicOperator
from nifty.operators.fft_operator import FFTOperator from nifty.operators.fft_operator import FFTOperator
from nifty.operators.smoothing_operator import smooth_util as su #from nifty.operators.smoothing_operator import smooth_util as su
from d2o import STRATEGIES from d2o import STRATEGIES
def smoothingHelper (x,sigma,dxmax=None):
""" Performs the 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.
The indices go from ibegin[i] (included) to ibegin[i]+nval[i](excluded).
wgt0: integer array of the same size as x
cumulative sum of nval
wgt: float array with sum(nval) entries containing the normalized smoothing
weights.
The weights for grid point i start at wgt[wgt0[i]].
"""
if dxmax==None:
dxmax=3.01*sigma
x=np.asarray(x)
ibegin=np.searchsorted(x,x-dxmax)
nval=np.searchsorted(x,x+dxmax)-ibegin
wgt0=np.empty(x.size,dtype=np.int)
wgt0[0]=0
wgt0[1:x.size]=np.cumsum(nval[0:x.size-1])
wgt=np.empty(nval.sum(),dtype=np.float64)
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[wgt0[i]:wgt0[i]+nval[i]]=t
return ibegin,nval,wgt0,wgt
def apply_kernel_along_array(power, startindex, endindex, distances,
smooth_length, smoothing_width,ibegin,nval,wgt0,wgt):
if smooth_length == 0.0:
return power[startindex:endindex]
p_smooth = np.empty(endindex-startindex, dtype=np.float64)
for i in xrange(startindex, endindex):
p_smooth[i-startindex]=np.sum(power[ibegin[i]:ibegin[i]+nval[i]]*wgt[wgt0[i]:wgt0[i]+nval[i]])
return p_smooth
def getShape(a):
return tuple(a.shape)
def apply_along_axis(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,wgt0,wgt=smoothingHelper(distances,smooth_length,smooth_length*smoothing_width)
ind = np.zeros(nd-1, dtype=np.int)
i = np.zeros(nd, dtype=object)
shape = getShape(arr)
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 = apply_kernel_along_array(slicedArr,
startindex,
endindex,
distances,
smooth_length,
smoothing_width, ibegin,nval,wgt0,wgt)
outshape = np.asarray(getShape(arr))
outshape[axis] = endindex - startindex
outarr = np.zeros(outshape, dtype=np.float64)
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 = apply_kernel_along_array(slicedArr,
startindex,
endindex,
distances,
smooth_length,
smoothing_width,ibegin,nval,wgt0,wgt)
outarr[tuple(i.tolist())] = res
k += 1
return outarr
class SmoothingOperator(EndomorphicOperator): class SmoothingOperator(EndomorphicOperator):
# ---Overwritten properties and methods--- # ---Overwritten properties and methods---
def __init__(self, domain=(), sigma=0, log_distances=False, def __init__(self, domain, sigma, log_distances=False,
default_spaces=None): default_spaces=None):
super(SmoothingOperator, self).__init__(default_spaces) super(SmoothingOperator, self).__init__(default_spaces)
...@@ -234,9 +357,10 @@ class SmoothingOperator(EndomorphicOperator): ...@@ -234,9 +357,10 @@ class SmoothingOperator(EndomorphicOperator):
else: else:
true_sigma = self.sigma true_sigma = self.sigma
#MR: do not split for complex arrays!
if data.dtype is np.dtype('float32'): if data.dtype is np.dtype('float32'):
distances = distances.astype(np.float32, copy=False) distances = distances.astype(np.float32, copy=False)
smoothed_data = su.apply_along_axis_f( smoothed_data = apply_along_axis(
data_axis, data, data_axis, data,
startindex=true_start, startindex=true_start,
endindex=true_end, endindex=true_end,
...@@ -245,7 +369,7 @@ class SmoothingOperator(EndomorphicOperator): ...@@ -245,7 +369,7 @@ class SmoothingOperator(EndomorphicOperator):
smoothing_width=self._direct_smoothing_width) smoothing_width=self._direct_smoothing_width)
elif data.dtype is np.dtype('float64'): elif data.dtype is np.dtype('float64'):
distances = distances.astype(np.float64, copy=False) distances = distances.astype(np.float64, copy=False)
smoothed_data = su.apply_along_axis( smoothed_data = apply_along_axis(
data_axis, data, data_axis, data,
startindex=true_start, startindex=true_start,
endindex=true_end, endindex=true_end,
......
...@@ -30,7 +30,7 @@ class PowerSpace(Space): ...@@ -30,7 +30,7 @@ class PowerSpace(Space):
# ---Overwritten properties and methods--- # ---Overwritten properties and methods---
def __init__(self, harmonic_domain=RGSpace((1,)), def __init__(self, harmonic_domain,
distribution_strategy='not', distribution_strategy='not',
log=False, nbin=None, binbounds=None): log=False, nbin=None, binbounds=None):
......
...@@ -18,7 +18,6 @@ ...@@ -18,7 +18,6 @@
from setuptools import setup, find_packages from setuptools import setup, find_packages
from Cython.Build import cythonize
import numpy import numpy
exec(open('nifty/version.py').read()) exec(open('nifty/version.py').read())
...@@ -33,9 +32,6 @@ setup(name="ift_nifty", ...@@ -33,9 +32,6 @@ setup(name="ift_nifty",
packages=find_packages(), packages=find_packages(),
package_dir={"nifty": "nifty"}, package_dir={"nifty": "nifty"},
zip_safe=False, zip_safe=False,
ext_modules=cythonize([
"nifty/operators/smoothing_operator/smooth_util.pyx"
]),
include_dirs=[numpy.get_include()], include_dirs=[numpy.get_include()],
dependency_links=[ dependency_links=[
'git+https://gitlab.mpcdf.mpg.de/ift/keepers.git#egg=keepers-0.3.7', 'git+https://gitlab.mpcdf.mpg.de/ift/keepers.git#egg=keepers-0.3.7',
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment