Ask Your Question
1

Suggestions for improving speed

asked 2015-11-22 23:34:58 +0100

jford1906 gravatar image

updated 2015-11-23 00:53:17 +0100

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 ...
(more)
edit retag flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted
2

answered 2015-11-23 00:03:28 +0100

tmonteil gravatar image

updated 2015-11-23 00:15:43 +0100

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 %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.

  2. 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).

  3. 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)

  4. Do not use pi, use pi_approx = RDF(pi)

  5. 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.

edit flag offensive delete link more

Comments

1

Here 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/thema...

Nathann gravatar imageNathann ( 2015-11-23 09:26:56 +0100 )edit
1

answered 2015-11-25 17:47:21 +0100

slelievre gravatar image

I 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)

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2015-11-22 23:34:58 +0100

Seen: 1,059 times

Last updated: Nov 25 '15