Piecewise and Codegen

[Nathan Woods]

 Hi.

I want to generate some C code that evaluates a discontinuous, symbolic, function for a given set of arguments. codegen works great, except for the discontinuous bit. I tried defining the function as a Piecewise, but no luck. Is there a better way to do this? Is this an easy fix? The best idea I've had so far is to use Heaviside functions for the discontinuities and write my own Heaviside implementation for inclusion by codegen.

Thanks for any help.

Nathan Woods

P.S. Thanks for a great product!

[Ondřej Čertík]

Hi Nathan,

On Mon, Oct 7, 2013 at 8:47 PM, Nathan Woods <charle...@gmail.com> wrote:
> Hi.
>
> I want to generate some C code that evaluates a discontinuous, symbolic,
> function for a given set of arguments. codegen works great, except for the

Can you post the function that you are trying to convert to a C code?

> discontinuous bit. I tried defining the function as a Piecewise, but no

That seems like the most natural way to do it. Can you post the error message?
Let's fix it.

> luck. Is there a better way to do this? Is this an easy fix? The best idea
> I've had so far is to use Heaviside functions for the discontinuities and
> write my own Heaviside implementation for inclusion by codegen.
>
> Thanks for any help.
>
> Nathan Woods
>
> P.S. Thanks for a great product!

I am glad it's useful!

Ondrej

[Nathan Woods]

The sympy function (the output of printing the expression) is:
"t*x*y + x**2 + y**2 + Piecewise((0, x < 0.5), (1, x >= 0.5)) + cos(t) - 1"

The error message is:

TypeError: unsupported operand type(s) for |=: 'set' and 'Piecewise'

The whole error output is:

Traceback (most recent call last):
  File "integration.py", line 201, in <module>
    print IntegrableFunction(integrands[0],**test).integrate()
  File "integration.py", line 52, in __init__
    sympy_function,self.sympy_variables,args=self.args)
  File "integration.py", line 167, in __init__
    argument_sequence=self.sympy_variables,to_files=True)
  File "/Users/woodscn/SCIPYTEST/lib/python2.7/site-packages/sympy/utilities/codegen.py", line 958, in codegen
    return code_gen.write(routines, prefix, to_files, header, empty)
  File "/Users/woodscn/SCIPYTEST/lib/python2.7/site-packages/sympy/utilities/codegen.py", line 445, in write
    dump_fn(self, routines, f, prefix, header, empty)
  File "/Users/woodscn/SCIPYTEST/lib/python2.7/site-packages/sympy/utilities/codegen.py", line 615, in dump_c
    self.dump_code(routines, f, prefix, header, empty)
  File "/Users/woodscn/SCIPYTEST/lib/python2.7/site-packages/sympy/utilities/codegen.py", line 489, in dump_code
    code_lines.extend(self._call_printer(routine))
  File "/Users/woodscn/SCIPYTEST/lib/python2.7/site-packages/sympy/utilities/codegen.py", line 593, in _call_printer
    constants, not_c, c_expr = ccode(result.expr, assign_to=assign_to, human=False)
  File "/Users/woodscn/SCIPYTEST/lib/python2.7/site-packages/sympy/printing/ccode.py", line 262, in ccode
    return CCodePrinter(settings).doprint(expr, assign_to)
  File "/Users/woodscn/SCIPYTEST/lib/python2.7/site-packages/sympy/printing/ccode.py", line 87, in doprint
    code0 = self._doprint_a_piece(expr, assign_to)
  File "/Users/woodscn/SCIPYTEST/lib/python2.7/site-packages/sympy/printing/codeprinter.py", line 30, in _doprint_a_piece
    d = get_contraction_structure(expr)
  File "/Users/woodscn/SCIPYTEST/lib/python2.7/site-packages/sympy/tensor/index_methods.py", line 401, in get_contraction_structure
    result[key] |= d[key]
TypeError: unsupported operand type(s) for |=: 'set' and 'Piecewise'

Thanks for the quick response.

--
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.
For more options, visit https://groups.google.com/groups/opt_out.

