Ask Your Question

Revision history [back]

In this case you could try to substitute the values in the exact solution. An example of this sort is given in pages 46-40 of Sage for power users (1).

Here is a minimal working example:

var('x, y, x0, y0, x1, y1, d0, d1')

# two circles
eq1 = (x-x0)^2 +(y-y0)^2 == d0^2
eq2 = (x-x1)^2 + (y-y1)^2 == d1^2

Let's ask Sage to solve exactly this system of equations:

%time solns = solve([eq1, eq2], x, y)
CPU times: user 7.28 s, sys: 8.01 ms, total: 7.28 s
Wall time: 7.28 s

It takes a couple of seconds, but it does a good job in finding the exact formulas,

sage: solns[0][0].rhs().variables()
(d0, d1, x0, x1, y0, y1)

Of course, the advantage is that the solutions obtained with this method are exact,

sage: solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1)
-1/8*sqrt(7) + 19/8

Let's see how much it takes to evaluate it numerically:

%timeit _ = solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1).n()
1000 loops, best of 3: 704 µs per loop

Now let's wrap the solution into a Cython function (1, 2):

%%cython
from math import sqrt
def sol1_x(double x0, double y0, double x1, double y1, double d0, double d1):
    return -0.5*(d0**2*x0 - d1**2*x0 - x0**3 + x0*x1**2 - x1**3 - x0*y0**2 - (x0 + x1)*y1**2 - (d0**2 - d1**2 - x0**2 + y0**2)*x1 + (2*x0*y0 + 2*x1*y0 + sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1))*y1 - sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1)*y0)/(x0**2 - 2*x0*x1 + x1**2 + y0**2 - 2*y0*y1 + y1**2)

Then,

%timeit _ = sol1_x(1.,2.,3.,4.,2.,1.)
100000 loops, best of 3: 7.03 µs per loop
sol1_x(1.,2.,3.,4.,2.,1.)
2.0442810861169263

That's a factor 100x of faster!

(1) if you are in the Jupyter notebook, create a new Code cell and in the first line use %%cython (no %cython).

(2) Two observations:

  • replace the operator ^ with **, for example with the command str(solns[0][0].rhs()).replace('^', '**')
  • use 0.5 instead of 1/2

Finally, Sage does provide a function to do this automatically, (eg. operating with Sage objects defined in the global scope), it is cython_lambda. I wasn't unable to use it properly in this case.. strange.

If you have to evaluate thousands of values, then it will be even more efficient to work with arrays, see the Cython documentation for more info or ask here :)

In this case you could try to substitute the values in the exact solution. An example of this sort is given in pages 46-40 46-50 of Sage for power users (1).

Here is a minimal working example:

var('x, y, x0, y0, x1, y1, d0, d1')

# two circles
eq1 = (x-x0)^2 +(y-y0)^2 == d0^2
eq2 = (x-x1)^2 + (y-y1)^2 == d1^2

Let's ask Sage to solve exactly this system of equations:

%time solns = solve([eq1, eq2], x, y)
CPU times: user 7.28 s, sys: 8.01 ms, total: 7.28 s
Wall time: 7.28 s

It takes a couple of seconds, but it does a good job in finding the exact formulas,

sage: solns[0][0].rhs().variables()
(d0, d1, x0, x1, y0, y1)

Of course, the advantage is that the solutions obtained with this method are exact,

sage: solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1)
-1/8*sqrt(7) + 19/8

Let's see how much it takes to evaluate it numerically:

%timeit _ = solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1).n()
1000 loops, best of 3: 704 µs per loop

Now let's wrap the solution into a Cython function (1, 2):

%%cython
from math import sqrt
def sol1_x(double x0, double y0, double x1, double y1, double d0, double d1):
    return -0.5*(d0**2*x0 - d1**2*x0 - x0**3 + x0*x1**2 - x1**3 - x0*y0**2 - (x0 + x1)*y1**2 - (d0**2 - d1**2 - x0**2 + y0**2)*x1 + (2*x0*y0 + 2*x1*y0 + sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1))*y1 - sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1)*y0)/(x0**2 - 2*x0*x1 + x1**2 + y0**2 - 2*y0*y1 + y1**2)

Then,

%timeit _ = sol1_x(1.,2.,3.,4.,2.,1.)
100000 loops, best of 3: 7.03 µs per loop
sol1_x(1.,2.,3.,4.,2.,1.)
2.0442810861169263

That's a factor 100x of faster!

(1) if you are in the Jupyter notebook, create a new Code cell and in the first line use %%cython (no %cython).

(2) Two observations:

  • replace the operator ^ with **, for example with the command str(solns[0][0].rhs()).replace('^', '**')
  • use 0.5 instead of 1/2

