Issue188 calc in out isoflux more flexible

_updateversion.py

@@ -23,9 +23,10 @@ def updateversion(path=_HERE):
             version_git = fh.read().strip().split("=")[-1].replace("'",'')
     version_git = version_git.lower().replace('v','')
 
-    version_msg = "# Do not edit this file, pipeline versioning is governed by git tags !"
+    version_msg = "# Do not edit, pipeline versioning governed by git tags!"
     with open(version_py,"w") as fh:
-        fh.write((version_msg + os.linesep + '__version__=').replace('',"") + "'%s'" % version_git)
+        msg = "{0}__version__ = '{1}'{0}".format(os.linesep, version_git)
+        fh.write(version_msg + msg)
     return version_git
 
 

tofu/geom/_core.py

@@ -4639,7 +4639,7 @@ def get_sample(self, res=None, resMode='abs', DL=None, method='sum', ind=None,
             k = np.split(k, lind, axis=-1)
         return k, reseff, lind
 
-    def _kInOut_IsoFlux_inputs(self, lPoly, lVIn=None):
+    def _kInOut_Isoflux_inputs(self, lPoly, lVIn=None):
 
         if self._method=='ref':
             D, u = np.ascontiguousarray(self.D), np.ascontiguousarray(self.u)
@@ -4668,52 +4668,101 @@ def _kInOut_IsoFlux_inputs(self, lPoly, lVIn=None):
             pass
         return largs, dkwd
 
-    def _kInOut_IsoFlux_inputs_usr(self, lPoly, lVIn=None):
+    def _kInOut_Isoflux_inputs_usr(self, lPoly, lVIn=None):
+        c0 = type(lPoly) in [np.ndarray, list, tuple]
 
         # Check lPoly
-        if type(lPoly) is np.ndarray:
-            lPoly = [lPoly]
-        lPoly = [np.ascontiguousarray(pp) for pp in lPoly]
-        msg = "Arg lPoly must be a list of (2,N) or (N,2) np.ndarrays !"
-        assert all([pp.ndim==2 and 2 in pp.shape for pp in lPoly]), msg
+        if c0 and type(lPoly) is np.ndarray:
+            c0 = c0 and lPoly.ndim in [2, 3]
+            if c0 and lPoly.ndim == 2:
+                c0 = c0 and lPoly.shape[0] == 2
+                if c0:
+                    lPoly = [np.ascontiguousarray(lPoly)]
+            elif c0:
+                c0 = c0 and lPoly.shape[1] == 2
+                if c0:
+                    lPoly = np.ascontiguousarray(lPoly)
+        elif c0:
+            lPoly = [np.ascontiguousarray(pp) for pp in lPoly]
+            c0 = all([pp.ndim == 2 and pp.shape[0] == 2 for pp in lPoly])
+        if not c0:
+            msg = "Arg lPoly must be either:\n"
+            msg += "    - a (2,N) np.ndarray (signle polygon of N points)\n"
+            msg += "    - a list of M polygons, each a (2,Ni) np.ndarray\n"
+            msg += "        - where Ni is the number of pts of each polygon\n"
+            msg += "    - a (M,2,N) np.ndarray where:\n"
+            msg += "        - M is the number of polygons\n"
+            msg += "        - N is the (common) number of points per polygon\n"
+            raise Exception(msg)
         nPoly = len(lPoly)
-        for ii in range(0,nPoly):
-            if lPoly[ii].shape[0]!=2:
-                lPoly[ii] = lPoly[ii].T
+
+        # Check anti-clockwise and closed
+        if type(lPoly) is list:
+            for ii in range(nPoly):
                 # Check closed and anti-clockwise
-                lPoly[ii] = _GG.Poly_Order(lPoly[ii], Clock=False, close=True)
+                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
+                    )
+        else:
+            for ii in range(nPoly):
+                # Check closed and anti-clockwise
+                if _GG.Poly_isClockwise(lPoly[ii]):
+                    lPoly[ii] = lPoly[ii][:, ::-1]
+            d = np.sum((lPoly[:, :, 0]-lPoly[:, :, -1])**2, axis=1)
+            if np.allclose(d, 0.):
+                pass
+            elif np.all(d > 0.):
+                lPoly = np.concatenate((lPoly, lPoly[:, :, 0:1]), axis=-1)
+            else:
+                msg = "All poly in lPoly should be closed or all non-closed!"
+                raise Exception(msg)
+
 
         # Check lVIn
         if lVIn is None:
             lVIn = []
             for pp in lPoly:
-                VIn = np.diff(pp, axis=1)
-                VIn = VIn/(np.sqrt(np.sum(VIn**2,axis=0))[np.newaxis,:])
-                VIn = np.ascontiguousarray([-VIn[1,:],VIn[0,:]])
-                lVIn.append(VIn)
+                vIn = np.diff(pp, axis=1)
+                vIn = vIn/(np.sqrt(np.sum(vIn**2, axis=0))[None, :])
+                vIn = np.ascontiguousarray([-vIn[1, :], vIn[0, :]])
+                lVIn.append(vIn)
         else:
