Repeated call of sTriangulation freezes

[kitchenknif]

Running this code seems to randomly hang after ~3000 iterations, and I have been unable to deduce why.
Using quadrature points from here.

import time
import numpy as np
import stripy

def get_lebedev_quadrature(order=29):
    phi, theta, weights = [], [], []
    with open("lebedev_{:03d}.txt".format(order), "r") as f:
        lines = f.readlines()
        for line in lines:
            ph, th, w = line.split()
            phi.append(float(ph) * np.pi / 180 + np.pi)
            theta.append(float(th) * np.pi / 180)
            weights.append(float(w))

    return np.asarray(phi), np.asarray(theta), np.asarray(weights)

if __name__ == "__main__":
    phi, theta, weights = get_lebedev_quadrature(29)
    for i in range(10000):
        a = time.time()
        tri = stripy.sTriangulation(lons=phi - np.pi, lats=-theta + np.pi / 2, permute=True)
        b = time.time() - a
        print(i, b)

lebedev_029.txt

[lmoresi]

My guess would be that you are seeing some failure in garbage collection which is occurring when f2py is used to wrap the underlying fortran code. It’s something of a black box and I suspect we will need to try a few things to figure out what is going on. I wonder if the behaviour would be more predictable if you explicitly delete the objects in the loop. We haven’t really battle-tested this use case - normally we build a small number of triangulations and keep them around for a while. However, we have recommended people re-triangulate rather than adding functionality to add / remove points so even our standard use-case is vulnerable to this problem. L On 28 Sep 2020, 8:23 PM +1000, Pavel Dmitriev <notifications@github.com>, wrote: Running this code seems to randomly hang after +- 3000 iterations, and I have been unable to deduce why. Using quadrature points from here<https://people.sc.fsu.edu/~jburkardt/datasets/sphere_lebedev_rule/sphere_lebedev_rule.html>. import time import numpy as np import stripy def get_lebedev_quadrature(order=29): phi, theta, weights = [], [], [] with open("lebedev_{:03d}.txt".format(order), "r") as f: lines = f.readlines() for line in lines: ph, th, w = line.split() phi.append(float(ph) * np.pi / 180 + np.pi) theta.append(float(th) * np.pi / 180) weights.append(float(w)) return np.asarray(phi), np.asarray(theta), np.asarray(weights) if __name__ == "__main__": phi, theta, weights = get_lebedev_quadrature(29) for i in range(10000): a = time.time() tri = stripy.sTriangulation(lons=phi - np.pi, lats=-theta + np.pi / 2, permute=True) b = time.time() - a print(i, b) lebedev_029.txt<https://github.com/underworldcode/stripy/files/5291519/lebedev_029.txt> — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#85>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADABPI4HQRMPGN37Q4MWDVDSIBPYXANCNFSM4R4KNCYQ>.

[kitchenknif]

I've been trying to forcibly run garbage collection at the beginning / end of the loop, but that doesn't seem to have any effect.

Explicitly deleting the triangulation after use seems to postpone the hang by a few thousand iterations, but it still ends up hanging.

[lmoresi]

I can confirm that I am also seeing this issue. At / around the time the code fails I get this error message: TRMESH - Fatal error! The first 3 nodes are collinear. Try reordering the data. This is not the same error that we get if we do not permute the mesh at all. I wonder if this is a case of the permutation loop not doing the right thing if it exhausts its attempts at permutation @brmather ? L On 28 Sep 2020, 8:42 PM +1000, Pavel Dmitriev <notifications@github.com>, wrote: I've been trying to forcibly run garbage collection at the beginning / end of the loop, but that doesn't seem to have any effect. — You are receiving this because you commented. Reply to this email directly, view it on GitHub<#85 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADABPI4IWTCPYCGJMJ3XKALSIBR77ANCNFSM4R4KNCYQ>.

[kitchenknif]

Yeah, that's another Issue that I have,
and it is weird that it only occurs sometimes, considering that the input points are exactly the same during each call to build the triangulation.

[lmoresi]

The permutation is random, and, of course, it is possible to generate another ordering that has co-circular points. It has a loop to keep trying but maybe we don’t try hard enough or test for failure correctly. The problem is easily reproducible from the code you sent so we can dig a little more. L On 28 Sep 2020, 9:35 PM +1000, Pavel Dmitriev <notifications@github.com>, wrote: Yeah, that's another Issue that I have, and it is weird that it only occurs sometimes, considering that the input points are exactly the same during each call to build the triangulation. — You are receiving this because you commented. Reply to this email directly, view it on GitHub<#85 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADABPI7QJK5YSMWAA7ZKBDLSIBYIXANCNFSM4R4KNCYQ>.

[brmather]

The "TRMESH - Fatal error!" message is printed at the fortran level if a permutation fails to produce non-collinear points in the first 3 entries of the lon / lat arrays, but stripy simply reshuffles them again before handing it back. In the past I have written initial checks for collinear points prior to constructing the mesh, but they occasionally terminated at the fortran TRMESH routine. It may just be floating point precision causing that particular issue.

I was able to reproduce the bug for a single case. Download this permutation.txt file and execute the following code:

lons, lats, p, ip = np.loadtxt("permutation.txt", unpack=True)
p  = p.astype(np.int)
ip = ip.astype(np.int)

lons0 = lons[p]
lats0 = lats[p]