Finally, Sage does provide a function to do this automatically, (eg. operating with Sage objects defined in the global scope), it is cython_lambda. I wasn't unable to use it properly in this case.. strange.

If you have to evaluate thousands of values, then it will be even more efficient to work with arrays, see the Cython documentation for more info or ask here :)

In this case you could try to substitute the values in the exact solution. An example of this sort is given in pages 46-50 of Sage for power users (1).

Here is a minimal working example:

var('x, y, x0, y0, x1, y1, d0, d1')

# two circles
eq1 = (x-x0)^2 +(y-y0)^2 == d0^2
eq2 = (x-x1)^2 + (y-y1)^2 == d1^2

Let's ask Sage to solve exactly this system of equations:

%time solns = solve([eq1, eq2], x, y)
CPU times: user 7.28 s, sys: 8.01 ms, total: 7.28 s
Wall time: 7.28 s

It takes a couple of seconds, but it does a good job in finding the exact formulas,

sage: solns[0][0].rhs().variables()
(d0, d1, x0, x1, y0, y1)

Of course, the advantage is that the solutions obtained with this method are exact,

sage: solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1)
-1/8*sqrt(7) + 19/8

Let's see how much it takes to evaluate it numerically:

%timeit _ = solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1).n()
1000 loops, best of 3: 704 µs per loop

Now let's wrap the solution into a Cython function (1, 2):(1), (2):

%%cython
from math import sqrt
def sol1_x(double x0, double y0, double x1, double y1, double d0, double d1):
    return -0.5*(d0**2*x0 - d1**2*x0 - x0**3 + x0*x1**2 - x1**3 - x0*y0**2 - (x0 + x1)*y1**2 - (d0**2 - d1**2 - x0**2 + y0**2)*x1 + (2*x0*y0 + 2*x1*y0 + sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1))*y1 - sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1)*y0)/(x0**2 - 2*x0*x1 + x1**2 + y0**2 - 2*y0*y1 + y1**2)

Then,

%timeit _ = sol1_x(1.,2.,3.,4.,2.,1.)
100000 loops, best of 3: 7.03 µs per loop
sol1_x(1.,2.,3.,4.,2.,1.)
2.0442810861169263

That's a factor 100x of faster!

(1) if you are in the Jupyter notebook, create a new Code cell and in the first line use %%cython (no (and not %cython).

(2) Two observations:

  • replace the operator ^ with **, for example with the command str(solns[0][0].rhs()).replace('^', '**')
  • use 0.5 instead of 1/2

Finally, Sage does provide a function to do this automatically, (eg. operating with Sage objects defined in the global scope), it is cython_lambda. I wasn't unable to use it properly in this case.. strange.

If you have to evaluate thousands of values, then it will be even more efficient to work with arrays, see the Cython documentation for more info or ask here :)

In this case you could try to substitute the values in the exact solution. An example of this sort is given in pages 46-50 of Sage for power users (1).

Here is a minimal working example:

var('x, y, x0, y0, x1, y1, d0, d1')

# two circles
eq1 = (x-x0)^2 +(y-y0)^2 == d0^2
eq2 = (x-x1)^2 + (y-y1)^2 == d1^2

Let's ask Sage to solve exactly this system of equations:

%time solns = solve([eq1, eq2], x, y)
CPU times: user 7.28 s, sys: 8.01 ms, total: 7.28 s
Wall time: 7.28 s

It takes a couple of seconds, but it does a good job in finding the exact formulas,

sage: solns[0][0].rhs().variables()
(d0, d1, x0, x1, y0, y1)

Of course, the advantage is that the solutions obtained with this method are exact,

sage: solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1)
-1/8*sqrt(7) + 19/8
sage: _.n()
2.04428108611693

Let's see how much it takes to evaluate it numerically:

%timeit _ = solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1).n()
1000 loops, best of 3: 704 µs per loop

Now let's wrap the solution into a Cython function (1), (2):

%%cython
from math import sqrt
def sol1_x(double x0, double y0, double x1, double y1, double d0, double d1):
    return -0.5*(d0**2*x0 - d1**2*x0 - x0**3 + x0*x1**2 - x1**3 - x0*y0**2 - (x0 + x1)*y1**2 - (d0**2 - d1**2 - x0**2 + y0**2)*x1 + (2*x0*y0 + 2*x1*y0 + sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1))*y1 - sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1)*y0)/(x0**2 - 2*x0*x1 + x1**2 + y0**2 - 2*y0*y1 + y1**2)

Then,

%timeit _ = sol1_x(1.,2.,3.,4.,2.,1.)
100000 loops, best of 3: 7.03 µs per loop
sol1_x(1.,2.,3.,4.,2.,1.)
2.0442810861169263

That's a factor 100x of faster!

(1) if you are in the Jupyter notebook, create a new Code cell and in the first line use %%cython (and not %cython).

