Another optimization request :-)

[jeffg]

If anyone wants to take this on... I would really really like to have
the spring_layout modified to support multi-threading if at all
possible.
My test data is 20,000, which makes this process 20,000 x 20,000 or
400,000,000 (400 million) calculations. This is taking somewhere
between 2-3 hours an iteration.
I plan to plot out over 1,000,000 data points or more, which would put
this at 1,000,000,000,000 (1 trillion) calculations. Any help in
making this more efficient would be much appreciated.

def spring_layout(G, iterations=50, dim=2, node_pos=None,
verbose=False):
"""Spring force model layout"""
if node_pos==None : # set the initial positions randomly in 1x1
box
vpos=random_layout(G, dim=dim)
else:
vpos=node_pos
if iterations==0:
return vpos
if G.order()==0:
k=1.0
else:
k=N.sqrt(1.0/G.order()) # optimal distance between nodes
disp={} # displacements

# initial "temperature" (about .1 of domain area)
# this is the largest step allowed in the dynamics
# linearly step down by dt on each iteration so
# on last iteration it is size dt.
t=0.1
dt=0.1/float(iterations+1)
for i in range(0,iterations):
for v in G:
if verbose==True:
print("Interation: " + str(i + 1) + ", Calculating: "
+ str(v.encode('iso-8859-15', "replace")))
disp[v]=N.zeros(dim)
for u in G:
delta=vpos[v]-vpos
dn=max(sqrt(N.dot(delta,delta)),0.01)
# repulsive force between all
deltaf=delta*k**2/dn**2
disp[v]=disp[v]+deltaf
# attractive force between neighbors
if G.has_edge(v,u):
deltaf=-delta*dn**2/(k*dn)
disp[v]=disp[v]+deltaf

# update positions
for v in G:
l=max(sqrt(N.dot(disp[v],disp[v])),0.01)
vpos[v]=vpos[v]+ disp[v]*t/l
t-=dt
return vpos

 

[jeffg]

sorry, that was stupid.  there is no inversion (apart from 1/m), just the
integration.

still, improving the integration would allow larger steps.

andrew

why are you dong this point by point?  surely you can express the physics
as a set of equations and invert the matrix?  wouldn't that be a lot
faster?  you'd replace the iteration over all combinations of points with
a faster matrix inversion.

see for example
http://www.medwelljournals.com/fulltext/ajit/2006/324-338.pdfpage 331
onwards.

there's a very nice into to the verlet integration mentioned here -
http://teknikus.dk/tj/gdc2001.htm

andrew

If anyone wants to take this on... I would really really like to have
the spring_layout modified to support multi-threading if at all
possible.
My test data is 20,000, which makes this process 20,000 x 20,000 or
400,000,000 (400 million) calculations.  This is taking somewhere
between 2-3 hours an iteration.
I plan to plot out over 1,000,000 data points or more, which would put
this at 1,000,000,000,000 (1 trillion) calculations.  Any help in
making this more efficient would be much appreciated.
def spring_layout(G, iterations=50, dim=2, node_pos=None,
verbose=False):
    """Spring force model layout"""
    if node_pos==None :  # set the initial positions randomly in 1x1
box
        vpos=random_layout(G, dim=dim)
    else:
        vpos=node_pos
    if iterations==0:
        return vpos
    if G.order()==0:
        k=1.0
    else:
        k=N.sqrt(1.0/G.order()) # optimal distance between nodes
    disp={}         # displacements
    # initial "temperature" (about .1 of domain area)
    # this is the largest step allowed in the dynamics
    # linearly step down by dt on each iteration so
    # on last iteration it is size dt.
    t=0.1
    dt=0.1/float(iterations+1)
    for i in range(0,iterations):
        for v in G:
            if verbose==True:
                print("Interation: " + str(i + 1) + ", Calculating: "
+ str(v.encode('iso-8859-15', "replace")))
            disp[v]=N.zeros(dim)
            for u in G:
                delta=vpos[v]-vpos
                dn=max(sqrt(N.dot(delta,delta)),0.01)
                # repulsive force between all
                deltaf=delta*k**2/dn**2
                disp[v]=disp[v]+deltaf
                # attractive force between neighbors
                if G.has_edge(v,u):
                    deltaf=-delta*dn**2/(k*dn)
                    disp[v]=disp[v]+deltaf
        # update positions
        for v in G:
            l=max(sqrt(N.dot(disp[v],disp[v])),0.01)
            vpos[v]=vpos[v]+ disp[v]*t/l
        t-=dt
    return vpos

To be honest, this is not my code and I'm new to python. It's part of
the open source project NetworkX, but I'm using this one call
extensively. I'm also not that familiar with the math behind the
physics. I'll read the documents and see if I can figure it
out. Thank you for the links and suggestions. I really need to
get this code performing at peak levels.

 

[Terry Reedy]

If anyone wants to take this on... I would really really like to have
the spring_layout modified to support multi-threading if at all
possible.
My test data is 20,000, which makes this process 20,000 x 20,000 or
400,000,000 (400 million) calculations. This is taking somewhere
between 2-3 hours an iteration.
I plan to plot out over 1,000,000 data points or more, which would put
this at 1,000,000,000,000 (1 trillion) calculations. Any help in
making this more efficient would be much appreciated.

def spring_layout(G, iterations=50, dim=2, node_pos=None,
verbose=False):
"""Spring force model layout"""
if node_pos==None : # set the initial positions randomly in 1x1
box
vpos=random_layout(G, dim=dim)
else:
vpos=node_pos
if iterations==0:
return vpos
if G.order()==0:
k=1.0
else:
k=N.sqrt(1.0/G.order()) # optimal distance between nodes
disp={} # displacements

# initial "temperature" (about .1 of domain area)
# this is the largest step allowed in the dynamics
# linearly step down by dt on each iteration so
# on last iteration it is size dt.
t=0.1
dt=0.1/float(iterations+1)
for i in range(0,iterations):
for v in G:
if verbose==True:
print("Interation: " + str(i + 1) + ", Calculating: "
+ str(v.encode('iso-8859-15', "replace")))
disp[v]=N.zeros(dim)
for u in G:
delta=vpos[v]-vpos

Let n = len(vpos). I believe there are only n-1 different deltas, dns,
and deltafs, where as n**2 calculations are made.

dn=max(sqrt(N.dot(delta,delta)),0.01)
# repulsive force between all
deltaf=delta*k**2/dn**2
disp[v]=disp[v]+deltaf
# attractive force between neighbors
if G.has_edge(v,u):
deltaf=-delta*dn**2/(k*dn)
disp[v]=disp[v]+deltaf

# update positions
for v in G:
l=max(sqrt(N.dot(disp[v],disp[v])),0.01)
vpos[v]=vpos[v]+ disp[v]*t/l
t-=dt
return vpos

That aside, array calculations like this should be *much* faster with
numpy. Or try compiling (with weave or Cython, for instance).

 

[bearophileHUGS]

If anyone wants to take this on... I would really really like to have
the spring_layout modified to support multi-threading if at all
possible.

You can start creating a compiled Python extension using ShedSkin, it
may be enough.

Bye,
bearophile