diff --git a/_updateversion.py b/_updateversion.py index cdcb1f415..b14f878f2 100644 --- a/_updateversion.py +++ b/_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 diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 88d8efa88..7bb2d9014 100644 --- a/tofu/geom/_core.py +++ b/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]kOutref[indok]) kIn[ind] = np.nan diff --git a/tofu/tests/tests01_geom/tests03_core.py b/tofu/tests/tests01_geom/tests03_core.py index 3781e27df..0e0c2ebdb 100644 --- a/tofu/tests/tests01_geom/tests03_core.py +++ b/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) diff --git a/tofu/version.py b/tofu/version.py index 6be347a9d..131f008cd 100644 --- a/tofu/version.py +++ b/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' \ No newline at end of file +# Do not edit, pipeline versioning governed by git tags! +__version__ = '1.4.1-49-g373fc3c'