(2) Two observations:

  • replace the operator ^ with **, for example with the command str(solns[0][0].rhs()).replace('^', '**')
  • use 0.5 instead of 1/2

Finally, Sage does provide a function to do this automatically, (eg. operating with Sage objects defined in the global scope), it is cython_lambda. I wasn't unable to use it properly in this case.. strange.

If you have to evaluate thousands of values, then it will be even more efficient to work with arrays, see the Cython documentation for more info or ask here :)

In this case you could try to substitute the values in the exact solution. An example of this sort is given in pages 46-50 of Sage for power users (1).

Here is a minimal working example:

var('x, y, x0, y0, x1, y1, d0, d1')

# two circles
eq1 = (x-x0)^2 +(y-y0)^2 == d0^2
eq2 = (x-x1)^2 + (y-y1)^2 == d1^2

Let's ask Sage to solve exactly this system of equations:

%time solns = solve([eq1, eq2], x, y)
CPU times: user 7.28 s, sys: 8.01 ms, total: 7.28 s
Wall time: 7.28 s

It takes a couple of seconds, but it does a good job in finding the exact formulas,

sage: solns[0][0].rhs().variables()
(d0, d1, x0, x1, y0, y1)

Of course, the advantage is that the solutions obtained with this method are exact,

sage: solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1)
-1/8*sqrt(7) + 19/8
sage: _.n()
2.04428108611693

Let's see how much it takes to evaluate it numerically:

%timeit _ = solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1).n()
1000 loops, best of 3: 704 µs per loop

Now let's wrap the solution into a Cython function (1), (2):

%%cython
from math import sqrt
def sol1_x(double x0, double y0, double x1, double y1, double d0, double d1):
    return -0.5*(d0**2*x0 - d1**2*x0 - x0**3 + x0*x1**2 - x1**3 - x0*y0**2 - (x0 + x1)*y1**2 - (d0**2 - d1**2 - x0**2 + y0**2)*x1 + (2*x0*y0 + 2*x1*y0 + sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1))*y1 - sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1)*y0)/(x0**2 - 2*x0*x1 + x1**2 + y0**2 - 2*y0*y1 + y1**2)

Then,

%timeit _ = sol1_x(1.,2.,3.,4.,2.,1.)
100000 loops, best of 3: 7.03 µs per loop
sol1_x(1.,2.,3.,4.,2.,1.)
2.0442810861169263

That's a factor 100x of faster!

(1) if you are in the Jupyter notebook, create a new Code cell and in the first line use %%cython (and not %cython).

(2) Two observations:

  • replace the operator ^ with **, for example with the command str(solns[0][0].rhs()).replace('^', '**')
  • use 0.5 instead of 1/2

Finally, Sage does provide a function to do this automatically, (eg. operating with Sage objects defined in the global scope), it is cython_lambda. I wasn't unable able to use apply it properly in this case.. here.. strange.

If you have to evaluate thousands of values, then it will be even more efficient to work with arrays, see the Cython documentation for more info or ask here :)

In this case you could try to substitute the values in the exact solution. An example of this sort is given in pages 46-50 of Sage for power users (1).

Here is a minimal working example:

var('x, y, x0, y0, x1, y1, d0, d1')

# two circles
eq1 = (x-x0)^2 +(y-y0)^2 == d0^2
eq2 = (x-x1)^2 + (y-y1)^2 == d1^2

Let's ask Sage to solve exactly this system of equations:

%time solns = solve([eq1, eq2], x, y)
CPU times: user 7.28 s, sys: 8.01 ms, total: 7.28 s
Wall time: 7.28 s

It takes a couple of seconds, but it does a good job in finding the exact formulas,

sage: solns[0][0].rhs().variables()
(d0, d1, x0, x1, y0, y1)

Of course, the advantage is that the solutions obtained with this method are exact,

sage: solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1)
-1/8*sqrt(7) + 19/8
sage: _.n()
2.04428108611693

Let's see how much it takes to evaluate it numerically:

%timeit _ = solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1).n()
1000 loops, best of 3: 704 µs per loop

Now let's wrap the solution into a Cython function (1), (2):

%%cython
from math import sqrt
def sol1_x(double x0, double y0, double x1, double y1, double d0, double d1):
    return -0.5*(d0**2*x0 - d1**2*x0 - x0**3 + x0*x1**2 - x1**3 - x0*y0**2 - (x0 + x1)*y1**2 - (d0**2 - d1**2 - x0**2 + y0**2)*x1 + (2*x0*y0 + 2*x1*y0 + sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1))*y1 - sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1)*y0)/(x0**2 - 2*x0*x1 + x1**2 + y0**2 - 2*y0*y1 + y1**2)

Then,

