Issue281 new default configs - AUG

Notes_Upgrades/meetings/nov_2019.md

@@ -32,7 +32,7 @@ Summary of meeting 27/11/2019
 - Documentation problems:
 
     - More and more users asking for doc (mostly on IMAS).
-    - docstring: 
+    - docstring:
       Didier will update docstring of main user functions
       **Assignee**: Didier
   - website:
@@ -41,7 +41,7 @@ Summary of meeting 27/11/2019
     - add tutorial for reflexions (number and types)
     - 2d radiations with ITER data (best case) or fake data -> update the tutorial called Computing a camera image with custom emissivity
     - sinogram: to be added to an existing tutorial
-    - how to build cameras: add more comments on how to create a camera (in the 5 minutes to Tofu) what are the different parameters. Maybe separate camera 1D and camera 2D and go into details on how it works. 
+    - how to build cameras: add more comments on how to create a camera (in the 5 minutes to Tofu) what are the different parameters. Maybe separate camera 1D and camera 2D and go into details on how it works.
 
     **Assignee**:  Laura + Didier + Florian ?
     Deadline: December ?

tofu/geom/_GG.pyx

@@ -188,30 +188,51 @@ 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)
     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 - 1
+        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])
+    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) + ".")
+        raise Exception(err_msg) # not working
     return res < 0.
 
 
@@ -274,8 +295,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)':

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 :

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

tofu/geom/_comp.py

@@ -34,9 +34,11 @@ 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
 
@@ -61,9 +63,13 @@ def _Struct_set_Poly(
         Vol, BaryV = None, None
     else:
         Vol, BaryV = _GG.Poly_VolAngTor(Poly)
-        msg = "Pb. with volume computation for Ves object of type 'Tor' !"
-        msg = "\n Here Volume = " + str(Vol)
-        assert Vol > 0.0, msg
+        if Vol <= 0.0:
+            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
     Vect = np.diff(Poly, n=1, axis=1)

tofu/geom/_core.py

@@ -436,19 +436,23 @@ def _checkformat_inputs_dgeom(
             Poly = Poly.T
 
         # Elimininate any double identical point
-        ind = np.sum(np.diff(Poly, axis=1) ** 2, axis=0) < 1.0e-12
+        ind = np.sum(np.diff(np.concatenate((Poly, Poly[:, 0:1]), axis=1),
+                             axis=1) ** 2, axis=0) < 1.0e-12
         if np.any(ind):
             npts = Poly.shape[1]
+            Poly = Poly[:, ~ind]
             msg = (
                 "%s instance: double identical points in Poly\n" % cls.__name__
             )
             msg += "  => %s points removed\n" % ind.sum()
             msg += "  => Poly goes from %s to %s points" % (
                 npts,
-                npts - ind.sum(),
+                Poly.shape[1],
             )
             warnings.warn(msg)
-            Poly = Poly[:, ~ind]
+            ind = np.sum(np.diff(np.concatenate((Poly, Poly[:, 0:1]), axis=1),
+                                 axis=1) ** 2, axis=0) < 1.0e-12
+            assert not np.any(ind), ind
 
         lC = [Lim is None, pos is None]
         if not any(lC):
@@ -2199,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
@@ -3130,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
@@ -3724,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"
@@ -4137,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
 
@@ -5029,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']
@@ -5123,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
@@ -5570,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
@@ -5589,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:
@@ -5659,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]])
@@ -6896,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