[Ondřej Čertík]

Nathan,

I see, here is how to reproduce it: https://gist.github.com/certik/6880300

This is definitely a bug. Franz, any idea why it would call the tensor module?

Ondrej

[F. B.]

/sympy/tensor/index_methods.pyc

I have never worked on this file. The real tensor module is /sympy/tensor/tensor.py, the other files in that directory are tools for code generation (I think), which are able to handle indexed data, but I never used them.

[Nathan Woods]

I was wondering if there were some kind of estimated time horizon on this. Is it just a bug fix with a few hours if work, or is it a new feature that would have to be added, potentially taking a few weeks? I'll need to work around the problem, if that's the case. 

N

—
Sent from Mailbox for iPhone

--

[Ondřej Čertík]

Nathan,

I am working on a fix.

Can you provide the exact C expression that you want Piecewise((0, x <
0.5), (1, x >= 0.5)) to generate?
I.e. can you fill in the right hand side here:

assert ccode(Piecewise((0, x < 0.5), (1, x >= 0.5))) == "piecewise(...)"

Ondrej

[Matthew Rocklin]

I could take a look at this problem but probably won't get to it for a day or two.  My understanding was that ccodegen supported piecewise though.  Can anyone identify what the issue is?  Looking at Ondrej's error printout my guess is that it's a simple issue. 

I more actively maintain the SymPy->Theano conversion which can provide efficient Theano functions that run at C speed but don't give nice code to use in external applications (so this may or may not be of use to you).

In general I respond to errors on that system within a day or so.

Some documentation on this is here
https://github.com/mrocklin/sympy/blob/e610a6898b67d43eb7adc5b014007225aeb1eb90/doc/src/modules/numeric-computation.rst
https://github.com/sympy/sympy/pull/2400

Halfway through writing this I see that Ondrej is working on it.  Great!

[Ondřej Čertík]

On Tue, Oct 8, 2013 at 1:58 PM, Ondřej Čertík <ondrej...@gmail.com> wrote:
> Nathan,
>
> I am working on a fix.
>
> Can you provide the exact C expression that you want Piecewise((0, x <
> 0.5), (1, x >= 0.5)) to generate?
> I.e. can you fill in the right hand side here:
>
> assert ccode(Piecewise((0, x < 0.5), (1, x >= 0.5))) == "piecewise(...)"

Ok, it actually works... It produces:

if (x < 0.5) {
0
}
else if (x >= 0.5) {
1
}

But obviously this will not work inside an expression. Can you provide
us the exact C code that you expect
for your expression?

Ondrej

[Nathan Woods]

I would be happy with either of the following implementations, one or the other of which might be preferred for other reasons. The immediate intended use is to wrap the resulting function in ctypes so that I can feed it to some existing code.

- An if/then construct, like what you mentioned. I don't quite understand why that doesn't work, though.

- 0 + (1-0)*(x>=0.5) (essentially a Heaviside implementation. C interprets a "false" as 0, and "true" as 1, correct?)

Ondrej

[Matthew Rocklin]

 A simple fix would be to put the entire expression within the Piecewise.  E.g.

Instead of

x + Piecewise((0, x < 0), (1, x > 0))

Try

Piecewise((x, x < 0), (x + 1, x > 0))

You shouldn't have to do this (ccodegen should be smart enough to handle this), but I suspect it will work in the short term.

[Ondřej Čertík]

Thanks Matthew!

That's right. Nathan, let me know if this works for you:

In [1]: e = t*x*y + x**2 + y**2 + Piecewise((0, x < 0.5), (1, x >=
0.5)) + cos(t) - 1

In [2]: ccode(piecewise_fold(e))
Out[2]: 'if (x < 0.5) {\n t*x*y + pow(x, 2) + pow(y, 2) + cos(t) -
1\n}\nelse if (x >= 0.5) {\n t*x*y + pow(x, 2) + pow(y, 2) +
cos(t)\n}'


Ondrej

[Nathan Woods]

That will actually work fine for my application, and it also more closely approximates how piecewise functions are used in mathematics, which is a plus. Let me test it out.

