ufuncify now creates actual ufuncs

[James Crist]

I still need to do some cleanups and add tests, but I finally have this working and thought I'd share. I'm really happy with this:

In [1]: from sympy import *

In [2]: a, b, c = symbols('a, b, c')

In [3]: expr = (sin(a) + sqrt(b)*c**2)/2

In [4]: from sympy.utilities.autowrap import ufuncify

In [5]: func = ufuncify((a, b, c), expr)

In [6]: func(1, 2, 3)
Out[6]: 6.7846965230828769

In [7]: func([1, 2, 3, 4, 5], [6, 7, 8, 9, 10], 3)
Out[7]: array([ 11.44343933,  12.36052961,  12.79848207,  13.12159875,  13.75078733])

In [8]: from numpy import arange

In [9]: a = arange(10).reshape((2, 5))

In [10]: c = arange(10, 20).reshape((2, 5))

In [11]: b = 25

In [12]: func(a, b, c)
Out[12]:
array([[ 250.        ,  302.92073549,  360.45464871,  422.57056   ,
         489.62159875],
       [ 562.02053786,  639.86029225,  722.8284933 ,  810.49467912,
         902.70605924]])

In [13]: type(func)
Out[13]: numpy.ufunc

This now does everything a numpy `ufunc` does normally, as it *is* a ufunc. Codegen is hooked up to numpy api. Type conversion and broadcasting are done automagically.

Caveats: only functions with a single output are accepted (this could be changed to accept multi-output without much effort though). Also, as with all unfuncs, input/outputs must all be scalars (no matrix/Indexed operations allowed).

[Matthew Rocklin]

Cool

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/76e0fbbe-5ce4-43b7-855b-6ac821f6b8ae%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jason Moore]

Awesome. I was working on this today but it looks like you've by passed what I had working. Do you have a PR with this?

Jason
moorepants.info
+01 530-601-9791

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAJ8oX-EHZXbd5aFFNRy7gJ0hcydpAsG2qxv7Py65DQ9cA9VUUA%40mail.gmail.com.

[James Crist]

Not yet. I wrote it this morning during an extremely boring meeting, and haven't had a chance to clean it up. This doesn't solve your problem about broadcasting a matrix calculation though...

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/azVZHLOv9Vc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAP7f1AieaeoOFtc_S4XPxWOX2jr2zmda9VCRpWpzHMTGLkmHPQ%40mail.gmail.com.

[Jason Moore]

Yeh, but if you simply create a ufunc for each expression in a matrix you still get substantial speedups. I wrote a bunch of test cases that I'll post to my blog tomorrow.

Jason
moorepants.info
+01 530-601-9791

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAJ2L7mfL_xO%3DO-ZRMx-zfpZzJKJ-%2BUdTzSCz5jYf%2B%3DdovR%2B_7Q%40mail.gmail.com.

[James Crist]

I was wondering about that. I wasn't sure if the overhead from looping through the inputs multiple times would outweigh improvements from fast C loops. Glad that in your case it does.

I've thrown a WIP PR up: https://github.com/sympy/sympy/pull/7929

For some reason, creating the functions in python with numpy calls still seems to be faster (for micro-benchmarks). This probably has something to do with function complexity (the example function above is simple), but I'd still think it'd be faster in pure C. I tried inlining the call, which was a small improvement, but it was still slower than the pure numpy-python version. Something to look into.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAP7f1AjcHrsopXjwK5uYdALeSrokxLMwA7xebTikHyhwL-%2BOVg%40mail.gmail.com.

[Jason Moore]

Jim and others,

Here are the benchmarks I made yesterday:

http://www.moorepants.info/blog/fast-matrix-eval.html

The working code is here: https://gist.github.com/moorepants/6ef8ab450252789a1411

Any feedback is welcome.

Jason
moorepants.info
+01 530-601-9791

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAJ2L7me73iJmkWm%3D_LiyWrsuOCZm%2B4OZbqD%2BkwwScWWx23HVdg%40mail.gmail.com.

[Tim Lahey]

On why Fortran is faster, Fortran semantics ensure that function arguments never alias, this allows the optimizer to make assumptions about the function and the arguments. This the main advantage of Fortran over C. But, because of this, it can lead to more memory usage. I know that the newer C++ standards have a keyword to mark arguments to indicate that they won't be aliased, but that requires that the code generator and the compiler support them.

