jah wrote:
> Long post, but I hoped to be thorough.
>
> I'm having some trouble generalizing a simple Cython implementation
> that made feasible a previously infeasible Python implementation. The
> problem can be cast as follows:
>
> Let f(x,p) be a function, for example f(x,p) = sin(p*x).
> Pick a range of values for p, call this: pvalues
> Then, in pseudocode:
>
> for p in pvalues:
> xn = x0
> for i in range(n):
> # this could benefit from ticket #474
> xn = f(xn, p)
> # increment count for bin at (xn, p)
> # ...
> return bincounts
>
> I'm following the technique in the Cython tutorial that was written
> for Scipy 2009, and I have the following extension class defined:
>
> ctypedef np.float_t FTYPE_t
> cdef class MyFunction:
> cpdef FTYPE_t evaluate(self, FTYPE_t x)
>
> This gives me the flexibility to define the function in Python space
> as well (which is desirable). So I've coded this function as above,
> making sure to provide types for everything, with the following call
> signature:
>
> @cython.boundscheck(False)
> @cython.wraparound(False)
> def histogram(MyFunction func not None,
> np.ndarray[FTYPE_t, ndim=1] pvalues,
> FTYPE_t x0,
> int n,
> int bins):
>
> Then the fundamental line which exists in the inner loop is:
>
> xn = func.evaluate(xn, p)
When turning this into using arrays, consider having an out-argument
instead of returning something, to avoid memory allocation on every
iteration.
I.e. func.evaluate(xn, p, xn)
>
> Everything works and is quick. Now, I'm having trouble generalising
> this so that both x and p can be vectors of arbitrary (but fixed)
> length. My understanding says that I must go from cpdef to just def
> (due to ticket #177):
>
> ctypedef np.float_t FTYPE_t
> cdef class MyFunction:
> def evaluate(self, np.ndarray[FTYPE_t, ndim=1, mode="c"] x,
> np.ndarray[FTYPE_t, ndim=1, mode="c"] p):
>
> Of course, the problem now is that I have making a Python call inside
> a doubly nested forloop and things get much slower. I can see that p
> would be fixed in the inner loop, but I'd still have to deal with the
> fact that xn is changing each time.
>
> So how can I implement this quickly and generally? Thanks!
You basically have to stoop down to using C pointers:
cpdef evaluate(FTYPE_t *x, FTYPE_t *p, FTYPE_t* out, Py_ssize_t len):
# Access x[i], p[i]
Then call like this:
func.evaluate(<FTYPE_t*>x.data, <FTYPE_t*>p.data, ...)
Note that this relies on mode='c' bit for correctness!
Dag Sverre
_______________________________________________
Cython-dev mailing list
[email protected]
http://codespeak.net/mailman/listinfo/cython-dev
