[bf #331] raising exceptions when poly ill defined

[lasofivec]

As discussed in #331 :
The polygon needs to be properly defined in order for the algorithm isClockwise works correctly.
If it finds an issue (null volume) it will raise an exception.

The exception is caught and raised whenever function used
Warning: I didn't manage to properly catch and raise the error. Can you see what I am doing wrong ? I left a print for the moment

I also realized in _core, isClockwise is defined sometimes before checking if the polygon is closed, however the algorithm needs for the polygon to be closed.

Lastly: some flake8 changes

[pep8speaks]

Hello @lasofivec! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2020-01-22 12:51:52 UTC

[lasofivec]

added this since I didnt manage to catch the Exception later

[lasofivec]

first we close the poly then we check if clockwise

[lasofivec]

is not being used, so I commented it

[Didou09]

What does this flag do ?

[lasofivec]

Cython normally checks if a memory view is initialized before using it. This turns the verification off (for speed purposes). I wanted to match the flags with the GG file.

[Didou09]

I think I understand what's going on:

The config creation shows that 4 polygons are not working,
This is one more than before, so the changes we've made so far made the problem worse somehow.

Let's consider the config creation with this fork and re-compiled files:

In [1]: import tofu as tf
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/imas2tofu/init.py:83: UserWarning:
You do not seem to be using the latest IMAS version:
'module list' vs 'module av IMAS' suggests:
- Current version: 3.23.1-4.0.3
- Latest version : 3.25.0-4.4.0
warnings.warn(msg)
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/init.py:95: UserWarning:
The following subpackages are not available:
- tofu.mag
=> see print(tofu.dsub[]) for details.
warnings.warn(msg)

In [2]: conf = tf.geom.utils.create_config('AUG')
In Poly_isClockwise : 
   Found lowest right point at index : 0, of coordinates :1.905, 0.868.
   The two neighboring points are : 22 and 1.
In Poly_isClockwise : 
   Found lowest right point at index : 0, of coordinates :1.905, 0.868.
   The two neighboring points are : 22 and 1.
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/geom/utils.py:838: UserWarning: Could not be loaded: TFG_PFC_ExpAUG_D2dBu1.txt
  warnings.warn(msg)
In Poly_isClockwise : 
   Found lowest right point at index : 0, of coordinates :2.002, 0.708.
   The two neighboring points are : 11 and 1.
In Poly_isClockwise : 
   Found lowest right point at index : 0, of coordinates :2.002, 0.708.
   The two neighboring points are : 11 and 1.
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/geom/utils.py:838: UserWarning: Could not be loaded: TFG_PFC_ExpAUG_D2dBu2.txt
  warnings.warn(msg)
In Poly_isClockwise : 
   Found lowest right point at index : 0, of coordinates :1.095, -0.76.
   The two neighboring points are : 8 and 1.
In Poly_isClockwise : 
   Found lowest right point at index : 0, of coordinates :1.095, -0.76.
   The two neighboring points are : 8 and 1.
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/geom/utils.py:838: UserWarning: Could not be loaded: TFG_PFC_ExpAUG_LIM09.txt
  warnings.warn(msg)
In Poly_isClockwise : 
   Found lowest right point at index : 0, of coordinates :1.1, -0.584.
   The two neighboring points are : 42 and 1.
In Poly_isClockwise : 
   Found lowest right point at index : 0, of coordinates :1.1, -0.584.
   The two neighboring points are : 42 and 1.
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/geom/utils.py:838: UserWarning: Could not be loaded: TFG_PFC_ExpAUG_SBi.txt
  warnings.warn(msg)

These 4 polygons have something in common: in all cases, the index of the lower-left corner in 0.
This means that the two neighbours are indices 1 and -1 (the last point).
But these polygons are all closed, so point -1 is identical to point 0.

Hence the algorithm fails, because it is applied on 3 points, but 2 of them are the same.

We can check that, contrary to what I thought, the polygons are not corrupted, and that the first point (index 0) i always the lowest point:

In [3]: import numpy as np
In [4]: import matplotlib.pyplot as plt

In [5]: p0 = np.loadtxt('tofu/geom/inputs/TFG_PFC_ExpAUG_D2dBu1.txt')[1:, :]

In [6]: p1 = np.loadtxt('tofu/geom/inputs/TFG_PFC_ExpAUG_D2dBu2.txt')[1:, :]

In [7]: p2 = np.loadtxt('tofu/geom/inputs/TFG_PFC_ExpAUG_LIM09.txt')[1:, :]

In [8]: p3 = np.loadtxt('tofu/geom/inputs/TFG_PFC_ExpAUG_SBi.txt')[1:, :]

In [9]: plt.figure()

In [10]: [plt.plot(pp[:, 0], pp[:, 1], '-k', marker='.', ms=4) for pp in [p3, p2, p1, p0]]

In [11]: [plt.plot(pp[0, 0], pp[0, 1], 'or', ms=4) for pp in [p3, p2, p1, p0]]

In [12]: plt.gca().set_aspect('equal')

Regarding the possibly corrupted polygon on the figure you showed in issue #331 , it is just a display effect: the index of each segment is printed on the outside of the polygon, at a certain distance from it.
Here the polygon is so thin that it looks like the two sides are inverted but in fact the first segments are on the right and the last ones are on the left, as we can see by printing the corresponding polygon and first point in red again and zooming:

I see 2 possible solutions:

• Either make sure the polygons are not closed when passed to Poly_isClockwise()
• Or make sure Poly_isClockwise() properly handles cases when the lowest point index is 0 (by setting one neighbourg to 1 and the other to -2 instead of -1)

Also, it seems that some of the correction you proposed, although relevant (make sure all polygons are closed / anti-clockwise) introduced another error: some of them are not C-contiguous anymore.

• Either we go back and remove these corrections (but loose the benefit)
• or we make sure all modified polygons are C-contiguous when necessary

What do you think ?

[lasofivec]

Ok, my bad. Forgot that it needed to be a "-2" instead of "-1".
This solved the issue, event hough I am not entirely sure why did this also solve the issue with the plotting of the structure ?

[Didou09]

It's because the plotting procedure also depends on the order of the polygon (clockwise or not).
The position of the labels (indices) is determined by the in/out vectors perpendicular to each segment, and these in/out vectors may be wrong if there was a bug in one of the preceding steps.

[Didou09]

I'm testing it locally, but it seems we are fine,
I just want to double-check because of the C-contiguous error I found last time

[Didou09]

So, locally it seems we are good,
But travis finds an error in the unit tests:

======================================================================
ERROR: tofu.tests.tests01_geom.tests01_GG.test02_Poly_CLockOrder
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/travis/miniconda/lib/python3.6/site-packages/nose/case.py", line 198, in runTest
    self.test(*self.arg)
  File "/home/travis/build/ToFuProject/tofu/tofu/tests/tests01_geom/tests01_GG.py", line 119, in test02_Poly_CLockOrder
    not GG.Poly_isClockwise(P), P.flags['C_CONTIGUOUS'],
  File "tofu/geom/_GG.pyx", line 223, in tofu.geom._GG.Poly_isClockwise
Exception: In Poly_isClockwise : 
   Found lowest right point at index : 1, of coordinates :0.0, 0.0.
   The two neighboring points are : 0 and 0.

Could you have a look ?

[lasofivec]

I am hunting this bug :)