Cheers,

Tim.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAP7f1Agdi_X-o0B%2B9mH2CGOSN-TyYGVwgZm4q8%3DYwxieBzZkzA%40mail.gmail.com.

[James Crist]

I was planning on going to bed, but ended up working on this instead. I have no self control...

Anyway, I've uncovered some things:

1. Addition of the restrict keyword to tell the compiler we're not aliasing offers marginal gains. Gain a couple microseconds here and there. This requires a c99 compiler, but it's 2014, everyone should have one by now.

2. Inlining the function call resulted in smaller gains than 1, but still *slightly* measurable. I suspect that for larger expression sizes this will be negligible to none.

3. Here's the big one: For small powers, pow(c, n) is considerably slower than c*c*c*c... Changing the ccode Pow handler to print all pows less than 5 (arbitrary number) out as multiplication I was able to match/beat (slightly) all of jason's benchmarks with the C + numpy ufuncs.

[Jason Moore]

Jason
moorepants.info
+01 530-601-9791

On Fri, Aug 29, 2014 at 2:38 AM, James Crist <cris...@umn.edu> wrote:

I was planning on going to bed, but ended up working on this instead. I have no self control...

Anyway, I've uncovered some things:

1. Addition of the restrict keyword to tell the compiler we're not aliasing offers marginal gains. Gain a couple microseconds here and there. This requires a c99 compiler, but it's 2014, everyone should have one by now.

2. Inlining the function call resulted in smaller gains than 1, but still *slightly* measurable. I suspect that for larger expression sizes this will be negligible to none.

3. Here's the big one: For small powers, pow(c, n) is considerably slower than c*c*c*c... Changing the ccode Pow handler to print all pows less than 5 (arbitrary number) out as multiplication I was able to match/beat (slightly) all of jason's benchmarks with the C + numpy ufuncs.

Oh yes! I knew that. In fact, I feel like I read in the current code somewhere. I forget, but that seems like a standard way we should be handling pows in C. Nice find!

 

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/6cfe63df-df00-4c36-a88a-6c477becc924%40googlegroups.com.

[Jason Moore]

Here is some work on the pow issue: https://github.com/sympy/sympy/pull/7519

Looks like it was merged so the ccode printer should print x*x*x... for less that 10 x's.

Jason
moorepants.info
+01 530-601-9791

[Jason Moore]

Sorry, it wasn't merged. He found that the --fast-math flag in the complier takes care of this.

Jason
moorepants.info
+01 530-601-9791

[Aaron Meurer]

I think we should print pow using repeated multiplication. People
might not know about --ffast-math, not realize that we are using pow
and that it is needed, or not want other optimizations that it
provides.

Is there a reason to put a limit on the power (5 was suggested here,
10 on the pull request)?

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CAP7f1AhNvQAwJ3V8y7uvSb2nDTpKDd8u8eiKVmjkOT-JZX4S2w%40mail.gmail.com.

[James Crist]

If I understand correctly, there is no cost in representing pow(x, n) as x*x*x*x... for any positive integer n, as long as it's done correctly. C compilers don't like to change how you write out calculations unless they're asked too. So x*x*x*x will not generate the same machine code as (x*x)*(x*x). The second case is preferable, as it will result in y = x*x, sol = y*y, rather than sol = x*x*x*x. This removes the need for one computation. Also, the second way apparently results in better precision for the end result (not sure why).

I am all for writing all positive integer powers as multiplication, provided we can get the parenthesis convention correct. So x**4 -> (x*x)*(x*x), or x**11 -> x*((x*x)*(x*x)*x)*((x*x)*(x*x)*x), etc... If others are supportive of this, I'll submit a PR.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6KdCQso6Nobdfeduyu395n-R5VjutEiGiZesTLL17siLA%40mail.gmail.com.

[Tim Lahey]

I recommend that you use the horner function in polys.

>> https://groups.google.com/d/msgid/sympy/CAJ2L7mfL_xO%3DO-ZRMx-zfpZzJKJ-%2BUdTzSCz5jYf%2B%3DdovR%2B_7Q%40mail.gmail.com
>> .