[Nathan Woods]

I'm getting a compiler error using "gcc -c integrand.c".

integrand.c: In function ‘integrand’:
integrand.c:13:11: error: expected expression before ‘if’

The generated C-code is :

/******************************************************************************
 *                      Code generated with sympy 0.7.2                       *
 *                                                                            *
 *              See http://www.sympy.org/ for more information.               *
 *                                                                            *
 *                       This file is part of 'project'                       *
 ******************************************************************************/
#include "integrand.h"
#include <math.h>

double integrand(double t, double x, double y) {

   return if (x < 0.5) {

      t*x*y + pow(x, 2) + pow(y, 2) + cos(t) - 1
   }

   else if (x >= 0.5) {

      t*x*y + pow(x, 2) + pow(y, 2) + cos(t)

   };

}

[Nathan Woods]

On a longer-term note,

Will it be a problem to add another required ordering to Sympy? For instance, I've found that Piecewise cannot wrap a Matrix expression, and now Piecewise must be the outermost layer in an expression, forcing the ordering of Matrix(Piecewise(other_stuff)). Again, this reflects mathematical practice, so it's probably not a bad thing, but it doesn't jive quite so well with programming practice, where users might expect something that's a little less stringent.

[Ondřej Čertík]

On Tue, Oct 8, 2013 at 2:32 PM, Nathan Woods <charle...@gmail.com> wrote:
> I'm getting a compiler error using "gcc -c integrand.c".
>
> integrand.c: In function ‘integrand’:
> integrand.c:13:11: error: expected expression before ‘if’
>
> The generated C-code is :
>
> /******************************************************************************
> * Code generated with sympy 0.7.2
> *
> *
> *
> * See http://www.sympy.org/ for more information.
> *
> *
> *
> * This file is part of 'project'
> *
>
> ******************************************************************************/
> #include "integrand.h"
> #include <math.h>
>
> double integrand(double t, double x, double y) {
>
> return if (x < 0.5) {

^^^ I see, obviously the printing of Piecewise() is wrong.

>
> t*x*y + pow(x, 2) + pow(y, 2) + cos(t) - 1

^^^ this is wrong too.

> }
> else if (x >= 0.5) {
> t*x*y + pow(x, 2) + pow(y, 2) + cos(t)
> };
>
> }

I started fixing it here:

https://github.com/sympy/sympy/pull/2524



Ondrej

[Matthew Rocklin]

If you don't need human readable code the Theano solution should handle all of this just fine by the way.

Ondrej

[Nathan Woods]

Oh, and I should point out that I used a different function (whatever I had in my code) and that the returned values look correct to me.

Ondrej

[Nathan Woods]

My only real requirement is that it generate something that can be compiled and then wrapped with ctypes. I haven't checked out theano yet, but I'll take a look. I'm hoping to eventually be able to loosen the ctypes requirement too, since the real underlying requirement is to have a function pointer i can pass to fortran code in place of a python callback. Does this sound like something theano would help with?

[Ondřej Čertík]

On Tue, Oct 8, 2013 at 2:42 PM, Nathan Woods <charle...@gmail.com> wrote:
> My only real requirement is that it generate something that can be compiled
> and then wrapped with ctypes. I haven't checked out theano yet, but I'll
> take a look. I'm hoping to eventually be able to loosen the ctypes
> requirement too, since the real underlying requirement is to have a function
> pointer i can pass to fortran code in place of a python callback. Does this
> sound like something theano would help with?

It's almost fixed at:

https://github.com/sympy/sympy/pull/2524

See the comments there.

Ondrej

[Ondřej Čertík]

OK, see this comment:

https://github.com/sympy/sympy/pull/2524#issuecomment-25927909

It produces something that compiles. There are a few bugs around, but
the comment shows
how to workaround them.

I don't have more time to work on this. Nathan, this PR shows how to
fix these kind of issues.
If you could finish it up, so that we can merge it, that would be
tremendous help.

Just get my branch, push more commits and send a new PR.

Ondrej