From ebd494b5949c16d562eccfda7dd9c4d439aecfe5 Mon Sep 17 00:00:00 2001 From: Laura Mendoza Date: Tue, 21 Jan 2020 20:37:55 +0100 Subject: [PATCH 1/5] [bf #331] raising exceptions when poly ill defined --- tofu/geom/_GG.pyx | 36 +++++++++++++++++++---------- tofu/geom/_basic_geom_tools.pxd | 1 + tofu/geom/_basic_geom_tools.pyx | 2 ++ tofu/geom/_comp.py | 9 +++++--- tofu/geom/_core.py | 40 +++++++++++++++++++-------------- 5 files changed, 56 insertions(+), 32 deletions(-) diff --git a/tofu/geom/_GG.pyx b/tofu/geom/_GG.pyx index 24dbd4dc2..26e5baeb8 100644 --- a/tofu/geom/_GG.pyx +++ b/tofu/geom/_GG.pyx @@ -188,16 +188,14 @@ def CoordShift(Pts, In='(X,Y,Z)', Out='(R,Z)', CrossRef=None): def Poly_isClockwise(np.ndarray[double,ndim=2] Poly): """ Assuming 2D closed Poly ! - TODO @LM : http://www.faqs.org/faqs/graphics/algorithms-faq/ - A slightly faster method is based on the observation that it isn't - necessary to compute the area. Find the lowest vertex (or, if - there is more than one vertex with the same lowest coordinate, - the rightmost of those vertices) and then take the cross product - of the edges fore and aft of it. Both methods are O(n) for n vertices, - but it does seem a waste to add up the total area when a single cross - product (of just the right edges) suffices. Code for this is - available at ftp://cs.smith.edu/pub/code/polyorient.C (2K). + Find the lowest vertex (or, if there is more than one vertex with + the same lowest coordinate, the rightmost of those vertices) and then + take the cross product of the edges before and after it. + Both methods are O(n) for n vertices, but it does seem a waste to add up + the total area when a single cross product (of just the right edges) + suffices. Code for this is available at + ftp://cs.smith.edu/pub/code/polyorient.C (2K). """ cdef double res cdef double[:,::1] mv_poly = np.ascontiguousarray(Poly) @@ -206,12 +204,23 @@ def Poly_isClockwise(np.ndarray[double,ndim=2] Poly): cdef double[::1] mvy = mv_poly[1,:] cdef int idmin = _bgt.find_ind_lowerright_corner(mvx, mvy, npts) cdef int idm1 = idmin - 1 - cdef int idp1 = (idmin + 1)%npts + cdef int idp1 = (idmin + 1) % npts + cdef str err_msg = "" if idmin == 0 : idm1 = npts - 1 res = mvx[idm1] * (mvy[idmin] - mvy[idp1]) + \ mvx[idmin] * (mvy[idp1] - mvy[idm1]) + \ mvx[idp1] * (mvy[idm1] - mvy[idmin]) + if abs(res) < _VSMALL: + err_msg += ("In Poly_isClockwise : \n" + + " Found lowest right point at index : " + + str(idmin) + + ", of coordinates :" + str(mvx[idmin]) + + ", " + str(mvy[idmin]) + ".\n" + + " The two neighboring points are : " + + str(idm1) + " and " + str(idp1) + ".") + print(err_msg) + raise Exception(err_msg) # not working return res < 0. @@ -274,8 +283,11 @@ def Poly_Order(np.ndarray[double,ndim=2] Poly, str order='C', Clock=False, if not np.allclose(poly[:,0],poly[:,-1], atol=_VSMALL): poly = np.concatenate((poly,poly[:,0:1]),axis=1) if poly.shape[0]==2 and not Clock is None: - if not Clock==Poly_isClockwise(poly): - poly = poly[:,::-1] + try: + if not Clock==Poly_isClockwise(poly): + poly = poly[:,::-1] + except Exception as excp: + raise excp if not close: poly = poly[:,:-1] if layout.lower()=='(n,cc)': diff --git a/tofu/geom/_basic_geom_tools.pxd b/tofu/geom/_basic_geom_tools.pxd index 99c5bb054..d61936605 100644 --- a/tofu/geom/_basic_geom_tools.pxd +++ b/tofu/geom/_basic_geom_tools.pxd @@ -1,6 +1,7 @@ # cython: boundscheck=False # cython: wraparound=False # cython: cdivision=True +# cython: initializedcheck=False # ################################################################################ # Utility functions for basic geometry : diff --git a/tofu/geom/_basic_geom_tools.pyx b/tofu/geom/_basic_geom_tools.pyx index edcb4f502..ba99e1dce 100644 --- a/tofu/geom/_basic_geom_tools.pyx +++ b/tofu/geom/_basic_geom_tools.pyx @@ -1,6 +1,8 @@ # cython: boundscheck=False # cython: wraparound=False # cython: cdivision=True +# cython: initializedcheck=False +# cimport cython from cython.parallel import prange diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index 5725d26ae..a00cb020d 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -34,9 +34,12 @@ def _Struct_set_Poly( """ Compute geometrical attributes of a Struct object """ # Make Poly closed, counter-clockwise, with '(cc,N)' layout and arrayorder - Poly = _GG.Poly_Order( - Poly, order="C", Clock=False, close=True, layout="(cc,N)", Test=True - ) + try: + Poly = _GG.Poly_Order( + Poly, order="C", Clock=False, close=True, layout="(cc,N)", Test=True + ) + except Exception as excp: + print(excp) assert Poly.shape[0] == 2, "Arg Poly must be a 2D polygon !" fPfmt = np.ascontiguousarray if arrayorder == "C" else np.asfortranarray diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 4bacd3755..fc4c2c8bf 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -2203,8 +2203,8 @@ def _get_largs_dsino(): def _checkformat_inputs_Struct(self, struct, err=True): assert issubclass(struct.__class__, Struct) - C0 = struct.Id.Exp==self.Id.Exp - C1 = struct.Id.Type==self.Id.Type + C0 = struct.Id.Exp == self.Id.Exp + C1 = struct.Id.Type == self.Id.Type C2 = struct.Id.Name.isidentifier() C2 = C2 and '_' not in struct.Id.Name msgi = None @@ -3134,7 +3134,6 @@ def plot_phithetaproj_dist(self, refpt=None, ntheta=None, nphi=None, tit=tit, wintit=wintit, invertx=invertx, draw=draw) - def isInside(self, pts, In="(X,Y,Z)", log="any"): """ Return a 2D array of bool @@ -3728,7 +3727,7 @@ def _checkformat_inputs_dgeom(self, dgeom=None): if (isinstance(dgeom, self._dcases[k]['type']) and all([kk in dgeom.keys() # noqa for kk in self._dcases[k]['lk']]))] - if not len(lC)==1: + if not len(lC) == 1: lstr = [v['lk'] for v in self._dcases.values()] msg = "Arg dgeom must be either:\n" msg += " - dict with keys:\n" @@ -4141,7 +4140,8 @@ def _complete_dX12(self, dgeom): k = np.sum(DDb*(u - np.sqrt(sca2)*dgeom['u'][:, 1:]), axis=0) k = k / (1.0-sca2) - if k[0] > 0 and np.allclose(k, k[0], atol=1.e-3, rtol=1.e-6): + if k[0] > 0 and np.allclose(k, k[0], atol=1.e-3, + rtol=1.e-6): pinhole = dgeom['D'][:, 0] + k[0]*u[:, 0] dgeom['pinhole'] = pinhole @@ -5033,8 +5033,8 @@ def D(self): @property def u(self): - if (self._dgeom['u'] is not None - and self._dgeom['u'].shape[1] == self._dgeom['nRays']): + if self._dgeom['u'] is not None \ + and self._dgeom['u'].shape[1] == self._dgeom['nRays']: u = self._dgeom['u'] elif self.isPinhole: u = self._dgeom['pinhole'][:, None] - self._dgeom['D'] @@ -5127,7 +5127,7 @@ def _isLOS(cls): def _check_indch(self, ind, out=int): if ind is not None: ind = np.asarray(ind) - assert ind.ndim==1 + assert ind.ndim == 1 assert ind.dtype in [np.int64, np.bool_, np.long] if ind.dtype == np.bool_: assert ind.size == self.nRays @@ -5574,17 +5574,17 @@ def _kInOut_Isoflux_inputs_usr(self, lPoly, lVIn=None): if type(lPoly) is list: for ii in range(nPoly): # Check closed and anti-clockwise - if _GG.Poly_isClockwise(lPoly[ii]): - lPoly[ii] = lPoly[ii][:, ::-1] if not np.allclose(lPoly[ii][:, 0], lPoly[ii][:, -1]): lPoly[ii] = np.concatenate( (lPoly[ii], lPoly[ii][:, 0:1]), axis=-1 ) + try: + if _GG.Poly_isClockwise(lPoly[ii]): + lPoly[ii] = lPoly[ii][:, ::-1] + except Exception as excp: + print("For structure ", ii, " : ", excp) else: - for ii in range(nPoly): - # Check closed and anti-clockwise - if _GG.Poly_isClockwise(lPoly[ii]): - lPoly[ii] = lPoly[ii][:, ::-1] + # Check closed and anti-clockwise d = np.sum((lPoly[:, :, 0]-lPoly[:, :, -1])**2, axis=1) if np.allclose(d, 0.): pass @@ -5593,6 +5593,12 @@ def _kInOut_Isoflux_inputs_usr(self, lPoly, lVIn=None): else: msg = "All poly in lPoly should be closed or all non-closed!" raise Exception(msg) + for ii in range(nPoly): + try: + if _GG.Poly_isClockwise(lPoly[ii]): + lPoly[ii] = lPoly[ii][:, ::-1] + except Exception as excp: + print("For structure ", ii, " : ", excp) # Check lVIn if lVIn is None: @@ -5663,14 +5669,14 @@ def calc_kInkOut_Isoflux(self, lPoly, lVIn=None, Lim=None, # Compute intersections assert(self._method in ['ref', 'optimized']) - if self._method=='ref': + if self._method == 'ref': for ii in range(0, nPoly): largs, dkwd = self._kInOut_Isoflux_inputs([lPoly[ii]], lVIn=[lVIn[ii]]) out = _GG.SLOW_LOS_Calc_PInOut_VesStruct(*largs, **dkwd) # PIn, POut, kin, kout, VperpIn, vperp, IIn, indout = out[] kIn[ii, :], kOut[ii, :] = out[2], out[3] - elif self._method=="optimized": + elif self._method == "optimized": for ii in range(0, nPoly): largs, dkwd = self._kInOut_Isoflux_inputs([lPoly[ii]], lVIn=[lVIn[ii]]) @@ -6900,7 +6906,7 @@ def get_summary( # Prepare kout = self._dgeom["kOut"] indout = self._dgeom["indout"] - lS = self._dconfig["Config"].lStruct + # lS = self._dconfig["Config"].lStruct angles = np.arccos(-np.sum(self.u*self.dgeom['vperp'], axis=0)) # ar0 From a0f432bb0c53b33a591c1db6be62df7a689d3d9a Mon Sep 17 00:00:00 2001 From: Laura Mendoza Date: Tue, 21 Jan 2020 20:45:31 +0100 Subject: [PATCH 2/5] [bf #331] black changes --- tofu/geom/_comp.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index a00cb020d..667aac2ae 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -35,9 +35,8 @@ def _Struct_set_Poly( # Make Poly closed, counter-clockwise, with '(cc,N)' layout and arrayorder try: - Poly = _GG.Poly_Order( - Poly, order="C", Clock=False, close=True, layout="(cc,N)", Test=True - ) + Poly = _GG.Poly_Order(Poly, order="C", Clock=False, close=True, + layout="(cc,N)", Test=True) except Exception as excp: print(excp) assert Poly.shape[0] == 2, "Arg Poly must be a 2D polygon !" From 9b83143d97bcf1c3af960d4c5b7cd4853655d9d4 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 22 Jan 2020 11:29:35 +0100 Subject: [PATCH 3/5] [Issue281_fork] Sligtly more explicit error message --- tofu/geom/_comp.py | 7 +++++-- tofu/version.py | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index 667aac2ae..a138d34d4 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -64,8 +64,11 @@ def _Struct_set_Poly( else: Vol, BaryV = _GG.Poly_VolAngTor(Poly) if Vol <= 0.0: - msg = "Pb. with volume computation for Ves object of type 'Tor' !" - msg = "\n Here Volume = " + str(Vol) + msg = ("Pb. with volume computation for Struct of type 'Tor' !\n" + + "\t- Vol = {}\n".format(Vol) + + "\t- Poly = {}\n\n".format(str(Poly)) + + " => Probably corrupted polygon\n" + + " => Please check polygon is not self-intersecting") raise Exception(msg) # Compute the non-normalized vector of each side of the Poly diff --git a/tofu/version.py b/tofu/version.py index e2eec41d8..2449f6818 100644 --- a/tofu/version.py +++ b/tofu/version.py @@ -1,2 +1,2 @@ # Do not edit, pipeline versioning governed by git tags! -__version__ = '1.4.2-a5-113-gcbcc34c7' +__version__ = '1.4.2b13-65-ga0f432bb' From 967e02a6ddbcc900db16b19979defbb5483924d0 Mon Sep 17 00:00:00 2001 From: Laura Mendoza Date: Wed, 22 Jan 2020 13:49:33 +0100 Subject: [PATCH 4/5] [bf #331] corrected last point index when lowest right is 0 --- tofu/geom/_GG.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tofu/geom/_GG.pyx b/tofu/geom/_GG.pyx index 26e5baeb8..ae17ef8c5 100644 --- a/tofu/geom/_GG.pyx +++ b/tofu/geom/_GG.pyx @@ -207,7 +207,7 @@ def Poly_isClockwise(np.ndarray[double,ndim=2] Poly): cdef int idp1 = (idmin + 1) % npts cdef str err_msg = "" if idmin == 0 : - idm1 = npts - 1 + idm1 = npts - 2 res = mvx[idm1] * (mvy[idmin] - mvy[idp1]) + \ mvx[idmin] * (mvy[idp1] - mvy[idm1]) + \ mvx[idp1] * (mvy[idm1] - mvy[idmin]) From 8d67f90d284657934882d588c8a4eea715b53ebd Mon Sep 17 00:00:00 2001 From: Laura Mendoza Date: Wed, 22 Jan 2020 16:08:00 +0100 Subject: [PATCH 5/5] [bf #331] checking if poly is in wrong order in isclockwise --- tofu/geom/_GG.pyx | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/tofu/geom/_GG.pyx b/tofu/geom/_GG.pyx index ae17ef8c5..c99f4488d 100644 --- a/tofu/geom/_GG.pyx +++ b/tofu/geom/_GG.pyx @@ -200,14 +200,27 @@ def Poly_isClockwise(np.ndarray[double,ndim=2] Poly): cdef double res cdef double[:,::1] mv_poly = np.ascontiguousarray(Poly) cdef int npts = mv_poly.shape[1] - cdef double[::1] mvx = mv_poly[0,:] - cdef double[::1] mvy = mv_poly[1,:] - cdef int idmin = _bgt.find_ind_lowerright_corner(mvx, mvy, npts) - cdef int idm1 = idmin - 1 - cdef int idp1 = (idmin + 1) % npts + cdef int ndim = mv_poly.shape[0] + cdef double[::1] mvx + cdef double[::1] mvy + cdef int idmin + cdef int idm1 + cdef int idp1 cdef str err_msg = "" + # Checking that Poly wasn't given in the shape (npts, ndim) + if ndim > npts: + mv_poly = np.ascontiguousarray(Poly.T) + npts = mv_poly.shape[1] + ndim = mv_poly.shape[0] + mvx = mv_poly[0,:] + mvy = mv_poly[1,:] + # Getting index of lower right corner and its neighbors + idmin = _bgt.find_ind_lowerright_corner(mvx, mvy, npts) + idm1 = idmin - 1 + idp1 = (idmin + 1) % npts if idmin == 0 : idm1 = npts - 2 + # Computing area of lower right triangle res = mvx[idm1] * (mvy[idmin] - mvy[idp1]) + \ mvx[idmin] * (mvy[idp1] - mvy[idm1]) + \ mvx[idp1] * (mvy[idm1] - mvy[idmin]) @@ -219,7 +232,6 @@ def Poly_isClockwise(np.ndarray[double,ndim=2] Poly): + ", " + str(mvy[idmin]) + ".\n" + " The two neighboring points are : " + str(idm1) + " and " + str(idp1) + ".") - print(err_msg) raise Exception(err_msg) # not working return res < 0.