>>>>>>>>
>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> You received this message because you are subscribed to a topic
>>>>>>>> in
>>>>>>>> the Google Groups "sympy" group.
>>>>>>>> To unsubscribe from this topic, visit
>>>>>>>> https://groups.google.com/d/topic/sympy/azVZHLOv9Vc/unsubscribe.
>>>>>>>> To unsubscribe from this group and all its topics, send an
>>>>>>>> email to
>>>>>>>> sympy+un...@googlegroups.com.
>>>>>>>> To post to this group, send email to sy...@googlegroups.com.
>>>>>>>> Visit this group at http://groups.google.com/group/sympy.
>>>>>>>> To view this discussion on the web visit
>>>>>>>>

>> https://groups.google.com/d/msgid/sympy/CAP7f1AjcHrsopXjwK5uYdALeSrokxLMwA7xebTikHyhwL-%2BOVg%40mail.gmail.com
>> .

>>>>>>>>
>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> You received this message because you are subscribed to the
>>>>>>>> Google
>>>>>>>> Groups "sympy" group.
>>>>>>>> To unsubscribe from this group and stop receiving emails from
>>>>>>>> it,
>>>>>>>> send an email to sympy+un...@googlegroups.com.
>>>>>>>> To post to this group, send email to sy...@googlegroups.com.
>>>>>>>> Visit this group at http://groups.google.com/group/sympy.
>>>>>>>> To view this discussion on the web visit
>>>>>>>>

>> https://groups.google.com/d/msgid/sympy/CAJ2L7me73iJmkWm%3D_LiyWrsuOCZm%2B4OZbqD%2BkwwScWWx23HVdg%40mail.gmail.com
>> .

>>>>>>>>
>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> You received this message because you are subscribed to the
>>>>>>>> Google
>>>>>>>> Groups "sympy" group.
>>>>>>>> To unsubscribe from this group and stop receiving emails from
>>>>>>>> it,
>>>>>>>> send an email to sympy+un...@googlegroups.com.
>>>>>>>> To post to this group, send email to sy...@googlegroups.com.
>>>>>>>> Visit this group at http://groups.google.com/group/sympy.
>>>>>>>> To view this discussion on the web visit
>>>>>>>>

>> https://groups.google.com/d/msgid/sympy/CAP7f1Agdi_X-o0B%2B9mH2CGOSN-TyYGVwgZm4q8%3DYwxieBzZkzA%40mail.gmail.com
>> .

>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>
>>>>>> --
>>>>>> You received this message because you are subscribed to the
>>>>>> Google
>>>>>> Groups "sympy" group.
>>>>>> To unsubscribe from this group and stop receiving emails from it,
>>>>>> send
>>>>>> an email to sympy+un...@googlegroups.com.
>>>>>> To post to this group, send email to sy...@googlegroups.com.
>>>>>> Visit this group at http://groups.google.com/group/sympy.
>>>>>> To view this discussion on the web visit
>>>>>>

>> https://groups.google.com/d/msgid/sympy/6cfe63df-df00-4c36-a88a-6c477becc924%40googlegroups.com
>> .

>>>>>>
>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>
>>>>>
>>>>
>>>
>>> --
>>> You received this message because you are subscribed to the Google
>>> Groups
>>> "sympy" group.
>>> To unsubscribe from this group and stop receiving emails from it,
>>> send an
>>> email to sympy+un...@googlegroups.com.
>>> To post to this group, send email to sy...@googlegroups.com.
>>> Visit this group at http://groups.google.com/group/sympy.
>>> To view this discussion on the web visit
>>>

>> https://groups.google.com/d/msgid/sympy/CAP7f1AhNvQAwJ3V8y7uvSb2nDTpKDd8u8eiKVmjkOT-JZX4S2w%40mail.gmail.com
>> .

>>>
>>> For more options, visit https://groups.google.com/d/optout.
>>
>> --
>> You received this message because you are subscribed to a topic in
>> the
>> Google Groups "sympy" group.
>> To unsubscribe from this topic, visit
>> https://groups.google.com/d/topic/sympy/azVZHLOv9Vc/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to
>> sympy+un...@googlegroups.com.
>> To post to this group, send email to sy...@googlegroups.com.
>> Visit this group at http://groups.google.com/group/sympy.
>> To view this discussion on the web visit

>> https://groups.google.com/d/msgid/sympy/CAKgW%3D6KdCQso6Nobdfeduyu395n-R5VjutEiGiZesTLL17siLA%40mail.gmail.com
>> .

