diff --git a/.gitignore b/.gitignore
index 9f55d0e4..7b57a1b1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,6 +5,10 @@ __pycache__/
# C extensions
*.so
+# Cython
+*.c
+*.html
+
# Distribution / packaging
.Python
env/
diff --git a/postpic/cythonfunctions.pyx b/postpic/cythonfunctions.pyx
new file mode 100644
index 00000000..ad2f9aae
--- /dev/null
+++ b/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 .
+#
+# 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[x] += 1.0
+ else:
+ ret[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 = (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
diff --git a/setup.py b/setup.py
index 75cef8dd..567b810b 100755
--- a/setup.py
+++ b/setup.py
@@ -16,6 +16,7 @@
# along with postpic. If not, see .
#
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',
diff --git a/test/test_cythonfunctions.py b/test/test_cythonfunctions.py
new file mode 100644
index 00000000..8489c135
--- /dev/null
+++ b/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 .
+#
+# 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()