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

[Brett Olsen]

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

    condition_1=p1*p2
    condition_2=p3*p4

    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.

~Brett

On Fri, Aug 23, 2013 at 5:09 PM, Josè Luis Mietta <
joseluismietta at 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?
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20130826/4c5d6ba5/attachment.html>