>> For more options, visit https://groups.google.com/d/optout.
>>
>
> --
> You received this message because you are subscribed to the Google
> Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/sympy/CAJ2L7mf3C0fB8Y5ynUxDui%3DOZBVyVPGdwK_NEWJqyfKzGSwizw%40mail.gmail.com.

[James Crist]

For handling Pow? horner(x**11) results in x**11. Or were you recommending applying horner to an entire expression tree?

.

an email to sympy+unsubscribe@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit

https://groups.google.com/d/msgid/sympy/6cfe63df-df00-4c36-a88a-6c477becc924%40googlegroups.com
.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an

email to sympy+unsubscribe@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit

https://groups.google.com/d/msgid/sympy/CAP7f1AhNvQAwJ3V8y7uvSb2nDTpKDd8u8eiKVmjkOT-JZX4S2w%40mail.gmail.com

.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the
Google Groups "sympy" group.
To unsubscribe from this topic, visit
https://groups.google.com/d/topic/sympy/azVZHLOv9Vc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to

sympy+unsubscribe@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/CAKgW%3D6KdCQso6Nobdfeduyu395n-R5VjutEiGiZesTLL17siLA%40mail.gmail.com
.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "sympy" group.

To unsubscribe from this group and stop receiving emails from it, send an email to sympy+unsubscribe@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAJ2L7mf3C0fB8Y5ynUxDui%3DOZBVyVPGdwK_NEWJqyfKzGSwizw%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/azVZHLOv9Vc/unsubscribe.

To unsubscribe from this group and all its topics, send an email to sympy+unsubscribe@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/9F9BC8F8-9C20-4623-A880-BA3E1DACE116%40gmail.com.

[Tim Lahey]

Try applying it to the expression tree. Horner's rule is supposed to be
optimal for evaluating polynomials.

>>>>>>>>>> https://groups.google.com/d/msgid/sympy/CAJ2L7mfL_xO%3DO-
>>>> ZRMx-zfpZzJKJ-%2BUdTzSCz5jYf%2B%3DdovR%2B_7Q%40mail.gmail.com

>>>>
>>>> .
>>>>
>>>>>
>>>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> You received this message because you are subscribed to a
>>>>>>>>>> topic in
>>>>>>>>>> the Google Groups "sympy" group.
>>>>>>>>>> To unsubscribe from this topic, visit
>>>>>>>>>> https://groups.google.com/d/topic/sympy/azVZHLOv9Vc/unsubscribe.
>>>>>>>>>> To unsubscribe from this group and all its topics, send an
>>>>>>>>>> email to
>>>>>>>>>> sympy+un...@googlegroups.com.
>>>>>>>>>> To post to this group, send email to sy...@googlegroups.com.
>>>>>>>>>> Visit this group at http://groups.google.com/group/sympy.
>>>>>>>>>> To view this discussion on the web visit
>>>>>>>>>>
>>>>>>>>>> https://groups.google.com/d/msgid/sympy/

>>>> CAP7f1AjcHrsopXjwK5uYdALeSrokxLMwA7xebTikHyhwL-%2BOVg%40mail.gmail.com

>>>>
>>>> .
>>>>
>>>>>
>>>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> You received this message because you are subscribed to the
>>>>>>>>>> Google
>>>>>>>>>> Groups "sympy" group.
>>>>>>>>>> To unsubscribe from this group and stop receiving emails from
>>>>>>>>>> it,
>>>>>>>>>> send an email to sympy+un...@googlegroups.com.
>>>>>>>>>> To post to this group, send email to sy...@googlegroups.com.
>>>>>>>>>> Visit this group at http://groups.google.com/group/sympy.
>>>>>>>>>> To view this discussion on the web visit
>>>>>>>>>>

>>>>>>>>>> https://groups.google.com/d/msgid/sympy/CAJ2L7me73iJmkWm%
>>>> 3D_LiyWrsuOCZm%2B4OZbqD%2BkwwScWWx23HVdg%40mail.gmail.com

>>>>
>>>> .
>>>>
>>>>>
>>>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> You received this message because you are subscribed to the
>>>>>>>>>> Google
>>>>>>>>>> Groups "sympy" group.
>>>>>>>>>> To unsubscribe from this group and stop receiving emails from
>>>>>>>>>> it,
>>>>>>>>>> send an email to sympy+un...@googlegroups.com.
>>>>>>>>>> To post to this group, send email to sy...@googlegroups.com.
>>>>>>>>>> Visit this group at http://groups.google.com/group/sympy.
>>>>>>>>>> To view this discussion on the web visit
>>>>>>>>>>

