Wiener Filter advanced

.gitlab-ci.yml

@@ -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
 

Makefile

-# This Makefile implements common tasks needed by developers
-# A list of implemented rules can be obtained by the command "make help"
-
-.DEFAULT_GOAL=build
-.PHONY .SILENT : help
-help :
-	echo
-	echo "    Implemented targets:"
-	echo
-	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"
-	echo
-
-.PHONY : clean
-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 setuptools-3.3.zip
-
-	# remove dist/
-	rm -rf dist
-
-.PHONY : build
-build : build2
-
-.PHONY : build2
-build2 :
-	python2 setup.py build_ext --inplace
-
-.PHONY :
-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 setup.py 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 \
-	; \
-	done
-
-.PHONY : sdist
-sdist :
-	python setup.py 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 setup.py build
-
-.PHONY : build-sdist3
-build-sdist3 : ./dist/pypmc*/NUL
-	cd dist/pypmc* && python3 setup.py 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

PKG-INFO

-Metadata-Version: 1.0
-Name: ift_nifty
-Version: 1.0.6
-Summary: Numerical Information Field Theory
-Home-page: http://www.mpa-garching.mpg.de/ift/nifty/
-Author: Theo Steininger
-Author-email: theos@mpa-garching.mpg.de
-License: GPLv3
-Description: UNKNOWN
-Platform: UNKNOWN

README.md

@@ -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:
 

ci/requirements.txt

 numpy
-cython
 mpi4py
 matplotlib
 plotly

nifty/operators/__init__.py

@@ -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 *
 

nifty/operators/linear_operator/linear_operator.py

@@ -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

nifty/operators/response_operator/response_operator.py

@@ -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)
 
@@ -80,27 +81,22 @@ class ResponseOperator(LinearOperator):
             shape_target = np.append(shape_target, self._domain[ii].shape)
 
         self._target = self._parse_domain(FieldArray(shape_target))
-        self._sigma = sigma
-        self._exposure = exposure
-
-        self._kernel = len(self._domain)*[None]
-
-        for ii in xrange(len(self._kernel)):
-            self._kernel[ii] = SmoothingOperator(self._domain[ii],
-                                                 sigma=self._sigma[ii])
-
-        self._composed_kernel = ComposedOperator(self._kernel)
-
-        self._exposure_op = len(self._domain)*[None]
-        if len(self._exposure_op) != len(self._kernel):
-            raise ValueError("Definition of kernel and exposure do not suit "
-                             "each other")
-        else:
-            for ii in xrange(len(self._exposure_op)):
-                self._exposure_op[ii] = DiagonalOperator(
-                                                self._domain[ii],
-                                                diagonal=self._exposure[ii])
-            self._composed_exposure = ComposedOperator(self._exposure_op)
+
+        kernel_smoothing = len(self._domain)*[None]
+        kernel_exposure = len(self._domain)*[None]
+
+        if len(sigma)!= len(exposure):
+            raise ValueError("Length of smoothing kernel and length of"
+                             "exposure do not match")
+
+        for ii in xrange(len(kernel_smoothing)):
+            kernel_smoothing[ii] = SmoothingOperator(self._domain[ii],
+                                                     sigma=sigma[ii])
+            kernel_exposure[ii] = DiagonalOperator(self._domain[ii],
+                                                   diagonal=exposure[ii])
+
+        self._composed_kernel = ComposedOperator(kernel_smoothing)
+        self._composed_exposure = ComposedOperator(kernel_exposure)
 
     @property
     def domain(self):

nifty/operators/smoothing_operator/__init__.py

@@ -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

nifty/operators/smoothing_operator/direct_smoothing_operator.py

+# -*- 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

nifty/operators/smoothing_operator/fft_smoothing_operator.py

+# -*- 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

nifty/operators/smoothing_operator/smooth_util.pyx

-#cython: nonecheck=False
-#cython: boundscheck=True