Ask Your Question
1

compile a symbolic function and use it in cython

asked 13 years ago

pang gravatar image

updated 13 years ago

Jason Grout gravatar image

I have written code to create plots of Newton fractals. Cython was key, pure python was way too slow.

But I'd like to create an interact: write the polynomial and watch the fractal. The problem is the user should enter the polynomial in a text box as a symbolic expression, then I should send that to cython and the cython code should do the rest, calling that function when necessary.

I've managed to do so: I create a fast_callable on the expression, and pass the output of fast_callable to cython, but it is very slow. Unfortunately, I cannot make the different data types for complex numbers match, and there might be something else.

Maybe you can think of a different approach...

Preview: (hide)

2 Answers

Sort by » oldest newest most voted
1

answered 13 years ago

Jason Grout gravatar image

updated 13 years ago

kcrisman gravatar image

You probably want to use the fast_callable's call_c function. I've posted some examples at http://sagenb.org/home/pub/4086/, but I repeat them below for convenience:

%cython

from sage.ext.interpreters.wrapper_rdf cimport Wrapper_rdf
from libc.math cimport sqrt

def benchmark(float n):
    cdef int i
    for i in range(1000):
        n = sqrt(n)
    return n

def iterate_c(Wrapper_rdf f,float n):
    cdef double *result = [n]
    cdef int i
    for i in range(1000):
        f.call_c(result, result)
    return result[0]

def iterate_py(f, float n):
    cdef int i
    for i in range(1000):
        n=f(n)
    return n

Then

g(x)=sqrt(x)
f=fast_callable(g,vars=[x],domain=float)

Timings:

sage: timeit('benchmark(1e23)')
625 loops, best of 3: 7.15 µs per loop
sage: timeit('iterate_c(f,1e23)')
625 loops, best of 3: 18.3 µs per loop
sage: timeit('iterate_py(f,1e23)')
625 loops, best of 3: 162 µs per loop
sage: timeit('iterate_py(g,1e23)')
5 loops, best of 3: 77.2 ms per loop
Preview: (hide)
link

Comments

I can't for the life of me figure out what you updated, kcrisman. Is there a place I can see the history? (Edit: oh, I see; just "edit" and then I can see the revisions...you just updated a link. I was afraid I made a coding error.)

Jason Grout gravatar imageJason Grout ( 13 years ago )

Awesome! Repeating the trick with "double complex", "double_complex", and a rusty knowledge of C was not easy, though.

pang gravatar imagepang ( 13 years ago )

@pang: if you end up with something even better (better looking, faster, or both) for the Newton basins, please add it to #11837. Sage will probably always be slower on this sort of thing, but it doesn't have to be ''too'' slow.

kcrisman gravatar imagekcrisman ( 13 years ago )
1

answered 13 years ago

kcrisman gravatar image

I'm not sure if this answers your actual question, but you might be interested in code that Simon King and I have at # 11837 for generating Newton basin fractals. You are right that it is slow. Hopefully not quite as slow as yours, though, maybe it's all you need.

Incidentally, even at David Joyce's awesome generator page he just takes the roots; trying to get that from a generic polynomial first will severely slow you down, I'm sure.

Preview: (hide)
link

Comments

Well, asking for the roots instead of the polynomial does the trick for polynomials, but the question remains: how can I get a symbolic function (or some string representing a function) from an interact, compile it, and use it in cython with similar speed to a function that was written in cython.

pang gravatar imagepang ( 13 years ago )

That's why I wasn't sure if it answered your actual question; I wasn't sure if you wanted Newton basins specifically, or the cython. I don't know how to use symbolic things in Cython, though Ginac is C++ so maybe there is something possible.

kcrisman gravatar imagekcrisman ( 13 years ago )

Well, both things are great to know. Thanks!

pang gravatar imagepang ( 13 years ago )

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 13 years ago

Seen: 983 times

Last updated: Jan 18 '12