>>>>>>>>>> https://groups.google.com/d/msgid/sympy/CAP7f1Agdi_X-o0B%
>>>> 2B9mH2CGOSN-TyYGVwgZm4q8%3DYwxieBzZkzA%40mail.gmail.com

>>>>
>>>> .
>>>>
>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> --
>>>>>>>> You received this message because you are subscribed to the
>>>>>>>> Google
>>>>>>>> Groups "sympy" group.
>>>>>>>> To unsubscribe from this group and stop receiving emails from
>>>>>>>> it,
>>>>>>>> send
>>>>>>>> an email to sympy+un...@googlegroups.com.
>>>>>>>> To post to this group, send email to sy...@googlegroups.com.
>>>>>>>> Visit this group at http://groups.google.com/group/sympy.
>>>>>>>> To view this discussion on the web visit
>>>>>>>>

>>>>>>>> https://groups.google.com/d/msgid/sympy/6cfe63df-df00-
>>>> 4c36-a88a-6c477becc924%40googlegroups.com

>>>> .
>>>>
>>>>>
>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>> --
>>>>> You received this message because you are subscribed to the Google
>>>>> Groups
>>>>> "sympy" group.
>>>>> To unsubscribe from this group and stop receiving emails from it,
>>>>> send
>>>>> an
>>>>> email to sympy+un...@googlegroups.com.
>>>>> To post to this group, send email to sy...@googlegroups.com.
>>>>> Visit this group at http://groups.google.com/group/sympy.
>>>>> To view this discussion on the web visit
>>>>>

>>>>> https://groups.google.com/d/msgid/sympy/CAP7f1AhNvQAwJ3V8y7uvSb2nDTpKD
>>>> d8u8eiKVmjkOT-JZX4S2w%40mail.gmail.com

>>>>
>>>> .
>>>>
>>>>>
>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>
>>>>
>>>> --
>>>> You received this message because you are subscribed to a topic in
>>>> the
>>>> Google Groups "sympy" group.
>>>> To unsubscribe from this topic, visit
>>>> https://groups.google.com/d/topic/sympy/azVZHLOv9Vc/unsubscribe.
>>>> To unsubscribe from this group and all its topics, send an email to
>>>> sympy+un...@googlegroups.com.
>>>> To post to this group, send email to sy...@googlegroups.com.
>>>> Visit this group at http://groups.google.com/group/sympy.
>>>> To view this discussion on the web visit

>>>> https://groups.google.com/d/msgid/sympy/CAKgW%3D6KdCQso6Nobdfeduyu395n-
>>>> R5VjutEiGiZesTLL17siLA%40mail.gmail.com

>>>> .
>>>> For more options, visit https://groups.google.com/d/optout.
>>>>
>>>>
>>> --
>>> You received this message because you are subscribed to the Google
>>> Groups
>>> "sympy" group.
>>> To unsubscribe from this group and stop receiving emails from it,
>>> send an
>>> email to sympy+un...@googlegroups.com.
>>> To post to this group, send email to sy...@googlegroups.com.
>>> Visit this group at http://groups.google.com/group/sympy.
>>> To view this discussion on the web visit
>>> https://groups.google.com/d/

>>> msgid/sympy/CAJ2L7mf3C0fB8Y5ynUxDui%3DOZBVyVPGdwK_NEWJqyfKzGSwizw%
>>> 40mail.gmail.com.

>>>
>>> For more options, visit https://groups.google.com/d/optout.
>>>
>>
>> --
>> You received this message because you are subscribed to a topic in
>> the
>> Google Groups "sympy" group.
>> To unsubscribe from this topic, visit https://groups.google.com/d/
>> topic/sympy/azVZHLOv9Vc/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to
>> sympy+un...@googlegroups.com.
>> To post to this group, send email to sy...@googlegroups.com.
>> Visit this group at http://groups.google.com/group/sympy.
>> To view this discussion on the web visit https://groups.google.com/d/

>> msgid/sympy/9F9BC8F8-9C20-4623-A880-BA3E1DACE116%40gmail.com.