-            if type(lVIn) is np.ndarray:
-                lVIn = [lVIn]
-            assert len(lVIn)==nPoly
-            lVIn = [np.ascontiguousarray(pp) for pp in lVIn]
-            msg = "Arg lVIn must be a list of (2,N) or (N,2) np.ndarrays !"
-            assert all([pp.ndim==2 and 2 in pp.shape for pp in lVIn]), msg
+            c0 = type(lVIn) in [np.ndarray, list, tuple]
+            if c0 and type(lVIn) is np.ndarray and lVIn.ndim == 2:
+                c0 = c0 and lVIn.shape == (2, lPoly[0].shape[1]-1)
+                if c0:
+                    lVIn = [np.ascontiguousarray(lVIn)]
+            elif c0 and type(lVIn) is np.ndarray:
+                c0 = c0 and lVIn.shape == (nPoly, 2, lPoly.shape[-1]-1)
+                if c0:
+                    lVIn = np.ascontiguousarray(lVIn)
+            elif c0:
+                c0 = c0 and len(lVIn) == nPoly
+                if c0:
+                    c0 = c0 and all([vv.shape == (2, pp.shape[1]-1)
+                                     for vv, pp in zip(lVIn, lPoly)])
+                    if c0:
+                        lVIn = [np.ascontiguousarray(vv) for vv in lVIn]
+
+            # Check normalization and direction
             for ii in range(0,nPoly):
-                if lVIn[ii].shape[0]!=2:
-                    lVIn[ii] = lVIn[ii].T
-                    lVIn[ii] = lVIn[ii]/(np.sqrt(np.sum(lVIn[ii]**2,axis=0))[np.newaxis,:])
-                    assert lVIn[ii].shape==(2,lPoly[ii].shape[1]-1)
-                    vect = np.diff(lPoly[ii],axis=1)
-                    det = vect[0,:]*lVIn[ii][1,:] - vect[1,:]*lVIn[ii][0,:]
-                    if not np.allclose(np.abs(det),1.):
-                        msg = "Each lVIn must be perp. to each lPoly segment !"
-                        raise Exception(msg)
-                    ind = np.abs(det+1)<1.e-12
-                    lVIn[ii][:,ind] = -lVIn[ii][:,ind]
+                lVIn[ii] = (lVIn[ii]
+                            / np.sqrt(np.sum(lVIn[ii]**2, axis=0))[None, :])
+                vect = np.diff(lPoly[ii], axis=1)
+                vect = vect / np.sqrt(np.sum(vect**2, axis=0))[None, :]
+                det = vect[0, :]*lVIn[ii][1, :] - vect[1, :]*lVIn[ii][0, :]
+                if not np.allclose(np.abs(det), 1.):
+                    msg = "Each lVIn must be perp. to each lPoly segment !"
+                    raise Exception(msg)
+                ind = np.abs(det+1) < 1.e-12
+                lVIn[ii][:, ind] = -lVIn[ii][:, ind]
 
         return nPoly, lPoly, lVIn
 