tri = stripy.sTriangulation(lons0, lats0, permute=False)

🕥

@lmoresi - I'm not sure what is wrong with this permutation for it to hang. The first 3 points are not collinear and they seem adequately spaced apart...

def is_left(x1,y1,z1, x2,y2,z2, x0,y0,z0):
    left = x0 * ( y1 * z2 - y2 * z1 ) \
         - y0 * ( x1 * z2 - x2 * z1 ) \
         + z0 * ( x1 * y2 - x2 * y1 ) >= 0.0
    return left

def is_collinear(lons, lats):
    from stripy.spherical import lonlat2xyz
    x, y, z = lonlat2xyz(lons[:3], lats[:3])
    
    if not is_left(x[0],y[0],z[0], x[1],y[1],z[1], x[2],y[2],z[2]):
        collinear = False
    elif not is_left(x[1],y[1],z[1], x[0],y[0],z[0], x[2],y[2],z[2]):
        collinear = False
    else:
        collinear = True
        
    cartesian_coords = np.c_[x,y,z]
    
    A = np.linalg.norm(cartesian_coords[1] - cartesian_coords[0])
    B = np.linalg.norm(cartesian_coords[2] - cartesian_coords[0])
    C = np.linalg.norm(cartesian_coords[2] - cartesian_coords[1])
    
    s = 0.5*(A+B+C)
    area = np.sqrt(s*(s-A)*(s-B)*(s-C))
    
    det = np.linalg.det(cartesian_coords.T)
    if np.isclose(det, 0.0):
        collinear = True

    print(area)

    return collinear

is_collinear(lons, lats)

> 0.10537979877033213
> False

[lmoresi]

Another data point: 1) I permuted the points by hand in the text file and turned off the permute option. One million non-hanging triangulations later … not even any memory leaks worth mentioning. 2) I used my permuted text file and turned permutation on again and sure enough: 0 50 100 150 200 Permutation failed, trying again Permutation failed, trying again 250 300 350 400 450 500 550 600 Permutation failed, trying again 650 700 750 800 and it hangs pretty quickly after that. On 29 Sep 2020, 3:14 PM +1000, Ben Mather <notifications@github.com>, wrote: The "TRMESH - Fatal error!" message is printed at the fortran level if a permutation fails to produce non-collinear points in the first 3 entries of the lon / lat arrays, but stripy simply reshuffles them again before handing it back. In the past I have written initial checks for collinear points prior to constructing the mesh, but they occasionally terminated at the fortran TRMESH routine. It may just be floating point precision causing that particular issue. I was able to reproduce the bug for a single case. Download this permutation.txt<https://github.com/underworldcode/stripy/files/5296438/permutation.txt> file and execute the following code: lons, lats, p, ip = np.loadtxt("permutation.txt", unpack=True) p = p.astype(np.int) ip = ip.astype(np.int) lons0 = lons[p] lats0 = lats[p] tri = stripy.sTriangulation(lons0, lats0, permute=False) 🕥 @lmoresi<https://github.com/lmoresi> - I'm not sure what is wrong with this permutation for it to hang. The first 3 points are not collinear and they seem adequately spaced apart... def is_left(x1,y1,z1, x2,y2,z2, x0,y0,z0): left = x0 * ( y1 * z2 - y2 * z1 ) \ - y0 * ( x1 * z2 - x2 * z1 ) \ + z0 * ( x1 * y2 - x2 * y1 ) >= 0.0 return left def is_collinear(lons, lats): from stripy.spherical import lonlat2xyz x, y, z = lonlat2xyz(lons[:3], lats[:3]) if not is_left(x[0],y[0],z[0], x[1],y[1],z[1], x[2],y[2],z[2]): collinear = False elif not is_left(x[1],y[1],z[1], x[0],y[0],z[0], x[2],y[2],z[2]): collinear = False else: collinear = True cartesian_coords = np.c_[x,y,z] A = np.linalg.norm(cartesian_coords[1] - cartesian_coords[0]) B = np.linalg.norm(cartesian_coords[2] - cartesian_coords[0]) C = np.linalg.norm(cartesian_coords[2] - cartesian_coords[1]) s = 0.5*(A+B+C) area = np.sqrt(s*(s-A)*(s-B)*(s-C)) det = np.linalg.det(cartesian_coords.T) if np.isclose(det, 0.0): collinear = True print(area) return collinear is_collinear(lons, lats)

0.10537979877033213 False

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#85 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADABPI3VZ6WSKUVRTW553UTSIFULHANCNFSM4R4KNCYQ>.

[brmather]

@kitchenknif for now you can fix your code by manually permuting the first 3 points in your list. Adding this line:

phi[1], phi[4] = phi[4], phi[1]
theta[1], theta[4] = theta[4], theta[1]

and setting permute=False was enough to get your loop working.

We will keep troubleshooting to determine the root cause of the problem.

[kitchenknif]

My current solution looks like this, manually permuting the points if something goes wrong.

    theta = np.roll(theta, 10)
    phi = np.roll(phi, 10)
    weights = np.roll(weights, 10)

    tri = None
    while tri is None:
        try:
            tri = stripy.sTriangulation(lons=phi - np.pi, lats=-theta + np.pi / 2, permute=False)
        except ValueError as e:
            print(e)
            phi = np.roll(phi, 1)
            theta = np.roll(theta, 1)
            weights = np.roll(weights, 1)