variable limit in ufuncify?

92 views
Skip to first unread message

David Shin

unread,
Sep 25, 2014, 3:42:59 PM9/25/14
to sy...@googlegroups.com
Hi, I recently began trying out sympy and am running into some difficulty.

I wrote the following script, called demo.py:

import sympy
from sympy.utilities.autowrap import ufuncify
import sys

N = int(sys.argv[1])
theta = []
values = []
for n in range(N):
    theta.append(sympy.symbols('x%s' % n))
    values.append(n)

summation = sum(theta)
f = ufuncify(theta, summation)
print f(*values)[0]

Running it for small N, it works fine:

$ python demo.py 21
210.0

But it fails for larger N. Can anyone advise? Thanks in advance.

$ python demo.py 22
Traceback (most recent call last):
  File "demo.py", line 13, in 
    f = ufuncify(theta, summation)
  File "/opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/sympy/utilities/autowrap.py", line 485, in ufuncify
    return autowrap(C.Equality(y[i], f(*args)), **kwargs)
  File "/opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/sympy/utilities/autowrap.py", line 403, in autowrap
    return code_wrapper.wrap_code(routine, helpers=helps)
  File "/opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/sympy/utilities/autowrap.py", line 139, in wrap_code
    self._process_files(routine)
  File "/opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/sympy/utilities/autowrap.py", line 158, in _process_files
    " ".join(command), e.output))
sympy.utilities.autowrap.CodeWrapError: Error while executing command: f2py -m wrapper_module_0 -c wrapped_code_0.f90. Command output is:
running build
running config_cc
unifing config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
build_src
building extension "wrapper_module_0" sources
f2py options: []
f2py:> /tmp/tmpKbJQuO/src.linux-x86_64-2.7/wrapper_module_0module.c
creating /tmp/tmpKbJQuO
creating /tmp/tmpKbJQuO/src.linux-x86_64-2.7
Reading fortran codes...
        Reading file 'wrapped_code_0.f90' (format:free)
Post-processing...
        Block: wrapper_module_0
                        Block: autofunc
Post-processing (stage 2)...
Building modules...
        Building module "wrapper_module_0"...
                Constructing wrapper function "autofunc"...
                  y_15 = autofunc(x_16,x1,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x2,x20,x21,x3,x4,x5,x6,x7,x8,x9,[m_17])
        Wrote C/API module "wrapper_module_0" to file "/tmp/tmpKbJQuO/src.linux-x86_64-2.7/wrapper_module_0module.c"
  adding '/tmp/tmpKbJQuO/src.linux-x86_64-2.7/fortranobject.c' to sources.
  adding '/tmp/tmpKbJQuO/src.linux-x86_64-2.7' to include_dirs.
copying /opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpKbJQuO/src.linux-x86_64-2.7
copying /opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpKbJQuO/src.linux-x86_64-2.7
build_src: building npy-pkg config files
running build_ext
customize UnixCCompiler
customize UnixCCompiler using build_ext
customize Gnu95FCompiler
Found executable /opt/user/x86_64/gcc-4.7.2/bin/gfortran
customize Gnu95FCompiler
customize Gnu95FCompiler using build_ext
building 'wrapper_module_0' extension
compiling C sources
C compiler: gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC

creating /tmp/tmpKbJQuO/tmp
creating /tmp/tmpKbJQuO/tmp/tmpKbJQuO
creating /tmp/tmpKbJQuO/tmp/tmpKbJQuO/src.linux-x86_64-2.7
compile options: '-I/tmp/tmpKbJQuO/src.linux-x86_64-2.7 -I/opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include -I/opt/user/x86_64/Python-2.7.3/include/python2.7 -c'
gcc: /tmp/tmpKbJQuO/src.linux-x86_64-2.7/wrapper_module_0module.c
In file included from /opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1728:0,
                 from /opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17,
                 from /opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:15,
                 from /tmp/tmpKbJQuO/src.linux-x86_64-2.7/fortranobject.h:13,
                 from /tmp/tmpKbJQuO/src.linux-x86_64-2.7/wrapper_module_0module.c:18:
/opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include/numpy/npy_deprecated_api.h:11:2: warning: #warning "Using deprecated NumPy API, disable it by #defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
/tmp/tmpKbJQuO/src.linux-x86_64-2.7/wrapper_module_0module.c:111:12: warning: âpy_sizeâefined but not used [-Wunused-function]
gcc: /tmp/tmpKbJQuO/src.linux-x86_64-2.7/fortranobject.c
In file included from /opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1728:0,
                 from /opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17,
                 from /opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:15,
                 from /tmp/tmpKbJQuO/src.linux-x86_64-2.7/fortranobject.h:13,
                 from /tmp/tmpKbJQuO/src.linux-x86_64-2.7/fortranobject.c:2:
/opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include/numpy/npy_deprecated_api.h:11:2: warning: #warning "Using deprecated NumPy API, disable it by #defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
compiling Fortran sources
Fortran f77 compiler: /opt/user/x86_64/gcc-4.7.2/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran f90 compiler: /opt/user/x86_64/gcc-4.7.2/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran fix compiler: /opt/user/x86_64/gcc-4.7.2/bin/gfortran -Wall -ffixed-form -fno-second-underscore -Wall -fno-second-underscore -fPIC -O3 -funroll-loops
compile options: '-I/tmp/tmpKbJQuO/src.linux-x86_64-2.7 -I/opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include -I/opt/user/x86_64/Python-2.7.3/include/python2.7 -c'
gfortran:f90: wrapped_code_0.f90
wrapped_code_0.f90:1.133:

4, x15, x16, x17, x18, x19, x2, x20, x21, x3, x4, x5, x6, x7, x8, x9, y_15
                                                                           1
Warning: Line truncated at (1)
wrapped_code_0.f90:1.132:

14, x15, x16, x17, x18, x19, x2, x20, x21, x3, x4, x5, x6, x7, x8, x9, y_15
                                                                           1
Error: Unexpected junk in formal argument list at (1)
wrapped_code_0.f90:33.3:

end subroutine
   1
Error: Expecting END PROGRAM statement at (1)
Error: Unexpected end of file in 'wrapped_code_0.f90'
wrapped_code_0.f90:1.133:

4, x15, x16, x17, x18, x19, x2, x20, x21, x3, x4, x5, x6, x7, x8, x9, y_15
                                                                           1
Warning: Line truncated at (1)
wrapped_code_0.f90:1.132:

14, x15, x16, x17, x18, x19, x2, x20, x21, x3, x4, x5, x6, x7, x8, x9, y_15
                                                                           1
Error: Unexpected junk in formal argument list at (1)
wrapped_code_0.f90:33.3:

end subroutine
   1
Error: Expecting END PROGRAM statement at (1)
Error: Unexpected end of file in 'wrapped_code_0.f90'
error: Command "/opt/user/x86_64/gcc-4.7.2/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops -I/tmp/tmpKbJQuO/src.linux-x86_64-2.7 -I/opt/user/x86_64/Python-2.7.3/lib/python2.7/site-packages/numpy/core/include -I/opt/user/x86_64/Python-2.7.3/include/python2.7 -c -c wrapped_code_0.f90 -o /tmp/tmpKbJQuO/wrapped_code_0.o" failed with exit status 1

Jason Moore

unread,
Sep 25, 2014, 4:30:12 PM9/25/14
to sy...@googlegroups.com
David,

This is because it wasn't wrapping lines correctly in the generated Fortran code. If you use the development version of SymPy it should work.

Here is the PR that fixed it: https://github.com/sympy/sympy/pull/7968

--
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/68af9759-5988-4e85-81ec-7cc303022e46%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Jason Moore

unread,
Sep 25, 2014, 4:40:25 PM9/25/14
to sy...@googlegroups.com
It works for the backend='f2py' in master but fails for the default backend which is 'numpy'. I get a segmentation fault on the 'numpy' backend for values greater than 20 or so.

This is an odd use of ufuncify, as summing would better be done in a loop with an array as input to the function. You can probably used the IndexBased class to set up a loop based sum. If you use the tempdir kwarg to ufuncify you can see the code it generates, and you'll basically get a Fortran function that has as many input args as your integer value which is not very efficient.

David Shin

unread,
Sep 25, 2014, 4:52:10 PM9/25/14
to sy...@googlegroups.com
Would you mind showing me how to convert my demo.py to have ufuncify generate code that has a single length-N array argument instead of N separate arguments?

Jason Moore

unread,
Sep 25, 2014, 5:00:04 PM9/25/14
to sy...@googlegroups.com
First, why do you need sympy to do this? Would NumPy be sufficient?

import numpy as np
values = np.random.random(100)
np.sum(values)

David Shin

unread,
Sep 25, 2014, 5:08:02 PM9/25/14
to sy...@googlegroups.com
I'm looking to minimize a complicated function of hundreds of variables. I want to ufuncify() my function so that I can pass it into scipy.optimize.minimize(). The summation func I gave is just a toy example, the actual function is much more complicated.

Please let me know if I should be doing this some other way.

David Shin

unread,
Sep 25, 2014, 5:09:39 PM9/25/14
to sy...@googlegroups.com
Also, I'd like to pass in the gradient and hessian of my function into scipy.optimize.minimize() to help it along. I don't want to do the calculus by hand, so I'm using sympy to do it for me.

Jason Moore

unread,
Sep 25, 2014, 5:21:32 PM9/25/14
to sy...@googlegroups.com
I see. What does your function look like? Does it have summations or things you want to iterate over? Or is it simply a scalar function?

David Shin

unread,
Sep 25, 2014, 5:46:07 PM9/25/14
to sy...@googlegroups.com
The function is a sum of thousands of summands, where each summand is a small function of a few of the input parameters. These small summand functions are composed of basic operations like multiplication, addition, and exp().

Jason Moore

unread,
Sep 25, 2014, 6:31:27 PM9/25/14
to sy...@googlegroups.com
I'm not sure this is supported. Ideally you'd create a sympy.Sum object representing your summation and then the code printers would print a loop that would look something like:

>>> sum = sympy.Sum(a, (a, 1, 5))
>>> sympy.ccode(sum, b)
double b = 0;
for (i = 0; i < 5; i++;){
   b = b + a[i]
}

I think Jim Crist is working on this functionality but it doesn't exist now.

David Shin

unread,
Sep 25, 2014, 7:03:05 PM9/25/14
to sy...@googlegroups.com
Thanks for the info, Jason.

In my opinion, what I am trying to do is not particularly exotic. Consider for example the task of implementing logistic regression optimization via sympy and scipy. You have an n-by-k data matrix X of independent variables, a length-n output vector Y of dependent variables, and wish to approximate each Y[i] as a function of X[i], parameterized by a set of k+1 parameters, theta. This boils down to writing an objective loss function to minimize, which is a symbolic function of theta involving n summands that are individually expressions in terms Y[i], X[i], and theta. As the minimum does not have an analytic solution, you need to employ approximation algorithms, hence scipy.optimize.

Of course, there are numerous packages available for vanilla logistic regression. However, if you want to customize it (for example by using an alternative loss function or regularization penalty function), then you kind of need to roll out your own implementation.

From what you are saying, it sounds like sympy is not well-suited for this type of task.

Jason Moore

unread,
Sep 25, 2014, 7:55:31 PM9/25/14
to sy...@googlegroups.com
I'm not saying that SymPy isn't well suited for this task. I'm just saying that we don't have an implementation for code generation for summations that would help with your problem. You could certainly add one.

I use SymPy for optimization problems myself and generate the symbolic Jacobians and Hessians of the objective and constraint functions for non-linear programming problems and it works great. So SymPy can likely do what you want to do, but you may need to write a code printer for sympy.Sum or something similar.

James Crist

unread,
Sep 25, 2014, 8:52:12 PM9/25/14
to sy...@googlegroups.com
I'd hazard that it's a limit in the internals of numpy for how they handle broadcasting, but I can't be certain on that.

However, we can handle this, you just need to frame your problem in a better way. You're trying to do optimization, so generally you'd frame your optimal condition as a vector. In SymPy you should do the same (using a MatrixSymbol). You also don't need ufuncify, because you're not doing any broadcasting. Instead, you can use `autowrap` to create a function. ufuncs can't broadcast arrays, so ufuncify won't work for functions that require a vector as an input for a single iteration. The resulting expression will be just as fast as a ufunc for a single iteration though, so no harm done.

Here's your toy example, done as described:

import sympy
from sympy.utilities.autowrap import autowrap
import sys

N = int(sys.argv[1
])

x = sympy.MatrixSymbol('x', N, 1)
summation = sum(x)
values = list(range(N))

f
= autowrap(summation, args=(x,)) print f(values)

This works regardless of how large N is. For your real problem, you can formulate your expression using symbols, and use subs to form them into a single vector input. Or you can work directly with the vector by using elements of it as variables. It's up to you how to formulate it.

-Jim

Björn Dahlgren

unread,
Sep 26, 2014, 9:22:54 AM9/26/14
to sy...@googlegroups.com

David Shin

unread,
Sep 26, 2014, 11:10:11 AM9/26/14
to sy...@googlegroups.com
James,

Your code doesn't seem to work as expected:

$ python demo.py 10

Traceback (most recent call last):
  File "demo.py", line 11, in <module>
    f = autowrap(summation, args=(x,))
  File "/opt/jump/x86_64/Python-2.7.3/lib/python2.7/site-packages/sympy/utilities/autowrap.py", line 387, in autowrap
    routine = Routine('autofunc', expr, args)
  File "/opt/jump/x86_64/Python-2.7.3/lib/python2.7/site-packages/sympy/utilities/codegen.py", line 229, in __init__
    new_args.append(InputArgument(symbol))
  File "/opt/jump/x86_64/Python-2.7.3/lib/python2.7/site-packages/sympy/utilities/codegen.py", line 350, in __init__
    Variable.__init__(self, name, datatype, dimensions, precision)
  File "/opt/jump/x86_64/Python-2.7.3/lib/python2.7/site-packages/sympy/utilities/codegen.py", line 300, in __init__
    raise TypeError("The first argument must be a sympy symbol.")
TypeError: The first argument must be a sympy symbol.

It appears to me that sum(x) has type int.

James Crist

unread,
Sep 26, 2014, 11:12:23 AM9/26/14
to sy...@googlegroups.com
@david:

This is all in sympy master. The codegen stuff has had major work done to it in the last development cycle.

--
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/kgzcR8Qu7s4/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.

David Shin

unread,
Oct 1, 2014, 3:11:10 PM10/1/14
to sy...@googlegroups.com
Just to confirm, the latest official release is 0.7.5 and does not work with your code. Is that right? Do you know when the next release will be?

James Crist

unread,
Oct 1, 2014, 3:16:50 PM10/1/14
to sy...@googlegroups.com
Yes. The next release should be soon. We're just finishing up the last couple of milestones. Can't but a date on it yet for sure though. If you want to try it now though, you can pull the latest master off of github.

-Jim
Reply all
Reply to author
Forward
0 new messages