>>
>> For more options, visit https://groups.google.com/d/optout.
>>
>
> --
> You received this message because you are subscribed to the Google
> Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/sympy/CAJ2L7meCyqxdLgi7cZ%3D4pnbUPstwne5HwVtAvttY3mVGhzppCw%40mail.gmail.com.

[Øyvind Jensen]

James,
this contribution is just awesome! Great work!

Assuming that no serious issues show up, I'd vote that from 0.7.6, ufuncify should use your new backend by default.

Øyvind

[Björn Dahlgren]

Cool! Great work!

On Friday, 29 August 2014 22:48:23 UTC+2, James Crist wrote:

If I understand correctly, there is no cost in representing pow(x, n) as x*x*x*x... for any positive integer n, as long as it's done correctly.

Compilation will be slower for large n.
 

C compilers don't like to change how you write out calculations unless they're asked too. So x*x*x*x will not generate the same machine code as (x*x)*(x*x). The second case is preferable, as it will result in y = x*x, sol = y*y, rather than sol = x*x*x*x. This removes the need for one computation. Also, the second way apparently results in better precision for the end result (not sure why).

I am all for writing all positive integer powers as multiplication, provided we can get the parenthesis convention correct. So x**4 -> (x*x)*(x*x), or x**11 -> x*((x*x)*(x*x)*x)*((x*x)*(x*x)*x), etc... If others are supportive of this, I'll submit a PR.

I prefer either to pass -ffast-math flag (setting compiler flags is already an issue since we need
to indicate optimization level, right?) or write a specialized callback to be inlined (tuned to an Intel i7):
https://github.com/bjodah/cinterpol/blob/master/cInterpol/power.c

and tell sympy to generate code using it:
https://github.com/sympy/sympy/pull/7517

[James Crist]

I prefer either to pass -ffast-math flag (setting compiler flags is already an issue since we need
to indicate optimization level, right?) or write a specialized callback to be inlined (tuned to an Intel i7):

I'm not a fan of using the -ffast-math flag. It does other things that may result in a loss of precision, and not all users may want to use it. Handling compiler options should definitely be added to the codewrappers though, but should be a separate issue.

I do like the specialized function idea, as it can easily be done using our current framework. Instead of linking to pow in math.h by default for integer exponents, we can use a specialized integer exponent function that is printed by the codeprinter (separate header file? In the header file? In the .c file? Not sure the best spot). This way the user can override our default as they normally would, and requires no additional kwargs. If pow is overridden by the user, it will not be included in the code printout. Easy, simple, and requires no user facing changes.

Thoughts?

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/8ad9fbf6-144b-44e7-a55c-47e065ee6853%40googlegroups.com.

[Frédéric Bastien]

Using explicit multiplicatin (even with the parenthesis) lower the precission of the result compared to pow.

That being said, we do that in Theano. Putting a limit on the exponent that you do that would allow to don't loose too much precission for big exponent.

Fred

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAJ2L7mf3C0fB8Y5ynUxDui%3DOZBVyVPGdwK_NEWJqyfKzGSwizw%40mail.gmail.com.

[Björn Dahlgren]

 

That being said, we do that in Theano. Putting a limit on the exponent that you do that would allow to don't loose too much precission for big exponent.

Yes, that sounds like a good idea. Did you make any investigation on what is a suitable limit?
 

I do like the specialized function idea, as it can easily be done using our current framework. Instead of linking to pow in math.h by default for integer exponents, we can use a specialized integer exponent function that is printed by the codeprinter (separate header file? In the header file? In the .c file? Not sure the best spot). This way the user can override our default as they normally would, and requires no additional kwargs. If pow is overridden by the user, it will not be included in the code printout. Easy, simple, and requires no user facing changes.

I agree that it's a good idea. I guess a single header which is included is easiest to start with?
Should the functions be implemented in C? or C and Fortran? Or C++ even?
Are there any other functions except pow that needs special implementation?

[Ondřej Čertík]

Hi James,

On Mon, Sep 1, 2014 at 2:37 PM, James Crist <cris...@umn.edu> wrote:
>>
>> I prefer either to pass -ffast-math flag (setting compiler flags is already an issue since we need
>> to indicate optimization level, right?) or write a specialized callback to be inlined (tuned to an Intel i7):
>
>
> I'm not a fan of using the -ffast-math flag. It does other things that may result in a loss of precision, and not all users may want to use it. Handling compiler options should definitely be added to the codewrappers though, but should be a separate issue.

