Revision: 4662
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4662&view=rev
Author: jdh2358
Date: 2007-12-06 20:43:36 -0800 (Thu, 06 Dec 2007)
Log Message:
-----------
updates to the pyrex examples
Modified Paths:
--------------
trunk/py4science/examples/pyrex/movavg/Makefile
trunk/py4science/examples/pyrex/movavg/movavg.pyx
trunk/py4science/examples/pyrex/movavg/ringbuf.pyx
trunk/py4science/examples/pyrex/movavg/ringbuf_demo.py
trunk/py4science/examples/pyrex/movavg/setup.py
Added Paths:
-----------
trunk/py4science/examples/pyrex/movavg/README
trunk/py4science/examples/pyrex/movavg/movavg_bruteforce.py
trunk/py4science/examples/pyrex/movavg/movavg_fast.py
trunk/py4science/examples/pyrex/movavg/movavg_ringbuf.py
trunk/py4science/examples/pyrex/movavg/ringbufnan.c
Removed Paths:
-------------
trunk/py4science/examples/pyrex/movavg/circbuffer.py
trunk/py4science/examples/pyrex/movavg/movavg_test.py
trunk/py4science/examples/pyrex/movavg/numpyx.pyx
trunk/py4science/examples/pyrex/movavg/run_test.py
Modified: trunk/py4science/examples/pyrex/movavg/Makefile
===================================================================
--- trunk/py4science/examples/pyrex/movavg/Makefile 2007-12-07 00:01:46 UTC
(rev 4661)
+++ trunk/py4science/examples/pyrex/movavg/Makefile 2007-12-07 04:43:36 UTC
(rev 4662)
@@ -6,4 +6,4 @@
.PHONY: clean
clean:
- rm -rf *~ *.so *.c *.o build
+ rm -rf *~ *.so *.o build ringbuf.c
Added: trunk/py4science/examples/pyrex/movavg/README
===================================================================
--- trunk/py4science/examples/pyrex/movavg/README
(rev 0)
+++ trunk/py4science/examples/pyrex/movavg/README 2007-12-07 04:43:36 UTC
(rev 4662)
@@ -0,0 +1,82 @@
+Introduction
+============
+
+This exercise introduces pyrex to wrap a C library for trailing
+statistics.
+
+Computation of trailing windowed statistics is common in many
+quantitative data driven disciplines, particularly where there is
+noisy data. Common uses of windowed statistics are the trailing
+moving average, standard deviation, minumum and maximum. Two common
+use cases which pose computational challenges for python: real time
+updating of trailing statistics as live data comes in, and posthoc
+computation of trailing statistics over a large data array. In the
+second case, for some statistics we can use convolution and related
+techniques for efficient computation, eg of the trailing 30 sample
+average
+
+ numpy.convolve(x, numpy.ones(30), mode=valid')[:len(x)]
+
+but for other statistics like the trailing 30 day maximum at each
+point, efficient routines like convolution are of no help.
+
+This exercise introduces pyrex to efficiently solve the problem of
+trailing statistics over arrays as well as for a live, incoming data
+stream. A pure C library, ringbuf, defines a circular C buffer and
+attached methods for efficiently computing trailing averages, and
+pyrex is used to provide a pythonic API on top of this extension code.
+The rigid segregation between the C library and the python wrappers
+insures that the C code can be used in other projects, be it a matlab
+(TM) extension or some other C library. The goal of the exericse is
+to compute the trailing statistics mean, median, stddev, min and max
+using three approaches:
+
+ - with brute force using numpy arrays, slices and methods
+ (movavg_bruteforce.py)
+
+ - with python bindings to the ringbuf code ringbuf.Rinbuf
+ (movavg_ringbuf.py). See ringbuf_demo.py for an example of how to
+ use the ringbuf module
+
+ - using a pyrex extension to the ringbuf runstats code
+ (movavg_fast.py)
+
+pyrex module support
+====================
+
+ - Makefile : simple interface to setup.py so you can just 'make'
+
+ - setup.py : configure and build the python modules
+
+ - c_numpy.pxd : the numpy C API for pyrex
+
+ - c_python.pxd : the python C API for pyrex
+
+ - c_ringbuf.pxi : the ringbuf C API for pyrex
+
+ - ringbuf.pyx : python interface to the ringbuf C API
+
+examples
+========
+
+ - ringbuf_demo.py : basic demo of the python bindings to basic
+ Ringbuf class
+
+ - movavg_bruteforce.py : do the trailing stats with brute force
+ numpy slices and methods
+
+ - movavg_ringbuf.py : do the trailing stats with the Ringbuf code
+
+ - movavg_fast.py : do the trailing stats with the ringbuf runstats wrapper
+
+ringbuf C code
+==============
+
+ - ringbuf.h : pure C ringbuf API headers
+
+ - ringbufnan.c : pure C ringbuf library
+
+Acknowledgements
+================
+
+Thanks to Eric Firing for the ringbuf and runstats code!
Deleted: trunk/py4science/examples/pyrex/movavg/circbuffer.py
===================================================================
--- trunk/py4science/examples/pyrex/movavg/circbuffer.py 2007-12-07
00:01:46 UTC (rev 4661)
+++ trunk/py4science/examples/pyrex/movavg/circbuffer.py 2007-12-07
04:43:36 UTC (rev 4662)
@@ -1,20 +0,0 @@
-import numpy as npy
-
-class CircularBuffer:
- def __init__(self, N):
- self.i = 0
- self.data = npy.zeros(N)
- self.N = N
-
- def add(self, val):
- self.data[self.i%self.N] = val
- self.i += 1
-
- def get_data(self):
- if self.i<self.N:
- return self.data[:self.i]
- else:
- return self.data
-
-
-
Modified: trunk/py4science/examples/pyrex/movavg/movavg.pyx
===================================================================
--- trunk/py4science/examples/pyrex/movavg/movavg.pyx 2007-12-07 00:01:46 UTC
(rev 4661)
+++ trunk/py4science/examples/pyrex/movavg/movavg.pyx 2007-12-07 04:43:36 UTC
(rev 4662)
@@ -1,41 +1,55 @@
-# -*- Mode: Python -*- Not really, but close enough
+'''
+This is kept as a module separate from pyringbuf so that the
+latter does not depend on numpy.
+'''
-cimport c_python
cimport c_numpy
import numpy
-# Numpy must be initialized
c_numpy.import_array()
+include "c_ringbuf.pxi"
-def sum_elements(c_numpy.ndarray arr):
- cdef int i
- cdef double x, val
+def runstats(data, nrb):
+ '''
+ Compute running stats on 1D array data for odd length nrb
+ '''
- x = 0.
- val = 0.
- for i from 0<=i<arr.dimensions[0]:
- val = (<double*>(arr.data + i*arr.strides[0]))[0]
- x = x + val
+ cdef c_numpy.ndarray c_data
+ cdef c_numpy.ndarray c_dmean
+ cdef c_numpy.ndarray c_dstd
+ cdef c_numpy.ndarray c_dmin
+ cdef c_numpy.ndarray c_dmax
+ cdef c_numpy.ndarray c_dmedian
+ cdef c_numpy.ndarray c_ng
- return x
+ data = numpy.asarray(data, dtype=numpy.float_)
+ if data.ndim != 1:
+ raise ValueError("data must be 1-D for now")
+ nd = data.shape[0]
+ dmean = numpy.empty_like(data)
+ dstd = numpy.empty_like(data)
+ dmin = numpy.empty_like(data)
+ dmax = numpy.empty_like(data)
+ dmedian = numpy.empty_like(data)
+ ng = numpy.empty(data.shape, dtype=numpy.int_)
-def sum_elements2(c_numpy.ndarray arr):
- cdef int i
- cdef double x, val
+ c_data = data
+ c_dmean = dmean
+ c_dstd = dstd
+ c_dmin = dmin
+ c_dmax = dmax
+ c_dmedian = dmedian
+ c_ng = ng
- arr = numpy.asarray(arr, numpy.float_)
-
- if arr.nd!=1:
- raise RuntimeError('only 1D arrays supported; found
shape=%s'%str(arr.shape))
- assert(arr.nd==1)
- x = 0.
- val = 0.
- for i from 0<=i<arr.dimensions[0]:
- val = (<double*>(arr.data + i*arr.strides[0]))[0]
- x = x + val
-
- return x
+ c_runstats(nrb, nd, <double *>c_data.data,
+ <double *>c_dmean.data,
+ <double *>c_dstd.data,
+ <double *>c_dmin.data,
+ <double *>c_dmax.data,
+ <double *>c_dmedian.data,
+ <int *>c_ng.data)
+ return dmean, dstd, dmin, dmax, dmedian, ng
Added: trunk/py4science/examples/pyrex/movavg/movavg_bruteforce.py
===================================================================
--- trunk/py4science/examples/pyrex/movavg/movavg_bruteforce.py
(rev 0)
+++ trunk/py4science/examples/pyrex/movavg/movavg_bruteforce.py 2007-12-07
04:43:36 UTC (rev 4662)
@@ -0,0 +1,14 @@
+"""
+Use a brute force, slow method
+"""
+import numpy
+
+x = numpy.random.rand(10000)
+
+data = []
+for i in range(len(x)):
+ slice = x[max(0, i-30):i+1] # slice out the recent history
+ data.append([slice.mean(), slice.std(), slice.min(), slice.max(),
numpy.median(slice), len(slice)])
+
+r = numpy.rec.fromarrays(data,names='dmean,dstd,dmin,dmax,dmedian,ngood')
+
Added: trunk/py4science/examples/pyrex/movavg/movavg_fast.py
===================================================================
--- trunk/py4science/examples/pyrex/movavg/movavg_fast.py
(rev 0)
+++ trunk/py4science/examples/pyrex/movavg/movavg_fast.py 2007-12-07
04:43:36 UTC (rev 4662)
@@ -0,0 +1,14 @@
+"""
+Use a brute force, slow method
+"""
+import numpy
+import ringbuf
+
+x = numpy.random.rand(10000)
+dmean, dstd, dmin, dmax, dmedian, ng = ringbuf.runstats(x, 30)
+
+r = numpy.rec.fromarrays([dmean, dstd, dmin, dmax, dmedian, ng],
+ names='dmean,dstd,dmin,dmax,dmedian,ngood')
+
+
+
Added: trunk/py4science/examples/pyrex/movavg/movavg_ringbuf.py
===================================================================
--- trunk/py4science/examples/pyrex/movavg/movavg_ringbuf.py
(rev 0)
+++ trunk/py4science/examples/pyrex/movavg/movavg_ringbuf.py 2007-12-07
04:43:36 UTC (rev 4662)
@@ -0,0 +1,16 @@
+"""
+Use a brute force, slow method
+"""
+import numpy
+import ringbuf
+
+r = ringbuf.Ringbuf(30)
+x = numpy.random.rand(10000)
+
+data = []
+for thisx in x:
+ r.add(thisx)
+ data.append([thisx, r.N_good(), r.min(), r.max(), r.mean(), r.sd()])
+
+r = numpy.rec.fromarrays(data,names='x,ngood,min30,max30,mean30,sd30')
+
Deleted: trunk/py4science/examples/pyrex/movavg/movavg_test.py
===================================================================
--- trunk/py4science/examples/pyrex/movavg/movavg_test.py 2007-12-07
00:01:46 UTC (rev 4661)
+++ trunk/py4science/examples/pyrex/movavg/movavg_test.py 2007-12-07
04:43:36 UTC (rev 4662)
@@ -1,3 +0,0 @@
-import numpy
-import movavg
-
Deleted: trunk/py4science/examples/pyrex/movavg/numpyx.pyx
===================================================================
--- trunk/py4science/examples/pyrex/movavg/numpyx.pyx 2007-12-07 00:01:46 UTC
(rev 4661)
+++ trunk/py4science/examples/pyrex/movavg/numpyx.pyx 2007-12-07 04:43:36 UTC
(rev 4662)
@@ -1,128 +0,0 @@
-# -*- Mode: Python -*- Not really, but close enough
-
-cimport c_python
-cimport c_numpy
-import numpy
-
-# Numpy must be initialized
-c_numpy.import_array()
-
-def print_array_info(c_numpy.ndarray arr):
- cdef int i
-
- print '-='*10
- print 'printing array info for ndarray at
0x%0lx'%(<c_python.Py_intptr_t>arr,)
- print 'print number of dimensions:',arr.nd
- print 'address of strides: 0x%0lx'%(<c_python.Py_intptr_t>arr.strides,)
- print 'strides:'
- for i from 0<=i<arr.nd:
- # print each stride
- print ' stride %d:'%i,<c_python.Py_intptr_t>arr.strides[i]
- print 'memory dump:'
- print_elements( arr.data, arr.strides, arr.dimensions,
- arr.nd, sizeof(double), arr.dtype )
- print '-='*10
- print
-
-
-
-def sum_elements(c_numpy.ndarray arr):
- cdef int i
- cdef double x, val
-
- x = 0.
- val = 0.
- for i from 0<=i<arr.dimensions[0]:
- val = (<double*>(arr.data + i*arr.strides[0]))[0]
- x = x + val
-
- return x
-
-
-def scale_elements(int N):
- cdef int i
- cdef double x, val
-
- x = 0.
- val = 0.
- for i from 0<=i<N:
- val = 2.5 * i
- x = x + val
- return x
-
-
-cdef print_elements(char *data,
- c_python.Py_intptr_t* strides,
- c_python.Py_intptr_t* dimensions,
- int nd,
- int elsize,
- object dtype):
- cdef c_python.Py_intptr_t i,j
- cdef void* elptr
-
- if dtype not in [numpy.dtype(numpy.object_),
- numpy.dtype(numpy.float64)]:
- print ' print_elements() not (yet) implemented for dtype
%s'%dtype.name
- return
-
- if nd ==0:
- if dtype==numpy.dtype(numpy.object_):
- elptr = (<void**>data)[0] #[0] dereferences pointer in Pyrex
- print ' ',<object>elptr
- elif dtype==numpy.dtype(numpy.float64):
- print ' ',(<double*>data)[0]
- elif nd == 1:
- for i from 0<=i<dimensions[0]:
- if dtype==numpy.dtype(numpy.object_):
- elptr = (<void**>data)[0]
- print ' ',<object>elptr
- elif dtype==numpy.dtype(numpy.float64):
- print ' ',(<double*>data)[0]
- data = data + strides[0]
- else:
- for i from 0<=i<dimensions[0]:
- print_elements(data, strides+1, dimensions+1, nd-1, elsize, dtype)
- data = data + strides[0]
-
-
-def test_methods(c_numpy.ndarray arr):
- """Test a few attribute accesses for an array.
-
- This illustrates how the pyrex-visible object is in practice a strange
- hybrid of the C PyArrayObject struct and the python object. Some
- properties (like .nd) are visible here but not in python, while others
- like flags behave very differently: in python flags appears as a separate,
- object while here we see the raw int holding the bit pattern.
-
- This makes sense when we think of how pyrex resolves arr.foo: if foo is
- listed as a field in the c_numpy.ndarray struct description, it will be
- directly accessed as a C variable without going through Python at all.
- This is why for arr.flags, we see the actual int which holds all the flags
- as bit fields. However, for any other attribute not listed in the struct,
- it simply forwards the attribute lookup to python at runtime, just like
- python would (which means that AttributeError can be raised for
- non-existent attributes, for example)."""
-
- print 'arr.any() :',arr.any()
- print 'arr.nd :',arr.nd
- print 'arr.flags :',arr.flags
-
-def test():
- """this function is pure Python"""
- arr1 = numpy.array(-1e-30,dtype=numpy.float64)
- arr2 = numpy.array([1.0,2.0,3.0],dtype=numpy.float64)
-
- arr3 = numpy.arange(9,dtype=numpy.float64)
- arr3.shape = 3,3
-
- four = 4
- arr4 = numpy.array(['one','two',3,four],dtype=numpy.object_)
-
- arr5 = numpy.array([1,2,3]) # int types not (yet) supported by
print_elements
-
- for arr in [arr1,arr2,arr3,arr4,arr5]:
- print_array_info(arr)
-
-
- arr = numpy.arange(10.0)
- print 'sum el', sum_elements(arr)
Modified: trunk/py4science/examples/pyrex/movavg/ringbuf.pyx
===================================================================
--- trunk/py4science/examples/pyrex/movavg/ringbuf.pyx 2007-12-07 00:01:46 UTC
(rev 4661)
+++ trunk/py4science/examples/pyrex/movavg/ringbuf.pyx 2007-12-07 04:43:36 UTC
(rev 4662)
@@ -1,11 +1,16 @@
-''' pyringbuf.pyx Pyrex interface to ringbuf.c
- It defines a single ring buffer class, Ringbuf,
- on which various statistics are calculated
- as each entry is added.
-'''
+"""
+ringbuf.pyx Pyrex interface to ringbuf.c It defines ring buffer class,
+Ringbuf, on which various statistics are calculated as each entry is
+added and a method runstats for computing a host of trailing
+statistics over a numpy array
+"""
include "c_ringbuf.pxi"
+# ringbuf is a C library, but we can give it a pythonic object
+# oriented API by managing the ringbuf pointers and function calls in
+# a python class
+
cdef class Ringbuf:
cdef ringbuf_t *rb_ptr
@@ -132,3 +137,72 @@
+# now we'll provide a numpy interface to the runstats function, which
+# fills predimensioned arrays with the trailing mean, std, max, etc,
+# for a data array
+cimport c_numpy
+import numpy
+c_numpy.import_array()
+
+def runstats(data, nrb):
+ """
+ Compute running stats on 1D array data over window length nrb
+ """
+
+ # we will be calling the ringbuf C function runstats, which
+ # expects a bunch of c data pointers to the data arrays as well as
+ # the output arrays for the mean, std, etc... Create a array for
+ # each C float array pointer the C function expects. The type of
+ # the array is c_numpy.ndarray
+ cdef c_numpy.ndarray c_data
+ cdef c_numpy.ndarray c_dmean
+ cdef c_numpy.ndarray c_dstd
+ cdef c_numpy.ndarray c_dmin
+ cdef c_numpy.ndarray c_dmax
+ cdef c_numpy.ndarray c_dmedian
+ cdef c_numpy.ndarray c_ng
+
+ # make sure that the input array is a 1D numpy array of floats.
+ # asarray is used to copy and cast a python sequence or array to
+ # the approriate type, with the advantage that if the data is
+ # already the right type, no copy is performed
+ data = numpy.asarray(data, dtype=numpy.float_)
+ if data.ndim != 1:
+ raise ValueError("data must be 1-D for now")
+ nd = data.shape[0]
+
+ # use numpy.empty_like to create and empty array of the same type
+ # and shape as data for each of the return stats (mean, std, ...)
+ # These are genuine numpy arrays we will be passing back to the
+ # python level. Remember the ng (the number of good elements) is
+ # an int, and everything else is a float
+ dmean = numpy.empty_like(data)
+ dstd = numpy.empty_like(data)
+ dmin = numpy.empty_like(data)
+ dmax = numpy.empty_like(data)
+ dmedian = numpy.empty_like(data)
+ ng = numpy.empty(data.shape, dtype=numpy.int_)
+
+ # now we have to assign the c_data structures and friends to their
+ # corresponding numpy arrays
+ c_data = data
+ c_dmean = dmean
+ c_dstd = dstd
+ c_dmin = dmin
+ c_dmax = dmax
+ c_dmedian = dmedian
+ c_ng = ng
+
+ # now we call the function and pass in the c data pointers to the
+ # arrays. The syntax <double *>c_data.data tells pyrex to pass
+ # the numpy data memory block as a pointer to a float array.
+ c_runstats(nrb, nd, <double *>c_data.data,
+ <double *>c_dmean.data,
+ <double *>c_dstd.data,
+ <double *>c_dmin.data,
+ <double *>c_dmax.data,
+ <double *>c_dmedian.data,
+ <int *>c_ng.data)
+
+ # all done, return the arrays
+ return dmean, dstd, dmin, dmax, dmedian, ng
Modified: trunk/py4science/examples/pyrex/movavg/ringbuf_demo.py
===================================================================
--- trunk/py4science/examples/pyrex/movavg/ringbuf_demo.py 2007-12-07
00:01:46 UTC (rev 4661)
+++ trunk/py4science/examples/pyrex/movavg/ringbuf_demo.py 2007-12-07
04:43:36 UTC (rev 4662)
@@ -1,2 +1,14 @@
+"""
+Example code showin how to use the ringbuf extension code from python
+"""
+import random
import ringbuf
+r = ringbuf.Ringbuf(30)
+
+for i in range(100):
+ r.add(random.random())
+ print 'Nadded=%d, Ngood=%d, min=%1.2f, max=%1.2f, mean=%1.2f, std=%1.2f'%(
+ r.N_added(), r.N_good(), r.min(), r.max(), r.mean(), r.sd())
+
+print 'slice', r[:10]
Added: trunk/py4science/examples/pyrex/movavg/ringbufnan.c
===================================================================
--- trunk/py4science/examples/pyrex/movavg/ringbufnan.c
(rev 0)
+++ trunk/py4science/examples/pyrex/movavg/ringbufnan.c 2007-12-07 04:43:36 UTC
(rev 4662)
@@ -0,0 +1,363 @@
+/*
+
+ Use a ring buffer to calculate running statistics of a data stream.
+
+
+
+ compile:
+ gcc -Wall -c -fPIC ringbufnan.c
+
+
+ 2003/07/28 EF
+
+*/
+
+#include <stdlib.h>
+#include <stdio.h> /* for debugging printf statements */
+#include <math.h>
+#include "ringbuf.h"
+
+static double NaN = 0.0;
+
+static void sort_ringbuf(ringbuf_t *rb_ptr);
+static void resum_ringbuf(ringbuf_t *rb_ptr);
+
+ringbuf_t *new_ringbuf(int N)
+{
+ if (!isnan(NaN))
+ {
+ NaN = strtod("NaN", NULL); /* May change for Matlab. */
+ }
+ ringbuf_t *rb_ptr;
+ rb_ptr = calloc(1, sizeof(ringbuf_t));
+ rb_ptr->i_sorted = calloc(N, sizeof(int));
+ rb_ptr->data = calloc(N, sizeof(double));
+ rb_ptr->N_size = N;
+ zero_ringbuf(rb_ptr);
+ return rb_ptr;
+}
+
+void zero_ringbuf(ringbuf_t *rb_ptr)
+{
+ rb_ptr->N_filled = 0;
+ rb_ptr->N_good = 0;
+ rb_ptr->N_added = 0;
+ rb_ptr->i_oldest = 0;
+ rb_ptr->i_next = 0;
+ rb_ptr->sum = 0.0;
+ rb_ptr->sumsq = 0.0;
+}
+
+void delete_ringbuf(ringbuf_t *rb_ptr)
+{
+ free(rb_ptr->data);
+ free(rb_ptr->i_sorted);
+ free(rb_ptr);
+}
+
+int ringbuf_index(ringbuf_t *rb_ptr, int i)
+/* i is a Python-style index; return the actual offset, or -1 for error */
+{
+ /* printf("incoming index: %d", i);*/
+ if (i >= rb_ptr->N_filled) return -1;
+ if (i < -rb_ptr->N_filled) return -1;
+ if (i < 0) i += rb_ptr->i_next;
+ else i += rb_ptr->i_oldest;
+ i %= rb_ptr->N_size;
+ if (i < 0) i += rb_ptr->N_size;
+ /*printf("changed to offset: %d\n", i);*/
+ return i;
+}
+
+int ringbuf_slice_i(ringbuf_t *rb_ptr, int i)
+/* For slices, out of range indices get clipped.*/
+/* This function is not needed for the pyrex interface. */
+{
+ /*printf("incoming index: %d", i);*/
+ if (i >= rb_ptr->N_filled) i = rb_ptr->N_filled;
+ if (i < -rb_ptr->N_filled) i = 0;
+ if (i < 0) i += rb_ptr->i_next;
+ else i += rb_ptr->i_oldest;
+ i %= rb_ptr->N_size;
+ if (i < 0) i += rb_ptr->N_size;
+ /*printf("changed to offset: %d\n", i);*/
+ return i;
+}
+
+double ringbuf_getitem(ringbuf_t *rb_ptr, int i)
+{
+ i = ringbuf_index(rb_ptr, i);
+ if (i < 0) return 1e38;
+ return rb_ptr->data[i];
+}
+
+
+void resum_ringbuf(ringbuf_t *rb_ptr)
+{
+ int i;
+ double d;
+ rb_ptr->sum = 0.0;
+ rb_ptr->sumsq = 0.0;
+ for (i = 0; i < rb_ptr->N_filled; i++)
+ {
+ d = rb_ptr->data[i];
+ if (!isnan(d))
+ {
+ rb_ptr->sum += d;
+ rb_ptr->sumsq += d*d;
+ }
+ }
+}
+
+void ringbuf_add(ringbuf_t *rb_ptr, double d)
+{
+
+ double d_old;
+ int i, i_new, good_new, N;
+
+ N = rb_ptr->N_size; /* We need this many times. */
+ i_new = rb_ptr->i_next;
+ /* Save the old value; otherwise, it will be overwritten. */
+ d_old = rb_ptr->data[rb_ptr->i_oldest];
+ rb_ptr->data[i_new] = d;
+ good_new = !isnan(d);
+#if 0
+ printf("new value: %lf good_new: %d\n", d, good_new);
+ printf("i_next: %d i_oldest: %d N_filled: %d N_good: %d\n",
+ rb_ptr->i_next, rb_ptr->i_oldest,
+ rb_ptr->N_filled, rb_ptr->N_good);
+#endif
+ rb_ptr->i_next++;
+ rb_ptr->i_next %= N;
+
+ if (rb_ptr->N_filled == rb_ptr->N_size) /* already full */
+ {
+ if (!isnan(d_old))
+ {
+ rb_ptr->sum -= d_old;
+ rb_ptr->sumsq -= d_old * d_old;
+ /* Remove the i_oldest entry from i_sorted,
+ and slide the remaining entries down.
+ */
+ for (i = 0; i < rb_ptr->N_good; i++)
+ {
+ if (rb_ptr->i_sorted[i] == rb_ptr->i_oldest)
+ break;
+ }
+ for ( ; i < rb_ptr->N_good - 1; i++)
+ {
+ rb_ptr->i_sorted[i] = rb_ptr->i_sorted[i+1];
+ }
+ rb_ptr->N_good--;
+ }
+ /* The new entry is temporarily put at the end of the array.
+ */
+ if (good_new)
+ {
+ rb_ptr->i_sorted[rb_ptr->N_good] = i_new;
+ rb_ptr->N_good++;
+ }
+ rb_ptr->i_oldest++;
+ rb_ptr->i_oldest %= N;
+ }
+ else /* not full before adding this element */
+ {
+ if (good_new)
+ {
+ rb_ptr->i_sorted[rb_ptr->N_good] = i_new;
+ rb_ptr->N_good++;
+ }
+ rb_ptr->N_filled++;
+ }
+ if (good_new)
+ {
+ rb_ptr->sum += d;
+ rb_ptr->sumsq += d*d;
+ sort_ringbuf(rb_ptr);
+ }
+ /* To prevent accumulation of truncation error, we
+ recalculate the sums periodically.
+ */
+ rb_ptr->N_added++;
+ if (rb_ptr->N_added % (rb_ptr->N_size + 10000) == 0)
+ {
+ resum_ringbuf(rb_ptr);
+ }
+#if 0
+ printf("i_next: %d i_oldest: %d N_filled: %d N_good: %d\n",
+ rb_ptr->i_next, rb_ptr->i_oldest,
+ rb_ptr->N_filled, rb_ptr->N_good);
+#endif
+}
+
+/* This is not a full sort--it assumes the list is
+ already sorted except for the last entry. It uses a single
+ pass of bubble sorting to put that entry in its place.
+ The code could be moved bodily into the function above.
+*/
+void sort_ringbuf(ringbuf_t *rb_ptr)
+{
+ int i, i_hold;
+ int *ip;
+ double *dp;
+
+ ip = rb_ptr->i_sorted;
+ dp = rb_ptr->data;
+
+ for (i = rb_ptr->N_good - 1; i > 0; i--)
+ {
+ if (dp[ip[i]] < dp[ip[i-1]])
+ {
+ i_hold = ip[i];
+ ip[i] = ip[i-1];
+ ip[i-1] = i_hold;
+ }
+ else
+ {
+ break;
+ }
+ }
+}
+
+double ringbuf_min(ringbuf_t *rb_ptr)
+{
+ return rb_ptr->data[rb_ptr->i_sorted[0]];
+}
+
+double ringbuf_max(ringbuf_t *rb_ptr)
+{
+ int i_end;
+ i_end = rb_ptr->N_good - 1;
+ return rb_ptr->data[rb_ptr->i_sorted[i_end]];
+}
+
+double ringbuf_median(ringbuf_t *rb_ptr)
+{
+ int i_mid, N;
+ N = rb_ptr->N_good;
+ if (N == 0) return NaN;
+ i_mid = N/2;
+ if (N % 2 == 1)
+ {
+ return rb_ptr->data[rb_ptr->i_sorted[i_mid]];
+ }
+ else
+ {
+ return 0.5 * (rb_ptr->data[rb_ptr->i_sorted[i_mid]]
+ + rb_ptr->data[rb_ptr->i_sorted[i_mid - 1]]);
+ }
+}
+
+int ringbuf_N_good(ringbuf_t *rb_ptr)
+{
+ return rb_ptr->N_good;
+}
+
+
+int ringbuf_N_filled(ringbuf_t *rb_ptr)
+{
+ return rb_ptr->N_filled;
+}
+
+
+
+int ringbuf_N_added(ringbuf_t *rb_ptr)
+{
+ return rb_ptr->N_added;
+}
+
+
+
+double ringbuf_mean(ringbuf_t *rb_ptr)
+{
+ int N;
+
+ N = rb_ptr->N_good;
+ if (N > 0)
+ {
+ return rb_ptr->sum / N;
+ }
+ else
+ {
+ return NaN;
+ }
+}
+
+double ringbuf_sd(ringbuf_t *rb_ptr)
+{
+ double m, s;
+ int N;
+
+ N = rb_ptr->N_good;
+ if (N == 0)
+ return NaN;
+ m = ringbuf_mean(rb_ptr);
+ s = rb_ptr->sumsq / N - m*m;
+ if (s > 0.0)
+ {
+ return sqrt(s);
+ }
+ else
+ {
+ return 0.0;
+ }
+}
+
+void c_runstats(int nrb, int nd, double *data, double *dmean, double *dstd,
+ double *dmin, double *dmax, double *dmed, int *ng)
+{
+ int i, j;
+ ringbuf_t *rb_ptr;
+
+ rb_ptr = new_ringbuf(nrb);
+
+ for (j=0; j<nd; i++, j++)
+ {
+ ringbuf_add(rb_ptr, data[j]);
+ dmean[j] = ringbuf_mean(rb_ptr);
+ dstd[j] = ringbuf_sd(rb_ptr);
+ dmin[j] = ringbuf_min(rb_ptr);
+ dmax[j] = ringbuf_max(rb_ptr);
+ dmed[j] = ringbuf_median(rb_ptr);
+ ng[j] = rb_ptr->N_good;
+ }
+ delete_ringbuf(rb_ptr);
+}
+
+
+void c_runstats2(int nrb, int nd, int step, int ofs,
+ double *data, double *dmean, double *dstd,
+ double *dmin, double *dmax, double *dmed, int *ng)
+{
+ int i, j;
+ int npad = (nrb - 1) / 2;
+ ringbuf_t *rb_ptr;
+
+ data += ofs;
+ dmean += ofs;
+ dstd += ofs;
+ dmin += ofs;
+ dmax += ofs;
+ dmed += ofs;
+ ng += ofs;
+
+ rb_ptr = new_ringbuf(nrb);
+ for (i = 0; i < npad; i++)
+ {
+ ringbuf_add(rb_ptr, data[i*step]);
+ }
+ for (j=0; j<nd; i++, j++)
+ {
+ if (i < nd) {ringbuf_add(rb_ptr, data[i*step]);}
+ else {ringbuf_add(rb_ptr, NaN);}
+ dmean[j*step] = ringbuf_mean(rb_ptr);
+ dstd[j*step] = ringbuf_sd(rb_ptr);
+ dmin[j*step] = ringbuf_min(rb_ptr);
+ dmax[j*step] = ringbuf_max(rb_ptr);
+ dmed[j*step] = ringbuf_median(rb_ptr);
+ ng[j*step] = rb_ptr->N_good;
+ }
+ delete_ringbuf(rb_ptr);
+}
+
+
Deleted: trunk/py4science/examples/pyrex/movavg/run_test.py
===================================================================
--- trunk/py4science/examples/pyrex/movavg/run_test.py 2007-12-07 00:01:46 UTC
(rev 4661)
+++ trunk/py4science/examples/pyrex/movavg/run_test.py 2007-12-07 04:43:36 UTC
(rev 4662)
@@ -1,3 +0,0 @@
-#!/usr/bin/env python
-from numpyx import test
-test()
Modified: trunk/py4science/examples/pyrex/movavg/setup.py
===================================================================
--- trunk/py4science/examples/pyrex/movavg/setup.py 2007-12-07 00:01:46 UTC
(rev 4661)
+++ trunk/py4science/examples/pyrex/movavg/setup.py 2007-12-07 04:43:36 UTC
(rev 4662)
@@ -18,16 +18,12 @@
ringbufmod = Extension('ringbuf',
['ringbuf.pyx',
- 'ringbufnan.c'])
+ 'ringbufnan.c'],
+ include_dirs = [numpy.get_include()])
-
-movavgmod = Extension('movavg',
- ['movavg.pyx'],
- include_dirs = [numpy.get_include()])
-
# Call the routine which does the real work
-setup(name = 'movavg',
- description = 'moveng averages, using a c ring buffer',
- ext_modules = [movavgmod, ringbufmod,],
+setup(name = 'ringbuf',
+ description = 'trailing stats using a c ring buffer',
+ ext_modules = [ringbufmod,],
cmdclass = {'build_ext': build_ext},
)
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
-------------------------------------------------------------------------
SF.Net email is sponsored by:
Check out the new SourceForge.net Marketplace.
It's the best place to buy or sell services for
just about anything Open Source.
http://sourceforge.net/services/buy/index.php
_______________________________________________
Matplotlib-checkins mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/matplotlib-checkins
