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, 25 Nov 2015 10:47:21 -0600Suggestions for improving speedhttp://ask.sagemath.org/question/30974/suggestions-for-improving-speed/ I'm working on visualizing a flow in R^3, and specifically focusing on how the flow behaves near a torus on inner radius 1 and outer radius 2, centered at the origin. The idea is to take a point (x,y,z), and look at this point using coordinates (r,z), in the plane containing (x,y,z) and the origin. Then I take equations in polar coordinates in this plane, and apply them to the point. Moving back into R^3, I take a rotation of the point, and add in the result of the flow restricted to the plane. Most of these conditions also require that I know the minimum distance from (x,y,z) to the torus mentioned.
This is taking a very long time, especially when I try to animate the flow. The last animation I tried took 20 hours, and that was before I found a flaw in the construction. This one works, but I'm hoping to optimize it before I try it again. Any suggestions on how to speed things up would be appreciated.
x,y,u,v = var('x y u v')
ep = 0.1 #step size
length = 350 #iterations
#initial conditions
omega_seeds = [(-1,-1,2),(1,-1,2),(-1,1,2),(1,1,2),(2,2,2),(-2,2,2),(2,-2,2),(-2,-2,2)]
twopi = 2 * RDF.pi()
c(r) = 2/3*(4**(r+1)-4)
def norm2((x,y)):
"""
Norm of `(x,y)`
"""
return float(sqrt((x^2 + y^2)))
def norm3((x,y,z)):
"""
Norm squared of `(x,y,z)`
"""
return x^2 + y^2 + z^2
def norm((x,y,z)):
"""
Norm of `(x,y,z)`
"""
return float(sqrt(norm3((x,y,z))))
def dist((x,y,z),(xx,yy,zz)):
"""
Distance from `(x,y,z)` to `(xx,yy,zz)`
"""
return norm((x - xx, y - yy,z-zz))
def atan2((x,y)): #generalized arctangent
if x == 0 and y < 0:
return 3*pi/2
elif x == 0 and y > 0:
return pi/2
else:
a = float(arctan(abs(y/x)))
if x > 0:
if y > 0:
return a
elif y < 0:
return -a
elif y == 0:
return 0
elif x < 0:
if y > 0:
return pi-a
elif y < 0:
return a-pi
#Functions to determine distances of a point from a specified torus
def Torus_dist((x,y,z),n):
A = vector([0,c(n),0])
if n%2 == 0:
B = vector([x,y,0])
else:
B = vector([0,y,z])
if A == B:
return 4^n
else:
t = (2*4^n)/(norm(A - B))
v = t*B+ (1-t)*A
C = vector([v[0],v[1],v[2]])
if dist(C,(x,y,z)) <= 4^n:
return 0
else:
D = vector([x,y,z])
t2 = 4^n/norm(C-D)
v2 = t2*D + (1-t2)*C
E = vector([v2[0],v2[1],v2[2]])
return dist((x,y,z),E)
def dz(x,y): #vertical component of the flow restricted to the plane
r = norm2((x,y))
if r <=1:
return 0
if r >=2:
return y
theta = atan2((x,y))
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
return float(dr*cos(theta)-r*sin(theta)*dtheta)
def dr(x,y): #horizontal component of the flow restricted to the plane
r = norm2((x,y))
if r <=1:
return 0
if r >=2:
return 0
theta = atan2((x,y))
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
return float(dr*sin(theta) + r*cos(theta)*dtheta)
def W1((x,y,z)): #dx for the points in R^3
r = 2 - norm2((x,y))
T = max(0,1-Torus_dist((x,y,z),0))
theta = atan2((x,y))
if 0<=z<=2:
return T*y + dr(r,z)*cos(theta)
elif -2<=z<=0:
return (-1)*T*y + dr(r,z)*cos(theta)
else:
return 0
def W2((x,y,z)): #dy for the points in R^3
r = 2 - norm2((x,y))
T = max(0,1-Torus_dist((x,y,z),0))
theta = atan2((x,y))
if 0<=z<=2:
return (-1)*T*x + dr(r,z)*sin(theta)
elif -2<=z<=0:
return T*x + dr(r,z)*sin(theta)
else:
return 0
def W3((x,y,z)): #dz for the points in R^3
r = 2 - norm2((x,y))
if -2<=z<=2:
return -dz(r,z)
else:
return 0
def Plugplot((x,y,z)): #takes a point in R^3 and returns a point epsilon distance away, in the direction of the gradient at that point
n = norm((x - W1((x,y,z)), y - W2((x,y,z)),z-W3((x,y,z))))*(1/ep)
if n == 0:
return ((x,y,z))
else:
return ((x + W1((x,y,z))/n), (y + W2((x,y,z))/n), (z + W3((x,y,z))/n))
def PlugFlow((x,y,z)): #creates a collection of lines through the points generated by Plugplot
L = []
ic = (x,y,z)
if z >= -2:
for i in range(length):
L.append(line((ic,Plugplot(ic)),color='blue'))
ic = Plugplot(ic)
return L
a,b = 2,1 # outer and inner radii
Tx = (a + b*cos(u))*cos(v)
Ty = (a + b*cos(u))*sin(v)
Tz = b*sin(u)
T = parametric_plot3d([Tx,Ty,Tz], (u,0,twopi),(v,0,twopi), opacity=.2,aspect_ratio=1,color='black')
#Adds all of the to an array with the torus T
L = [T]
for i in omega_seeds:
L = L + PlugFlow(i)
show(sum(L) ,frame=true,aspect_ratio=1)Sun, 22 Nov 2015 16:34:58 -0600http://ask.sagemath.org/question/30974/suggestions-for-improving-speed/Answer by slelievre for <div class="snippet"><p>I'm working on visualizing a flow in R^3, and specifically focusing on how the flow behaves near a torus on inner radius 1 and outer radius 2, centered at the origin. The idea is to take a point (x,y,z), and look at this point using coordinates (r,z), in the plane containing (x,y,z) and the origin. Then I take equations in polar coordinates in this plane, and apply them to the point. Moving back into R^3, I take a rotation of the point, and add in the result of the flow restricted to the plane. Most of these conditions also require that I know the minimum distance from (x,y,z) to the torus mentioned.</p>
<p>This is taking a very long time, especially when I try to animate the flow. The last animation I tried took 20 hours, and that was before I found a flaw in the construction. This one works, but I'm hoping to optimize it before I try it again. Any suggestions on how to speed things up would be appreciated.</p>
<pre><code>x,y,u,v = var('x y u v')
ep = 0.1 #step size
length = 350 #iterations
#initial conditions
omega_seeds = [(-1,-1,2),(1,-1,2),(-1,1,2),(1,1,2),(2,2,2),(-2,2,2),(2,-2,2),(-2,-2,2)]
twopi = 2 * RDF.pi()
c(r) = 2/3*(4**(r+1)-4)
def norm2((x,y)):
"""
Norm of `(x,y)`
"""
return float(sqrt((x^2 + y^2)))
def norm3((x,y,z)):
"""
Norm squared of `(x,y,z)`
"""
return x^2 + y^2 + z^2
def norm((x,y,z)):
"""
Norm of `(x,y,z)`
"""
return float(sqrt(norm3((x,y,z))))
def dist((x,y,z),(xx,yy,zz)):
"""
Distance from `(x,y,z)` to `(xx,yy,zz)`
"""
return norm((x - xx, y - yy,z-zz))
def atan2((x,y)): #generalized arctangent
if x == 0 and y < 0:
return 3*pi/2
elif x == 0 and y > 0:
return pi/2
else:
a = float(arctan(abs(y/x)))
if x > 0:
if y > 0:
return a
elif y < 0:
return -a
elif y == 0:
return 0
elif x < 0:
if y > 0:
return pi-a
elif y < 0:
return a-pi
#Functions to determine distances of a point from a specified torus
def Torus_dist((x,y,z),n):
A = vector([0,c(n),0])
if n%2 == 0:
B = vector([x,y,0])
else:
B = vector([0,y,z])
if A == B:
return 4^n
else:
t = (2*4^n)/(norm(A - B))
v = t*B+ (1-t)*A
C = vector([v[0],v[1],v[2]])
if dist(C,(x,y,z)) <= 4^n:
return 0
else:
D = vector([x,y,z])
t2 = 4^n/norm(C-D)
v2 = t2*D + (1-t2)*C
E = vector([v2[0],v2[1],v2[2]])
return dist((x,y,z),E)
def dz(x,y): #vertical component of the flow restricted to the plane
r = norm2((x,y))
if r <=1 ...</code></pre><span class="expander"> <a>(more)</a></span></div> http://ask.sagemath.org/question/30974/suggestions-for-improving-speed/?answer=31058#post-id-31058I fully agree with @tmonteil's answer, in particular
- avoiding Sage's `SymbolicRing` should be a huge gain
- using `RDF` should be much faster than `RR` too.
Two small additional comments:
- for pi, I would suggest using `RDF.pi()`
- instead of the norm, using the square of the norm
saves extracting square roots, eg instead of testing
whether `sqrt(x^2+y^2) < 2`, test whether `x^2 + y^2 < 4`.
(You should test whether this really saves you time...
it could be that the `norm` method calls faster code
than code you might write yourself for `norm_squared`)Wed, 25 Nov 2015 10:47:21 -0600http://ask.sagemath.org/question/30974/suggestions-for-improving-speed/?answer=31058#post-id-31058Answer by tmonteil for <div class="snippet"><p>I'm working on visualizing a flow in R^3, and specifically focusing on how the flow behaves near a torus on inner radius 1 and outer radius 2, centered at the origin. The idea is to take a point (x,y,z), and look at this point using coordinates (r,z), in the plane containing (x,y,z) and the origin. Then I take equations in polar coordinates in this plane, and apply them to the point. Moving back into R^3, I take a rotation of the point, and add in the result of the flow restricted to the plane. Most of these conditions also require that I know the minimum distance from (x,y,z) to the torus mentioned.</p>
<p>This is taking a very long time, especially when I try to animate the flow. The last animation I tried took 20 hours, and that was before I found a flaw in the construction. This one works, but I'm hoping to optimize it before I try it again. Any suggestions on how to speed things up would be appreciated.</p>
<pre><code>x,y,u,v = var('x y u v')
ep = 0.1 #step size
length = 350 #iterations
#initial conditions
omega_seeds = [(-1,-1,2),(1,-1,2),(-1,1,2),(1,1,2),(2,2,2),(-2,2,2),(2,-2,2),(-2,-2,2)]
twopi = 2 * RDF.pi()
c(r) = 2/3*(4**(r+1)-4)
def norm2((x,y)):
"""
Norm of `(x,y)`
"""
return float(sqrt((x^2 + y^2)))
def norm3((x,y,z)):
"""
Norm squared of `(x,y,z)`
"""
return x^2 + y^2 + z^2
def norm((x,y,z)):
"""
Norm of `(x,y,z)`
"""
return float(sqrt(norm3((x,y,z))))
def dist((x,y,z),(xx,yy,zz)):
"""
Distance from `(x,y,z)` to `(xx,yy,zz)`
"""
return norm((x - xx, y - yy,z-zz))
def atan2((x,y)): #generalized arctangent
if x == 0 and y < 0:
return 3*pi/2
elif x == 0 and y > 0:
return pi/2
else:
a = float(arctan(abs(y/x)))
if x > 0:
if y > 0:
return a
elif y < 0:
return -a
elif y == 0:
return 0
elif x < 0:
if y > 0:
return pi-a
elif y < 0:
return a-pi
#Functions to determine distances of a point from a specified torus
def Torus_dist((x,y,z),n):
A = vector([0,c(n),0])
if n%2 == 0:
B = vector([x,y,0])
else:
B = vector([0,y,z])
if A == B:
return 4^n
else:
t = (2*4^n)/(norm(A - B))
v = t*B+ (1-t)*A
C = vector([v[0],v[1],v[2]])
if dist(C,(x,y,z)) <= 4^n:
return 0
else:
D = vector([x,y,z])
t2 = 4^n/norm(C-D)
v2 = t2*D + (1-t2)*C
E = vector([v2[0],v2[1],v2[2]])
return dist((x,y,z),E)
def dz(x,y): #vertical component of the flow restricted to the plane
r = norm2((x,y))
if r <=1 ...</code></pre><span class="expander"> <a>(more)</a></span></div> http://ask.sagemath.org/question/30974/suggestions-for-improving-speed/?answer=30976#post-id-30976I am not sure if you spend more time in just drawing pictures or in the computations. If the plot is guilty, you should look for another Python library that will be faster, but i am not specialist about pictures and animations. You can use `%timeit` to see which part of the program takes the most time.
If the computations take too much time, you should get rid of symbolic computations, and use floating-point arithmetic instead, using `RDF` where needed. Be careful that the coercion goes toward the symbolic ring, so each time you use it, the follogin computations use it as well. I did not run your code, but you can start with the following and tell if this solves your problem and if you need more hints:
1. First step remove the `x,y,u,v = var('x y u v')` at the beginning, and see for which function there is a problem, then replace those symbolic functions to PYthon function.
1. replace `ep = 0.1` with `ep = RDF(0.1)` (`Real Double Field` is faster than `Real Field with 53 bits of precision`, and the coercion goes from the first to the second).
1. replace `c(r) = 2/3*(4**(r+1)-4)` which defines a symbolic function with the following Python function : `c = lambda r : 2/3*(4**(r+1)-4)`
1. Do not use `pi`, use `pi_approx = RDF(pi)`
1. note that vectors have a `.norm()` method already defined (note also that you can specify the base ring of the vector),
for example:
sage: v = vector([1,2])
sage: v.norm()
sqrt(5) # bad
sage: v = vector(RDF, [3,4])
sage: v.norm()
5.0
I think you get the idea.
Sun, 22 Nov 2015 17:03:28 -0600http://ask.sagemath.org/question/30974/suggestions-for-improving-speed/?answer=30976#post-id-30976Comment by Nathann for <p>I am not sure if you spend more time in just drawing pictures or in the computations. If the plot is guilty, you should look for another Python library that will be faster, but i am not specialist about pictures and animations. You can use <code>%timeit</code> to see which part of the program takes the most time.</p>
<p>If the computations take too much time, you should get rid of symbolic computations, and use floating-point arithmetic instead, using <code>RDF</code> where needed. Be careful that the coercion goes toward the symbolic ring, so each time you use it, the follogin computations use it as well. I did not run your code, but you can start with the following and tell if this solves your problem and if you need more hints:</p>
<ol>
<li><p>First step remove the <code>x,y,u,v = var('x y u v')</code> at the beginning, and see for which function there is a problem, then replace those symbolic functions to PYthon function.</p></li>
<li><p>replace <code>ep = 0.1</code> with <code>ep = RDF(0.1)</code> (<code>Real Double Field</code> is faster than <code>Real Field with 53 bits of precision</code>, and the coercion goes from the first to the second).</p></li>
<li><p>replace <code>c(r) = 2/3*(4**(r+1)-4)</code> which defines a symbolic function with the following Python function : <code>c = lambda r : 2/3*(4**(r+1)-4)</code></p></li>
<li><p>Do not use <code>pi</code>, use <code>pi_approx = RDF(pi)</code></p></li>
<li><p>note that vectors have a <code>.norm()</code> method already defined (note also that you can specify the base ring of the vector), </p></li>
</ol>
<p>for example:</p>
<pre><code>sage: v = vector([1,2])
sage: v.norm()
sqrt(5) # bad
sage: v = vector(RDF, [3,4])
sage: v.norm()
5.0
</code></pre>
<p>I think you get the idea.</p>
http://ask.sagemath.org/question/30974/suggestions-for-improving-speed/?comment=30982#post-id-30982Here is a tutorial on profiling in Sage, in case you want to figure out what exactly is going too slow: http://doc.sagemath.org/html/en/thematic_tutorials/profiling.html#profilingMon, 23 Nov 2015 02:26:56 -0600http://ask.sagemath.org/question/30974/suggestions-for-improving-speed/?comment=30982#post-id-30982