Can you point me to some examples from your experience when
-ffast-math may result in loss of precision? I use it all the time in
all my production codes. So I would like to learn more about the
pitfalls.

Ondrej

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAJ2L7mf2YbCPXcrkBtdOQXRrOBiM_twUyd47Ne-TboECGX%3DFew%40mail.gmail.com.

[James Crist]

I agree that it's a good idea. I guess a single header which is included is easiest to start with?
Should the functions be implemented in C? or C and Fortran? Or C++ even?
Are there any other functions except pow that needs special implementation?

As we're generating C compatible code, they should be in C. There is no need for a Fortran variant; Fortan already handles Pow(double, int) in an efficient way. I can't think of any other functions that should be implemented currently. Adding to the already created header sounds fine, although it may be unexpected. Right now I'm kind of leaning to the creation and inclusion of a separate header-only library that includes SymPy's specially defined functions.

Can you point me to some examples from your experience when
-ffast-math may result in loss of precision? I use it all the time in
all my production codes. So I would like to learn more about the
pitfalls.

Personally, I've never had an issue arise from using it, but a number of people more knowledgable than myself have warned against it. -ffast-math is just an easy way of setting a number of compiler flags that may result in more optimized math operations. Most of these result in code that is functional, and does what is intended, but no longer follows the IEEE floating point standard exactly. As the GCC manual says:

This option should never be turned on by any -O option since it can result in incorrect output for programs which depend on an exact implementation of IEEE or ISO rules/specifications for math functions. 

Looking through the flags that -ffast-math sets, it seems that -funsafe-math-operations is the big one that may cause issues. This allows a larger set of transformations on the computations to be used, which may result in changes to summation error, over/underflow behavior, and differences in rounding behavior from the code *as written*. For many things these may not matter, but for some they may.

- Jim

 

[Ondřej Čertík]

On Wed, Sep 3, 2014 at 1:21 PM, James Crist <cris...@umn.edu> wrote:
> Can you point me to some examples from your experience when
>> -ffast-math may result in loss of precision? I use it all the time in
>> all my production codes. So I would like to learn more about the
>> pitfalls.
>
>
> Personally, I've never had an issue arise from using it, but a number of
> people more knowledgable than myself have warned against it. -ffast-math is
> just an easy way of setting a number of compiler flags that may result in
> more optimized math operations. Most of these result in code that is
> functional, and does what is intended, but no longer follows the IEEE
> floating point standard exactly. As the GCC manual says:
>
>> This option should never be turned on by any -O option since it can result
>> in incorrect output for programs which depend on an exact implementation of
>> IEEE or ISO rules/specifications for math functions.
>
>
> Looking through the flags that -ffast-math sets, it seems that
> -funsafe-math-operations is the big one that may cause issues. This allows a
> larger set of transformations on the computations to be used, which may
> result in changes to summation error, over/underflow behavior, and
> differences in rounding behavior from the code *as written*. For many things
> these may not matter, but for some they may.

Yes, I am aware of this. Note that the Fortran language gives
compilers some freedom anyway in rearranging your expression when
evaluating.

In practice, I never had any issue with -ffast-math, because a good
algorithm in my opinion shouldn't really depend on details of IEEE
implementation. I could be wrong of course, that's why I was curious
if you actually have an example.

Ondrej

[Björn Dahlgren]

On Monday, 1 September 2014 22:37:31 UTC+2, James Crist wrote:

/.../ Handling compiler options should definitely be added to the codewrappers though, but should be a separate issue.

Yes, I even went so far to write my own package to handle compilaton (github.com/bjodah/pycompilation). Looking back
maybe I should have turned to WAF or something similar instead of reinventing the wheel (an inferior wheel maybe even).

When (if) we tackle this, it would be great if we allow sympy to cache e.g. object files from compiled cython wrapper
code - at least for me cythonizing/compiling the cython generated C file has been a bottle neck more than once.

On Wednesday, 3 September 2014 21:21:06 UTC+2, James Crist wrote:

Adding to the already created header sounds fine, although it may be unexpected. Right now I'm kind of leaning to the creation and inclusion of a separate header-only library that includes SymPy's specially defined functions.

header only library sounds fine to me, though if there is only one single function "sympy_power_double_int" or
whatever it will be named it may be overkill.