-    def calc_kInkOut_IsoFlux(self, lPoly, lVIn=None, Lim=None,
+    def calc_kInkOut_Isoflux(self, lPoly, lVIn=None, Lim=None,
                              kInOut=True):
         """ Calculate the intersection points of each ray with each isoflux
 
@@ -4733,34 +4782,33 @@ def calc_kInkOut_IsoFlux(self, lPoly, lVIn=None, Lim=None,
         """
 
         # Preformat input
-        nPoly, lPoly, lVIn = self._kInOut_IsoFlux_inputs_usr(lPoly, lVIn=lVIn)
+        nPoly, lPoly, lVIn = self._kInOut_Isoflux_inputs_usr(lPoly, lVIn=lVIn)
 
         # Prepare output
-        kIn = np.full((self.nRays,nPoly), np.nan)
-        kOut = np.full((self.nRays,nPoly), np.nan)
+        kIn = np.full((nPoly, self.nRays), np.nan)
+        kOut = np.full((nPoly, self.nRays), np.nan)
 
         # Compute intersections
         assert(self._method in ['ref', 'optimized'])
         if self._method=='ref':
             for ii in range(0,nPoly):
-                largs, dkwd = self._kInOut_IsoFlux_inputs([lPoly[ii]],
+                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]
+                kIn[ii, :], kOut[ii, :] = out[2], out[3]
         elif self._method=="optimized":
             for ii in range(0,nPoly):
-                largs, dkwd = self._kInOut_IsoFlux_inputs([lPoly[ii]],
+                largs, dkwd = self._kInOut_Isoflux_inputs([lPoly[ii]],
                                                           lVIn=[lVIn[ii]])
 
-                out = _GG.LOS_Calc_PInOut_VesStruct(*largs, **dkwd)
-                kin, kout, _, _ = out
-                kIn[:,ii], kOut[:,ii] = kin, kout
+                out = _GG.LOS_Calc_PInOut_VesStruct(*largs, **dkwd)[:2]
+                kIn[ii, :], kOut[ii, :] = out
         if kInOut:
             indok = ~np.isnan(kIn)
-            ind = np.zeros((self.nRays,nPoly), dtype=bool)
-            kInref = np.tile(self.kIn[:,np.newaxis],nPoly)
-            kOutref = np.tile(self.kOut[:,np.newaxis],nPoly)
+            ind = np.zeros((nPoly, self.nRays), dtype=bool)
+            kInref = np.tile(self.kIn, (nPoly, 1))
+            kOutref = np.tile(self.kOut, (nPoly, 1))
             ind[indok] = (kIn[indok]<kInref[indok]) | (kIn[indok]>kOutref[indok])
             kIn[ind] = np.nan
 

tofu/tests/tests01_geom/tests03_core.py

@@ -791,8 +791,7 @@ def test08_get_sample(self):
                         assert np.all((k[ii][ind]>=obj.kIn[ii]-res[ii])
                                       & (k[ii][ind]<=obj.kOut[ii]+res[ii]))
 
-
-    def test09_calc_kInkOut_IsoFlux(self):
+    def test09_calc_kInkOut_Isoflux(self):
         nP = 10
         r = np.linspace(0.1,0.4,nP)
         theta = np.linspace(0.,2*np.pi,100)
@@ -801,33 +800,30 @@ def test09_calc_kInkOut_IsoFlux(self):
         for typ in self.dobj.keys():
             for c in self.dobj[typ].keys():
                 obj = self.dobj[typ][c]
-                kIn, kOut = obj.calc_kInkOut_IsoFlux(lp2D)
-                assert kIn.shape==(obj.nRays, nP)
-                assert kOut.shape==(obj.nRays, nP)
-                for ii in range(0,nP):
-                    ind = ~np.isnan(kIn[:,ii])
-                    if not np.all((kIn[ind,ii]>=obj.kIn[ind])
-                                  & (kIn[ind,ii]<=obj.kOut[ind])):
+                kIn, kOut = obj.calc_kInkOut_Isoflux(lp2D)
+                assert kIn.shape == (nP, obj.nRays)
+                assert kOut.shape == (nP, obj.nRays)
+                for ii in range(0, nP):
+                    ind = ~np.isnan(kIn[ii, :])
+                    if not np.all((kIn[ii, ind] >= obj.kIn[ind])
+                                  & (kIn[ii, ind] <= obj.kOut[ind])):
                         msg = typ+' '+c+' '+str(ii)
-                        msg += "\n {0} {1}".format(obj.kIn[ind],obj.kOut[ind])
-                        msg += "\n {0}".format(str(kIn[ind,ii]))
-                        print(msg)
+                        msg += "\n {0} {1}".format(obj.kIn[ind], obj.kOut[ind])
+                        msg += "\n {0}".format(str(kIn[ii, ind]))
                         raise Exception(msg)
 
-                    ind = ~np.isnan(kOut[:,ii])
-                    if not np.all((kOut[ind,ii]>=obj.kIn[ind])
-                                  & (kOut[ind,ii]<=obj.kOut[ind])):
+                    ind = ~np.isnan(kOut[ii, :])
+                    if not np.all((kOut[ii, ind] >= obj.kIn[ind])
+                                  & (kOut[ii, ind] <= obj.kOut[ind])):
                         msg = typ+' '+c+' '+str(ii)
-                        msg += "\n {0} {1}".format(obj.kIn[ind],obj.kOut[ind])
-                        msg += "\n {0}".format(str(kOut[ind,ii]))
-                        print(msg)
+                        msg += "\n {0} {1}".format(obj.kIn[ind], obj.kOut[ind])
+                        msg += "\n {0}".format(str(kOut[ii, ind]))
                         raise Exception(msg)
-                    ind = (~np.isnan(kIn[:,ii])) & (~np.isnan(kOut[:,ii]))
-                    if not np.all(kIn[ind,ii]<=kOut[ind,ii]):
+                    ind = (~np.isnan(kIn[ii, :])) & (~np.isnan(kOut[ii, :]))
+                    if not np.all(kIn[ii, ind] <= kOut[ii, ind]):
                         msg = typ+' '+c+' '+str(ii)
-                        msg += "\n {0}".format(str(kIn[ind,ii]))
-                        msg += "\n {0}".format(str(kOut[ind,ii]))
-                        print(msg)
+                        msg += "\n {0}".format(str(kIn[ii, ind]))
+                        msg += "\n {0}".format(str(kOut[ii, ind]))
                         raise Exception(msg)
 
 

tofu/version.py

@@ -1,2 +1,2 @@
-# Do not edit this file, pipeline versioning is governed by git tags !
-__version__='1.4.1-4-g7cfbe9f'
+# Do not edit, pipeline versioning governed by git tags!
+__version__ = '1.4.1-49-g373fc3c'