From d452351e5fd19f76193b8a5973666c6bdc52812a Mon Sep 17 00:00:00 2001 From: Stephan Kuschel Date: Tue, 10 Mar 2015 16:06:42 +0100 Subject: [PATCH] add histogram function written in cython simple tests are added as well --- .gitignore | 4 ++ postpic/cythonfunctions.pyx | 73 ++++++++++++++++++++++++++++++++++++ setup.py | 4 +- test/test_cythonfunctions.py | 65 ++++++++++++++++++++++++++++++++ 4 files changed, 145 insertions(+), 1 deletion(-) create mode 100644 postpic/cythonfunctions.pyx create mode 100644 test/test_cythonfunctions.py 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()