ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 13 Apr 2011 06:17:04 -0500Using Sage Symbolic Functions in Scipy fsolvehttp://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/I need to find the roots of a system of multivariate non-linear (algebraic) functions and Sage's ``solve`` function is running indefinitely. (Where "indefinitely" means "more than 30 mins".) I only require numerical solutions so I was hoping to use Scipy's ``fsolve`` function instead. However, the functions that I'm generating are Sage symbolic functions, which have been mightly convenient to work with, and I'm having trouble getting them into a format that ``fsolve`` will understand. Does anyone have any suggestions?
The documentation for Scipy's ``fsolve`` can be found [here](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve). It looks like a Python function is sufficient.
As a toy example, I've tried the following:
sage: from scipy.optimize import fsolve
sage: var('x,y')
sage: f = x^2 + y^2
sage: def ff(v):
return f(x=v[0],y=v[0])
sage: fsolve(ff,[1,1])
It seems like the naive approach but I'm receiving the error ``TypeError: no canonical coercion from <type 'numpy.int64'> to Symbolic Ring``. Perhaps this is an issue with Sage properly dealing with Numpy/Scipy data types?Sun, 10 Apr 2011 14:52:57 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/Answer by Kelvin Li for <p>I need to find the roots of a system of multivariate non-linear (algebraic) functions and Sage's <code>solve</code> function is running indefinitely. (Where "indefinitely" means "more than 30 mins".) I only require numerical solutions so I was hoping to use Scipy's <code>fsolve</code> function instead. However, the functions that I'm generating are Sage symbolic functions, which have been mightly convenient to work with, and I'm having trouble getting them into a format that <code>fsolve</code> will understand. Does anyone have any suggestions?</p>
<p>The documentation for Scipy's <code>fsolve</code> can be found <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve">here</a>. It looks like a Python function is sufficient.</p>
<p>As a toy example, I've tried the following:</p>
<pre><code>sage: from scipy.optimize import fsolve
sage: var('x,y')
sage: f = x^2 + y^2
sage: def ff(v):
return f(x=v[0],y=v[0])
sage: fsolve(ff,[1,1])
</code></pre>
<p>It seems like the naive approach but I'm receiving the error <code>TypeError: no canonical coercion from <type 'numpy.int64'> to Symbolic Ring</code>. Perhaps this is an issue with Sage properly dealing with Numpy/Scipy data types?</p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?answer=12284#post-id-12284The original scenario involves two variables and only one equation. Because `ff` accepts two variables but does not return a list (or tuple, or array, ...) of length 2, the "shape mismatch" error is thrown. Here is a working example (with two equations and minimal modifications):
import scipy.optimize
var('x,y')
f = x^2 + y^2 - 1
g = x - y^2
def ff(v):
x, y = map(float, v)
return [f(x=x, y=y), g(x=x, y=y)]
print scipy.optimize.fsolve(ff, [1, 1])
Here is a more compact version using lambda functions and the trick described link:[here](http://www.sagemath.org/doc/reference/sage/modules/vector_callable_symbolic_dense.html "test")
import scipy.optimize
var('x,y')
f(x,y) = (x^2 + y^2 - 1, x - y^2)
print scipy.optimize.fsolve(lambda v: f(*map(float, v)), (1, 1))
Since you said you have a large number of variables and equations, here is a slightly more "scale-friendly" approach (but somewhat kludgy):
import scipy.optimize
# Create symbolic variables "x_0" and "x_1".
# This is a pretty messy kludge
x = tuple(
SR.var(
"".join(['x_', str(t)])
)
for t in range(2)
)
# Represent system of equations using "x" as defined above.
# I use the intermediate "exprs" variable to allow breaking up the tuple
# over multiple lines.
exprs = (
x[0]^2 + x[1]^2 - 1,
x[0] - x[1]^2
)
f(*x) = exprs
# Guess values
guess = (1, 1)
print scipy.optimize.fsolve(lambda v: f(*map(float, v)), guess)
In all cases, the result should look like:
`[ 0.61803399 0.78615138]`
Sun, 10 Apr 2011 21:47:24 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?answer=12284#post-id-12284Comment by cswiercz for <p>The original scenario involves two variables and only one equation. Because <code>ff</code> accepts two variables but does not return a list (or tuple, or array, ...) of length 2, the "shape mismatch" error is thrown. Here is a working example (with two equations and minimal modifications):</p>
<pre><code>import scipy.optimize
var('x,y')
f = x^2 + y^2 - 1
g = x - y^2
def ff(v):
x, y = map(float, v)
return [f(x=x, y=y), g(x=x, y=y)]
print scipy.optimize.fsolve(ff, [1, 1])
</code></pre>
<p>Here is a more compact version using lambda functions and the trick described link:<a href="http://www.sagemath.org/doc/reference/sage/modules/vector_callable_symbolic_dense.html" title="test">here</a></p>
<pre><code>import scipy.optimize
var('x,y')
f(x,y) = (x^2 + y^2 - 1, x - y^2)
print scipy.optimize.fsolve(lambda v: f(*map(float, v)), (1, 1))
</code></pre>
<p>Since you said you have a large number of variables and equations, here is a slightly more "scale-friendly" approach (but somewhat kludgy):</p>
<pre><code>import scipy.optimize
# Create symbolic variables "x_0" and "x_1".
# This is a pretty messy kludge
x = tuple(
SR.var(
"".join(['x_', str(t)])
)
for t in range(2)
)
# Represent system of equations using "x" as defined above.
# I use the intermediate "exprs" variable to allow breaking up the tuple
# over multiple lines.
exprs = (
x[0]^2 + x[1]^2 - 1,
x[0] - x[1]^2
)
f(*x) = exprs
# Guess values
guess = (1, 1)
print scipy.optimize.fsolve(lambda v: f(*map(float, v)), guess)
</code></pre>
<p>In all cases, the result should look like:
<code>[ 0.61803399 0.78615138]</code></p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21866#post-id-21866Thank you very much for the tips!Mon, 11 Apr 2011 07:50:48 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21866#post-id-21866Comment by cswiercz for <p>The original scenario involves two variables and only one equation. Because <code>ff</code> accepts two variables but does not return a list (or tuple, or array, ...) of length 2, the "shape mismatch" error is thrown. Here is a working example (with two equations and minimal modifications):</p>
<pre><code>import scipy.optimize
var('x,y')
f = x^2 + y^2 - 1
g = x - y^2
def ff(v):
x, y = map(float, v)
return [f(x=x, y=y), g(x=x, y=y)]
print scipy.optimize.fsolve(ff, [1, 1])
</code></pre>
<p>Here is a more compact version using lambda functions and the trick described link:<a href="http://www.sagemath.org/doc/reference/sage/modules/vector_callable_symbolic_dense.html" title="test">here</a></p>
<pre><code>import scipy.optimize
var('x,y')
f(x,y) = (x^2 + y^2 - 1, x - y^2)
print scipy.optimize.fsolve(lambda v: f(*map(float, v)), (1, 1))
</code></pre>
<p>Since you said you have a large number of variables and equations, here is a slightly more "scale-friendly" approach (but somewhat kludgy):</p>
<pre><code>import scipy.optimize
# Create symbolic variables "x_0" and "x_1".
# This is a pretty messy kludge
x = tuple(
SR.var(
"".join(['x_', str(t)])
)
for t in range(2)
)
# Represent system of equations using "x" as defined above.
# I use the intermediate "exprs" variable to allow breaking up the tuple
# over multiple lines.
exprs = (
x[0]^2 + x[1]^2 - 1,
x[0] - x[1]^2
)
f(*x) = exprs
# Guess values
guess = (1, 1)
print scipy.optimize.fsolve(lambda v: f(*map(float, v)), guess)
</code></pre>
<p>In all cases, the result should look like:
<code>[ 0.61803399 0.78615138]</code></p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21869#post-id-21869I realize now that my toy example doesn't accurately reflect my real situation where I have as many nonlinear equations as I have unknowns. I'll try using your suggestion. Mon, 11 Apr 2011 06:00:14 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21869#post-id-21869Comment by cswiercz for <p>The original scenario involves two variables and only one equation. Because <code>ff</code> accepts two variables but does not return a list (or tuple, or array, ...) of length 2, the "shape mismatch" error is thrown. Here is a working example (with two equations and minimal modifications):</p>
<pre><code>import scipy.optimize
var('x,y')
f = x^2 + y^2 - 1
g = x - y^2
def ff(v):
x, y = map(float, v)
return [f(x=x, y=y), g(x=x, y=y)]
print scipy.optimize.fsolve(ff, [1, 1])
</code></pre>
<p>Here is a more compact version using lambda functions and the trick described link:<a href="http://www.sagemath.org/doc/reference/sage/modules/vector_callable_symbolic_dense.html" title="test">here</a></p>
<pre><code>import scipy.optimize
var('x,y')
f(x,y) = (x^2 + y^2 - 1, x - y^2)
print scipy.optimize.fsolve(lambda v: f(*map(float, v)), (1, 1))
</code></pre>
<p>Since you said you have a large number of variables and equations, here is a slightly more "scale-friendly" approach (but somewhat kludgy):</p>
<pre><code>import scipy.optimize
# Create symbolic variables "x_0" and "x_1".
# This is a pretty messy kludge
x = tuple(
SR.var(
"".join(['x_', str(t)])
)
for t in range(2)
)
# Represent system of equations using "x" as defined above.
# I use the intermediate "exprs" variable to allow breaking up the tuple
# over multiple lines.
exprs = (
x[0]^2 + x[1]^2 - 1,
x[0] - x[1]^2
)
f(*x) = exprs
# Guess values
guess = (1, 1)
print scipy.optimize.fsolve(lambda v: f(*map(float, v)), guess)
</code></pre>
<p>In all cases, the result should look like:
<code>[ 0.61803399 0.78615138]</code></p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21867#post-id-21867My other question is whether or not this is the most efficient way to approach the problem or if there is (hopefully) a one-liner that will Python function wrap my symbolic functions. I suppose I could use lambda notation but at this point I suppose I'm just talking to myself. :) Thanks for your response!Mon, 11 Apr 2011 06:36:32 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21867#post-id-21867Answer by DSM for <p>I need to find the roots of a system of multivariate non-linear (algebraic) functions and Sage's <code>solve</code> function is running indefinitely. (Where "indefinitely" means "more than 30 mins".) I only require numerical solutions so I was hoping to use Scipy's <code>fsolve</code> function instead. However, the functions that I'm generating are Sage symbolic functions, which have been mightly convenient to work with, and I'm having trouble getting them into a format that <code>fsolve</code> will understand. Does anyone have any suggestions?</p>
<p>The documentation for Scipy's <code>fsolve</code> can be found <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve">here</a>. It looks like a Python function is sufficient.</p>
<p>As a toy example, I've tried the following:</p>
<pre><code>sage: from scipy.optimize import fsolve
sage: var('x,y')
sage: f = x^2 + y^2
sage: def ff(v):
return f(x=v[0],y=v[0])
sage: fsolve(ff,[1,1])
</code></pre>
<p>It seems like the naive approach but I'm receiving the error <code>TypeError: no canonical coercion from <type 'numpy.int64'> to Symbolic Ring</code>. Perhaps this is an issue with Sage properly dealing with Numpy/Scipy data types?</p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?answer=12291#post-id-12291The "minimize" wrapper has better Sage integration (e.g. symbolic gradients). The equivalent of the above example would be something like
sage: var("x, y")
(x, y)
sage: f = x^2 + y^2 - 1
sage: g = x - y^2
sage: eqs = f, g
sage: minimize(norm(vector(eqs)), [1.0, 1.0])
Divide-by-zero encountered: rhok assumed large
Warning: Desired error not necessarily achieveddue to precision loss
Current function value: 0.000000
Iterations: 57
Function evaluations: 124
Gradient evaluations: 123
(0.61803398875, 0.786151377757)
So if you wanted a one-liner, something along those lines might do it.Tue, 12 Apr 2011 18:51:48 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?answer=12291#post-id-12291Comment by cswiercz for <p>The "minimize" wrapper has better Sage integration (e.g. symbolic gradients). The equivalent of the above example would be something like</p>
<pre><code>sage: var("x, y")
(x, y)
sage: f = x^2 + y^2 - 1
sage: g = x - y^2
sage: eqs = f, g
sage: minimize(norm(vector(eqs)), [1.0, 1.0])
Divide-by-zero encountered: rhok assumed large
Warning: Desired error not necessarily achieveddue to precision loss
Current function value: 0.000000
Iterations: 57
Function evaluations: 124
Gradient evaluations: 123
(0.61803398875, 0.786151377757)
</code></pre>
<p>So if you wanted a one-liner, something along those lines might do it.</p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21847#post-id-21847Of course. D'oh.Wed, 13 Apr 2011 06:17:04 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21847#post-id-21847Comment by DSM for <p>The "minimize" wrapper has better Sage integration (e.g. symbolic gradients). The equivalent of the above example would be something like</p>
<pre><code>sage: var("x, y")
(x, y)
sage: f = x^2 + y^2 - 1
sage: g = x - y^2
sage: eqs = f, g
sage: minimize(norm(vector(eqs)), [1.0, 1.0])
Divide-by-zero encountered: rhok assumed large
Warning: Desired error not necessarily achieveddue to precision loss
Current function value: 0.000000
Iterations: 57
Function evaluations: 124
Gradient evaluations: 123
(0.61803398875, 0.786151377757)
</code></pre>
<p>So if you wanted a one-liner, something along those lines might do it.</p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21848#post-id-21848That's why I put in the "norm" call.Wed, 13 Apr 2011 04:21:57 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21848#post-id-21848Comment by cswiercz for <p>The "minimize" wrapper has better Sage integration (e.g. symbolic gradients). The equivalent of the above example would be something like</p>
<pre><code>sage: var("x, y")
(x, y)
sage: f = x^2 + y^2 - 1
sage: g = x - y^2
sage: eqs = f, g
sage: minimize(norm(vector(eqs)), [1.0, 1.0])
Divide-by-zero encountered: rhok assumed large
Warning: Desired error not necessarily achieveddue to precision loss
Current function value: 0.000000
Iterations: 57
Function evaluations: 124
Gradient evaluations: 123
(0.61803398875, 0.786151377757)
</code></pre>
<p>So if you wanted a one-liner, something along those lines might do it.</p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21849#post-id-21849Thanks for the one-liner suggestion! Since I'm looking for roots and my functions aren't necessarily non-negative I might have to take absolute values for this method to work. I'll try it out in my real problem.Wed, 13 Apr 2011 03:10:05 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21849#post-id-21849Answer by Felix Lawrence for <p>I need to find the roots of a system of multivariate non-linear (algebraic) functions and Sage's <code>solve</code> function is running indefinitely. (Where "indefinitely" means "more than 30 mins".) I only require numerical solutions so I was hoping to use Scipy's <code>fsolve</code> function instead. However, the functions that I'm generating are Sage symbolic functions, which have been mightly convenient to work with, and I'm having trouble getting them into a format that <code>fsolve</code> will understand. Does anyone have any suggestions?</p>
<p>The documentation for Scipy's <code>fsolve</code> can be found <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve">here</a>. It looks like a Python function is sufficient.</p>
<p>As a toy example, I've tried the following:</p>
<pre><code>sage: from scipy.optimize import fsolve
sage: var('x,y')
sage: f = x^2 + y^2
sage: def ff(v):
return f(x=v[0],y=v[0])
sage: fsolve(ff,[1,1])
</code></pre>
<p>It seems like the naive approach but I'm receiving the error <code>TypeError: no canonical coercion from <type 'numpy.int64'> to Symbolic Ring</code>. Perhaps this is an issue with Sage properly dealing with Numpy/Scipy data types?</p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?answer=12268#post-id-12268This is one of the problems you're having:
sage: ff(scipy.array([1,1]))
TypeError: no canonical coercion from <type 'numpy.int64'> to Symbolic Ring
Which can be remedied as such:
sage: ff(scipy.array([1.,1.]))
2
But this merely exposes another problem:
sage: scipy.optimize.fsolve(ff, scipy.array([1.,1.]))
TypeError: There is a mismatch between the input and output shape of ff.
And I don't know how to solve that one.Sun, 10 Apr 2011 20:36:18 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?answer=12268#post-id-12268Comment by cswiercz for <p>This is one of the problems you're having:</p>
<pre><code>sage: ff(scipy.array([1,1]))
TypeError: no canonical coercion from <type 'numpy.int64'> to Symbolic Ring
</code></pre>
<p>Which can be remedied as such:</p>
<pre><code>sage: ff(scipy.array([1.,1.]))
2
</code></pre>
<p>But this merely exposes another problem:</p>
<pre><code>sage: scipy.optimize.fsolve(ff, scipy.array([1.,1.]))
TypeError: There is a mismatch between the input and output shape of ff.
</code></pre>
<p>And I don't know how to solve that one.</p>
http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21868#post-id-21868I think Kelvin Li appropriately identified the issue with that last TypeError. Exploring some more of the Scipy documentation shows no examples for "underdetermined" systems; systems which may still have local minima / maxima / solutions.Mon, 11 Apr 2011 06:33:52 -0500http://ask.sagemath.org/question/8069/using-sage-symbolic-functions-in-scipy-fsolve/?comment=21868#post-id-21868