%timeit _ = sol1_x(1.,2.,3.,4.,2.,1.)
100000 loops, best of 3: 7.03 µs per loop
sol1_x(1.,2.,3.,4.,2.,1.)
2.0442810861169263

That's a factor 100x of faster!

(1) if you are in the Jupyter notebook, create a new Code cell and in the first line use %%cython (and not %cython).

(2) Two observations:

  • replace the operator ^ with **, for example with the command str(solns[0][0].rhs()).replace('^', '**')
  • use 0.5 instead of 1/2

Finally, Sage does provide a function to do this automatically, (eg. operating with Sage objects defined in the global scope), it is cython_lambda. I wasn't able to apply it here.. strange.

If you have to evaluate thousands of values, then it will be even more efficient to work with arrays, see the Cython documentation for more info or ask here :)

Without passing through the exact solution obtained with solve, here is the same answer with SciPy's fsolve interface:

from scipy.optimize import fsolve

def circles(variables, *args):
    x, y = variables
    x0, y0, x1, y1, d0, d1 = args
    eq1 = (x-x0)^2 +(y-y0)^2 - d0^2
    eq2 = (x-x1)^2 + (y-y1)^2 - d1^2
    return (eq1, eq2)


%timeit _ = fsolve(circles, (-2.0, 2.0), args=(1.,2.,3.,4.,2.,1.))
1000 loops, best of 3: 1.24 ms per loop
array([ 2.04428109,  3.70571891])

In this case you could try to substitute the values in the exact solution. An example of this sort is given in pages 46-50 of Sage for power users (1).

Here is a minimal working example:

var('x, y, x0, y0, x1, y1, d0, d1')

# two circles
eq1 = (x-x0)^2 +(y-y0)^2 == d0^2
eq2 = (x-x1)^2 + (y-y1)^2 == d1^2

Let's ask Sage to solve exactly this system of equations:

%time solns = solve([eq1, eq2], x, y)
CPU times: user 7.28 s, sys: 8.01 ms, total: 7.28 s
Wall time: 7.28 s

It takes a couple of seconds, but it does a good job in finding the exact formulas,

sage: solns[0][0].rhs().variables()
(d0, d1, x0, x1, y0, y1)

Of course, the advantage is that the solutions obtained with this method are exact,

sage: solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1)
-1/8*sqrt(7) + 19/8
sage: _.n()
2.04428108611693

Let's see how much it takes to evaluate it numerically:

%timeit _ = solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1).n()
1000 loops, best of 3: 704 µs per loop

Now let's wrap the solution into a Cython function (1), (2):

%%cython
from math import sqrt
def sol1_x(double x0, double y0, double x1, double y1, double d0, double d1):
    return -0.5*(d0**2*x0 - d1**2*x0 - x0**3 + x0*x1**2 - x1**3 - x0*y0**2 - (x0 + x1)*y1**2 - (d0**2 - d1**2 - x0**2 + y0**2)*x1 + (2*x0*y0 + 2*x1*y0 + sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1))*y1 - sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1)*y0)/(x0**2 - 2*x0*x1 + x1**2 + y0**2 - 2*y0*y1 + y1**2)

Then,

%timeit _ = sol1_x(1.,2.,3.,4.,2.,1.)
100000 loops, best of 3: 7.03 µs per loop
sol1_x(1.,2.,3.,4.,2.,1.)
2.0442810861169263

That's a factor 100x of faster!

(1) if you are in the Jupyter notebook, create a new Code cell and in the first line use %%cython (and not %cython).

(2) Two observations:

  • replace the operator ^ with **, for example with the command str(solns[0][0].rhs()).replace('^', '**')
  • use 0.5 instead of 1/2

Finally, Sage does provide a function to do this automatically, (eg. operating with Sage objects defined in the global scope), it is cython_lambda. I wasn't able to apply it here.. strange.

If you have to evaluate thousands of values, then it will be even more efficient to work with arrays, see the Cython documentation for more info or ask here :)

Without passing through the exact solution obtained with solve, here is the same answer with SciPy's fsolve interface:

from scipy.optimize import fsolve

def circles(variables, *args):
    x, y = variables
    x0, y0, x1, y1, d0, d1 = args
    eq1 = (x-x0)^2 +(y-y0)^2 - d0^2
    eq2 = (x-x1)^2 + (y-y1)^2 - d1^2
    return (eq1, eq2)


%timeit _ = fsolve(circles, (-2.0, 2.0), args=(1.,2.,3.,4.,2.,1.))
1000 loops, best of 3: 1.24 ms per loop

sage: fsolve(circles, (-2.0, 2.0), args=(1.,2.,3.,4.,2.,1.))
array([ 2.04428109,  3.70571891])
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
about | faq | help | privacy policy | terms of service
Powered by Askbot version 0.7.59