[lasofivec]

Ok it took me a while to understand the bug but it is not really a bug.....
This is what happens:

• A polygon is created in the right order: (number_of_dims, number_of_points) = (2, npts)
• It is passed to the Poly_Order function with the parameter order=(N,cc) meaning that the polygon is going to be translated so that its shape is (npts, 2)
• Then we check if the polygon is not clockwise, however in this algorithm we suppose that the polygon is in the order (2, npts).
• Original version of the algorithm (coded by you, see below) and the version in devel, the computation was done and because the algorithm will return False since it's being computed on a segment instead of the polygon.
• Now we raise an exception if it finds that the volume is 0 (which is the case here), thus the assert assert not GG.Poly_isClockwise(P) fails...

So, I would suggest either:

• Taking out the assert in the test function and maybe adding a test to check if the number of points in GG.Poly_isClockwise is > 2.
• Checking if we find npts < 3, and in that case translate the polygon, and let the algorithm run as is. If for example we send a polygon with shape (2,2), it will still fail (I can add in the exception message the shape of the polygon received).

I wouldn't suggest calling GG.Poly_Order in GG.Poly_isClockwise since GG.Poly_Order already calls GG.Poly_isClockwise.

FTR, original version of the function

@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
def Poly_isClockwise(cnp.ndarray[double,ndim=2] Poly):
    """ Assuming 2D closed Poly ! """
    cdef int ii, NP=Poly.shape[1]
    cdef double Sum=0.
    for ii in range(0,NP-1):
        Sum += Poly[0,ii]*Poly[1,ii+1]-Poly[0,ii+1]*Poly[1,ii]
    return Sum < 0.

[lasofivec]

I think the second solution is the most elegant, so I made changes in that direction. Let me know what you think

[Didou09]

Hum, I think I was originally translating from (2, npts) to (npts, 2) because before I was using external library Polygon2 (or Polygon3) to computed some properies (area...).
Since your optimizations, we got rid of this dependency, so we can probably stop translating from (2, npts) to (npts, 2) and stick with (2, npts) all the time.

hat should solve the issue, wouldn't it ?

[Didou09]

Also, I agree your solution 2 is the most elegant

[lasofivec]

Yes, but this would mean we should make some major changes to PolyOrder:

• Take out the option of changing the order
• Renaming the function :)
• Checking all tests and code that call it

[lasofivec]

Probably we can open a new issue about this if you prefer

[Didou09]

Hum, good point,

What do you think ?
I'll let you decide which option you prefer, either a hotfix (your option 2) or a full-blown rewriting of the incriminated parts.

Given that we are a bit late on time I would probably go for a hotfix, but I understand if you prefer something cleaner and more complete. your move ;-)

[lasofivec]

we can go for the hotfix (in this branch and later in devel) and open an issue for the cleaner solution.

[Didou09]

OK for me, this is the hotfix right ?
If so we can merge, I guess we can also delete the branch ?