add histogram function written in cython

simple tests are added as well

.gitignore

@@ -5,6 +5,10 @@ __pycache__/
 # C extensions
 *.so
 
+# Cython
+*.c
+*.html
+
 # Distribution / packaging
 .Python
 env/

postpic/cythonfunctions.pyx

@@ -0,0 +1,73 @@
+#
+# This file is part of postpic.
+#
+# postpic is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# postpic is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with postpic. If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright Stephan Kuschel 2015
+'''
+This file adds some has some runtime critical funtions implemented in cython.
+'''
+cimport cython
+
+import numpy as np
+cimport numpy as np
+
+
+def histogram(np.ndarray[np.double_t, ndim=1] data, range=None, int bins=20,
+           np.ndarray[np.double_t, ndim=1] weights=None, int order=0):
+    '''
+    Mimics numpy.histogram.
+    Additional Arguments:
+        - order = 0:
+            sets the order of the particle shapes.
+            order = 0 returns a normal histogram.
+            order = 1 uses top hat particle shape.
+    '''    
+    cdef np.ndarray[np.double_t, ndim=1] ret = np.zeros(bins, dtype=np.double);
+    cdef double min, max
+    if range is None:
+        min = np.min(data)
+        max = np.max(data)
+    else:
+        min = range[0]
+        max = range[1]
+    bin_edges = np.linspace(min, max, bins+1)
+    cdef int n = len(data)
+    cdef double tmp = 1.0 / (max - min) * bins
+    cdef double x
+    cdef int xr
+    if order == 0:
+        # normal Histogram
+        for i in xrange(n):
+            x = (data[i] - min) * tmp;
+            if x > 0.0 and x < bins:
+                if weights is None:
+                    ret[<int>x] += 1.0
+                else:
+                    ret[<int>x] += weights[i]
+    elif order == 1:
+        # Particle shape is spline of order 1 = TopHat
+        for i in xrange(n):
+            x = (data[i] - min) * tmp;
+            xr = <int>(x + 0.5);
+            if (xr >= 0.0 and xr < bins):
+                if weights is None:
+                    ret[xr] += (0.5 + x - xr) * 1.0
+                    if (xr > 1.0):
+                        ret[xr - 1] += (0.5 - x + xr) * 1.0
+                else:
+                    ret[xr] += (0.5 + x - xr) * weights[i]
+                    if (xr > 1.0):
+                        ret[xr - 1] += (0.5 - x + xr) * weights[i]
+    return ret, bin_edges

setup.py

@@ -16,6 +16,7 @@
 # along with postpic. If not, see <http://www.gnu.org/licenses/>.
 #
 from setuptools import setup
+from Cython.Build import cythonize
 
 setup(name='postpic',
       version='0.1.1',
@@ -24,8 +25,9 @@
       description='The open source particle-in-cell post processor.',
       url='http://github.com/skuschel/postpic',
       packages=['postpic'],
+      ext_modules = cythonize("postpic/cythonfunctions.pyx"),
       license='GPLv3+',
-      install_requires=['matplotlib', 'numpy>=1.7', 'scipy'],
+      install_requires=['matplotlib', 'numpy>=1.7', 'scipy', 'cython'],
       classifiers=[
           'Development Status :: 3 - Alpha',
           'Intended Audience :: Science/Research',

test/test_cythonfunctions.py

@@ -0,0 +1,65 @@
+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+#
+# This file is part of postpic.
+#
+# postpic is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# postpic is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with postpic. If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright Stephan Kuschel, 2015
+
+import unittest
+import postpic.cythonfunctions as cf
+import numpy as np
+
+class TestCythonfunctions(unittest.TestCase):
+
+    def setUp(self):
+        self.data = np.random.random(1e3)
+        self.weights = np.random.random(1e3)
+
+    def test_histogram(self):
+        # in this special case, cf.histogram and np.histogram must yield equal results
+        # check hist and edges
+        cfh, cfe = cf.histogram(self.data,bins=20, range=(0,1))
+        nph, npe = np.histogram(self.data, bins=20, range=(0,1))
+        self.assertListEqual(list(nph), list(cfh))
+        self.assertListEqual(list(npe), list(cfe))
+
+    def test_histogramw(self):
+        # in this special case, cf.histogram and np.histogram must yield equal results
+        # check hist and edges
+        cfh, cfe = cf.histogram(self.data,bins=100, range=(0,1), weights=self.weights)
+        nph, npe = np.histogram(self.data, bins=100, range=(0,1), weights=self.weights)
+        diff = np.abs(cfh - nph)
+        self.assertAlmostEqual(np.sum(diff), 0)
+        self.assertListEqual(list(npe), list(cfe))
+
+    def test_histogramo1(self):
+        cfh, cfe = cf.histogram(self.data,bins=100, range=(-0.1,1.1), order=1)
+        nph, npe = np.histogram(self.data, bins=100, range=(-0.1,1.1))
+        # just check that no mass is lost
+        diff = np.sum(cfh) - np.sum(nph)
+        self.assertAlmostEqual(diff, 0)
+        self.assertListEqual(list(npe), list(cfe))
+
+    def test_histogramo1w(self):
+        cfh, cfe = cf.histogram(self.data,bins=100, range=(-0.1,1.1), order=1, weights=self.weights)
+        nph, npe = np.histogram(self.data, bins=100, range=(-0.1,1.1), weights=self.weights)
+        # just check that no mass is lost
+        diff = np.sum(cfh) - np.sum(nph)
+        self.assertAlmostEqual(diff, 0)
+        self.assertListEqual(list(npe), list(cfe))
+
+if __name__ == '__main__':
+    unittest.main()