[Numpy-discussion] Stick (line segments) percolation algorithm - graph theory?

Josè Luis Mietta joseluismietta@yahoo.com...
Thu Aug 29 17:55:17 CDT 2013

Thanks a lot!!
José Luis

 De: Brett Olsen <brett.olsen@gmail.com>
Para: Discussion of Numerical Python <numpy-discussion@scipy.org> 
Enviado: lunes, 26 de agosto de 2013 14:08
Asunto: Re: [Numpy-discussion] Stick (line segments) percolation algorithm - graph theory?

I can see a couple opportunities for improvements in your algorithm.  Running your code on a single experiment, I get about 2.9 seconds to run. I get this down to about 1.0 seconds by (1) exploiting the symmetry of the M matrix and (2) avoiding the costly inner loop over k in favor of array operations:

def check_segments(j, others, data):
    x1, y1, x2, y2 = data
    x_A1B1 = x2[j]-x1[j]
    y_A1B1 = y2[j]-y1[j]
    x_A1A2 = x1[others]-x1[j]
    y_A1A2 = y1[others]-y1[j]      
    x_A2A1 = -1*x_A1A2
    y_A2A1 = -1*y_A1A2
    x_A2B2 = x2[others]-x1[others]
    y_A2B2 = y2[others]-y1[others]
    x_A1B2 = x2[others]-x1[j]
    y_A1B2 = y2[others]-y1[j]
    x_A2B1 = x2[j]-x1[others]
    y_A2B1 = y2[j]-y1[others]

    p1 = x_A1B1*y_A1A2 - y_A1B1*x_A1A2
    p2 = x_A1B1*y_A1B2 - y_A1B1*x_A1B2
    p3 = x_A2B2*y_A2B1 - y_A2B2*x_A2B1
    p4 = x_A2B2*y_A2A1 - y_A2B2*x_A2A1


    return (p1 * p2 <= 0) & (p3 * p4 <= 0)

for j in xrange(1, N):
    valid = check_segments(j, range(j), (x1, y1, x2, y2))
    M[j,0:j] = valid
    M[0:j,j] = valid

I don't see any other particularly simple ways to improve this.  You could probably add an interval check to ensure that the x and y intervals for the segments of interest overlap before doing the full check, but how much that would help would depend on the implementations.


On Fri, Aug 23, 2013 at 5:09 PM, Josè Luis Mietta <joseluismietta@yahoo.com.ar> wrote:

I wrote an algorithm for study stick percolation (i.e.: networks between line segments that intersect between them). In my algorithm N sticks (line segments) are created inside a rectangular box of sides 'b' and 'h' and then, one by one, the algorithm explores the intersection between all line segments. This is a Monte Carlo simulation, so the 'experiment' is executed many times (no less than 100 times). Written like that, very much RAM is consumed:  Here, the element Mij=1 if stick i intersects stick j and Mij=0 if not.
>How can I optimize my algorithm? Graph theory is useful in this case?

NumPy-Discussion mailing list
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20130829/ae9a3ba9/attachment-0001.html 

More information about the NumPy-Discussion mailing list