From 0e57448edd6bc5169a0f73444e763aa0eef8a4fe Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Mon, 26 Aug 2019 12:42:16 +0200 Subject: [PATCH 01/19] [Issue 158] Started Config._reflect_Types() and Config._reflect_geom() and Rays.add_reflections() --- tofu/geom/_core.py | 111 +++++++++++++++++++++++++++++++++++++++++++++ tofu/version.py | 2 +- 2 files changed, 112 insertions(+), 1 deletion(-) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index f36b95f7d..72d8a919b 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -2544,6 +2544,64 @@ def fdistfromwall(self, r, z, phi): return -min(distRZ,distPhi) + # Method handling reflections + + def _reflect_Types(self, indout=None, Type=None): + """ Return an array indicating the Type of reflection for each LOS + + Return a (nRays,) np.ndarray of int indices, each index corresponds to: + - 0: specular reflections + - 1: diffusive reflections + - 2: corner cube reflections + + If indout is provided, the Types are computed according to the + information stored in each corresponding Struct + + If Type is provided, the Type is forced (user-defined) for all LOS + + """ + nRays = u.shape[1] + if Type is not None: + assert Type in ['specular', 'diffusive', 'corner cube'] + Types = np.full((nRays,), _DREFLECT[Type],, dtype=int) + else: + Types = None + return Types + + + def _reflect_geom(self, u, vperp, indout=None, Type=None): + + # Get Types of relection for each Ray + Types = self._reflect_Types(indout=indout, Type=Type) + + # Get indices of each Type + indspec = Types == 0 + inddiff = Types == 1 + indcorn = Types == 2 + + if np.any(np.logical_or(indspec,inddiff)): + sca = np.sum(u*vperp, axis=0) + vpar = np.array([u[1,:]*vperp[2,:] - u[2,:]*vpepr[1,:], + u[2,:]*vperp[0,:] - u[0,:]*vperp[2,:], + u[0,:]*vperp[1,:] - u[1,:]*vperp[0,:]]) + + if np.any(indspec): + # Compute u2 for specular + sca = np.sum(u[:,indspec]*vperp[:,indspec],axis=0,keepdims=True) + sca2 = np.sum(u[:,indspec]*vpar[:,indspec],axis=0,keepdims=True) + u2[:,indspec] = - sca*vperp[:,indspec] + sca2*vpar[:,indspec] + + if np.any(inddiff): + # Compute u2 for diffusive + sca = np.random.random((1,inddiff.sum())) + u2[:,inddiff] = (sca * vperp[:,inddiff] + + np.sqrt(1.-sca**2) * vpar[:,inddiff]) + + if np.any(indcorn): + u2[:,indcorn] = -u[:,indcorn] + return u2, Types + + def plot(self, lax=None, proj='all', element='P', dLeg=_def.TorLegd, indices=False, Lim=None, Nstep=None, @@ -3731,6 +3789,59 @@ def _set_color(self, color=None): def _set_dmisc(self, color=None): self._set_color(color) + + ########### + # Reflections + ########### + + def add_reflections(self, Type=None, nb=None, coefs=None): + """ Add relfected LOS to the camera + + Reflected LOS can be of 3 types: + - 'speculiar' + - 'diffusive' + - 'corer cube' + + """ + + # Check inputs + if Type is not None: + assert Type in ['speculiar', 'diffusive', 'corner cube'] + if nb is None: + nb = 1 + nb = int(nb) + assert nb > 0 + if coefs is not None: + coefs = np.atleast_1d(coefs).astype(float).ravel() + if coefs.size == 1: + coefs = np.repeat(coefs, self.nRays) + assert coefs.shape == (self.nRays,) + assert np.all(coefs >= 0.) and np.all(coefs <= 1.) + + # Get info from config + Types = [] + # Ds = + # us = + # indouts = + # kouts = + # vperp = + # coefs = + + Ds[0,:,:] = self.D + self._dgeom['kOut'][None,:] * self.u + us[0,:,:] = self.config._reflect_geom(Ds[ii-1,:,:]) + kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. + coefs[ii,:] = self.config._get_reflections_coefs() + + for ii in range(1,nb): + Ds[ii,:,:], us[ii,:,:] = self.config._reflect_geom(Ds[ii-1,:,:]) + kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. + coefs[ii,:] = self.config._reflect_coefs(coefs=coefs, + indout=indouts[ii,:,:]) + + self._dgeom['dreflect'] = None + + + ########### # strip dictionaries ########### diff --git a/tofu/version.py b/tofu/version.py index 26e3e8425..cb34e60f5 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.0-592-gb47bb5a' \ No newline at end of file +__version__='1.4.0-594-gea5ba4f' \ No newline at end of file From 27e2383faf22be87464e4ee9f429e88afaf08f66 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Mon, 26 Aug 2019 14:22:23 +0200 Subject: [PATCH 02/19] [Issue 158] Config._reflect_Types() and Config._reflect_geom() both working. Now bug in Rays._complete_dX12() for CamLOS1D from reflections --- tofu/geom/_core.py | 36 +++++++++++++++++++++--------------- tofu/version.py | 2 +- 2 files changed, 22 insertions(+), 16 deletions(-) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 72d8a919b..d8361e9bd 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -49,7 +49,7 @@ _PHITHETAPROJ_NTHETA = 1000 _RES = 0.005 _NTHREADS = 16 - +_DREFLECT = {'specular':0, 'diffusive':1, 'ccube':2} """ ############################################################################### @@ -2546,13 +2546,13 @@ def fdistfromwall(self, r, z, phi): # Method handling reflections - def _reflect_Types(self, indout=None, Type=None): + def _reflect_Types(self, indout=None, Type=None, nRays=None): """ Return an array indicating the Type of reflection for each LOS Return a (nRays,) np.ndarray of int indices, each index corresponds to: - 0: specular reflections - 1: diffusive reflections - - 2: corner cube reflections + - 2: ccube reflections (corner cube) If indout is provided, the Types are computed according to the information stored in each corresponding Struct @@ -2560,35 +2560,43 @@ def _reflect_Types(self, indout=None, Type=None): If Type is provided, the Type is forced (user-defined) for all LOS """ - nRays = u.shape[1] if Type is not None: - assert Type in ['specular', 'diffusive', 'corner cube'] - Types = np.full((nRays,), _DREFLECT[Type],, dtype=int) + assert Type in ['specular', 'diffusive', 'ccube'] + Types = np.full((nRays,), _DREFLECT[Type], dtype=int) else: Types = None return Types def _reflect_geom(self, u, vperp, indout=None, Type=None): + assert u.shape == vperp.shape and u.shape[0] == 3 + if indout is not None: + assert indout.shape == (3,u.shape[1]) # Get Types of relection for each Ray - Types = self._reflect_Types(indout=indout, Type=Type) + Types = self._reflect_Types(indout=indout, Type=Type, nRays=u.shape[1]) # Get indices of each Type indspec = Types == 0 inddiff = Types == 1 indcorn = Types == 2 + u2 = np.full(u.shape, np.nan) if np.any(np.logical_or(indspec,inddiff)): - sca = np.sum(u*vperp, axis=0) - vpar = np.array([u[1,:]*vperp[2,:] - u[2,:]*vpepr[1,:], - u[2,:]*vperp[0,:] - u[0,:]*vperp[2,:], - u[0,:]*vperp[1,:] - u[1,:]*vperp[0,:]]) + vpar = np.array([vperp[1,:]*u[2,:] - vperp[2,:]*u[1,:], + vperp[2,:]*u[0,:] - vperp[0,:]*u[2,:], + vperp[0,:]*u[1,:] - vperp[1,:]*u[0,:]]) + vpar = np.array([vpar[1,:]*vperp[2,:] - vpar[2,:]*vperp[1,:], + vpar[2,:]*vperp[0,:] - vpar[0,:]*vperp[2,:], + vpar[0,:]*vperp[1,:] - vpar[1,:]*vperp[0,:]]) + vpar = vpar / np.sqrt(np.sum(vpar**2, axis=0))[None,:] if np.any(indspec): # Compute u2 for specular sca = np.sum(u[:,indspec]*vperp[:,indspec],axis=0,keepdims=True) sca2 = np.sum(u[:,indspec]*vpar[:,indspec],axis=0,keepdims=True) + assert np.all(sca<=0.) and np.all(sca>=-1.) + assert np.all(sca2>=0.) and np.all(sca<=1.) u2[:,indspec] = - sca*vperp[:,indspec] + sca2*vpar[:,indspec] if np.any(inddiff): @@ -3805,8 +3813,6 @@ def add_reflections(self, Type=None, nb=None, coefs=None): """ # Check inputs - if Type is not None: - assert Type in ['speculiar', 'diffusive', 'corner cube'] if nb is None: nb = 1 nb = int(nb) @@ -3829,12 +3835,12 @@ def add_reflections(self, Type=None, nb=None, coefs=None): Ds[0,:,:] = self.D + self._dgeom['kOut'][None,:] * self.u us[0,:,:] = self.config._reflect_geom(Ds[ii-1,:,:]) - kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. + # kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. coefs[ii,:] = self.config._get_reflections_coefs() for ii in range(1,nb): Ds[ii,:,:], us[ii,:,:] = self.config._reflect_geom(Ds[ii-1,:,:]) - kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. + # kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. coefs[ii,:] = self.config._reflect_coefs(coefs=coefs, indout=indouts[ii,:,:]) diff --git a/tofu/version.py b/tofu/version.py index cb34e60f5..e00d567b4 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.0-594-gea5ba4f' \ No newline at end of file +__version__='1.4.0-595-g0e57448' \ No newline at end of file From 8b9012bdfdfbb8a805451560d63833894e715d65 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Mon, 26 Aug 2019 15:27:34 +0200 Subject: [PATCH 03/19] [Issue 158] Bug in Rays._complete_dX12() fixed for not aligned CamLOS1D --- tofu/geom/_core.py | 3 ++- tofu/version.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index d8361e9bd..b3631f098 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -3389,8 +3389,8 @@ def _complete_dX12(self, dgeom): if dgeom['dX12'] is None: dgeom['dX12'] = {} dgeom['dX12'].update({'e1':e1, 'x1':x1, 'n1':x1.size}) - else: + elif self._is2D(): ind = np.nanargmax(vect2) v1 = van[:,ind] nn = np.cross(v0,v1) @@ -3401,6 +3401,7 @@ def _complete_dX12(self, dgeom): # check nIn orientation sca = np.sum(self.u*nn[:,np.newaxis],axis=0) lc = [np.all(sca>=0.), np.all(sca<=0.)] + assert any(lc) nIn = nn if lc[0] else -nn e1 = v0 diff --git a/tofu/version.py b/tofu/version.py index e00d567b4..acf0b567c 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.0-595-g0e57448' \ No newline at end of file +__version__='1.4.0-596-g27e2383' \ No newline at end of file From ac65ced90716f008ace7062626a9058ccc9c0233 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Mon, 26 Aug 2019 15:41:36 +0200 Subject: [PATCH 04/19] [Issue 158] Rays.get_reflections_as_cam() operational (returns CamLOS1D) --- tofu/geom/_core.py | 12 ++++++++++++ tofu/version.py | 2 +- 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index b3631f098..59c62fb74 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -3803,6 +3803,18 @@ def _set_dmisc(self, color=None): # Reflections ########### + def get_reflections_as_cam(self, Type=None, Name=None): + Ds = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u + us, Types = self.config._reflect_geom(self.u, self._dgeom['vperp'], + Type=Type) + clas = Rays if self.__class__.__name__ == Rays else CamLOS1D + if Name is None: + Name = self.Id.Name + '_Reflect%s'%str(Type) + return clas(dgeom=(Ds,us), config=self.config, + Exp=self.Id.Exp, Diag=self.Id.Diag, + Name=Name, shot=self.Id.shot) + + def add_reflections(self, Type=None, nb=None, coefs=None): """ Add relfected LOS to the camera diff --git a/tofu/version.py b/tofu/version.py index acf0b567c..568fd226e 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.0-596-g27e2383' \ No newline at end of file +__version__='1.4.0-597-g8b9012b' \ No newline at end of file From 4eb2ce5471c13b9838bccbe30024d2311d69b3de Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Mon, 26 Aug 2019 15:52:31 +0200 Subject: [PATCH 05/19] [Issue 158] Corrected bug in Config._reflect_geom(Type='diffusive') (used to return only posivitive parallel components) --- tofu/geom/_core.py | 6 +++--- tofu/version.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 59c62fb74..1ae897b5e 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -2601,9 +2601,9 @@ def _reflect_geom(self, u, vperp, indout=None, Type=None): if np.any(inddiff): # Compute u2 for diffusive - sca = np.random.random((1,inddiff.sum())) - u2[:,inddiff] = (sca * vperp[:,inddiff] - + np.sqrt(1.-sca**2) * vpar[:,inddiff]) + sca = 2.*(np.random.random((1,inddiff.sum()))-0.5) + u2[:,inddiff] = (np.sqrt(1.-sca**2) * vperp[:,inddiff] + + sca * vpar[:,inddiff]) if np.any(indcorn): u2[:,indcorn] = -u[:,indcorn] diff --git a/tofu/version.py b/tofu/version.py index 568fd226e..6e339ed9c 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.0-597-g8b9012b' \ No newline at end of file +__version__='1.4.0-598-gac65ced' \ No newline at end of file From c1bd190dd67f542ee36c8e76ed1630c8ed4d0bd2 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Mon, 26 Aug 2019 16:15:04 +0200 Subject: [PATCH 06/19] [Issue 158] Rays.add_reflections() needs to be finished / tested / debugged --- tofu/geom/_core.py | 59 +++++++++++++++++++++++++++++----------------- tofu/version.py | 2 +- 2 files changed, 39 insertions(+), 22 deletions(-) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 1ae897b5e..63c04ffbb 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -3804,12 +3804,23 @@ def _set_dmisc(self, color=None): ########### def get_reflections_as_cam(self, Type=None, Name=None): + """ Return a camera made of reflected LOS + + Reflected LOS can be of 3 types: + - 'speculiar': standard mirror-like reflection + - 'diffusive': random reflection + - 'ccube': corner-cube reflection (ray goes back its way) + + As opposed to self.add_reflections(), the reflected rays are + return as an independent camera (CamLOS1D) + + """ Ds = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u us, Types = self.config._reflect_geom(self.u, self._dgeom['vperp'], Type=Type) - clas = Rays if self.__class__.__name__ == Rays else CamLOS1D if Name is None: Name = self.Id.Name + '_Reflect%s'%str(Type) + clas = Rays if self.__class__.__name__ == Rays else CamLOS1D return clas(dgeom=(Ds,us), config=self.config, Exp=self.Id.Exp, Diag=self.Id.Diag, Name=Name, shot=self.Id.shot) @@ -3819,9 +3830,12 @@ def add_reflections(self, Type=None, nb=None, coefs=None): """ Add relfected LOS to the camera Reflected LOS can be of 3 types: - - 'speculiar' - - 'diffusive' - - 'corer cube' + - 'speculiar': standard mirror-like reflection + - 'diffusive': random reflection + - 'ccube': corner-cube reflection (ray goes back its way) + + As opposed to self.get_reflections_as_cam(), the reflected rays are + stored in the camera object """ @@ -3838,26 +3852,29 @@ def add_reflections(self, Type=None, nb=None, coefs=None): assert np.all(coefs >= 0.) and np.all(coefs <= 1.) # Get info from config - Types = [] - # Ds = - # us = - # indouts = - # kouts = - # vperp = - # coefs = - - Ds[0,:,:] = self.D + self._dgeom['kOut'][None,:] * self.u - us[0,:,:] = self.config._reflect_geom(Ds[ii-1,:,:]) - # kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. - coefs[ii,:] = self.config._get_reflections_coefs() + nRays = self.Rays + Types = np.full((nb,nRays), 0, dtype=int) + us = np.full((nb,3,nRays), np.nan, dtype=float) + kouts = np.full((nb,nRays), np.nan, dtype=float) + indouts = np.full((nb,3,nRays), 0, dtype=int) + vperp = np.full((nb,3,nRays), np.nan, dtype=float) + + # Prepare inputs + largs, dkwd = self._prepare_inputs_kInOut() + + Ds[0,:,:] = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u + us[0,:,:], Types[0,:] = self.config._reflect_geom(self.u, + self._dgeom['vperp'], + Type=Type) + outi = _GG.LOS_Calc_PInOut_VesStruct(*largs, **dkwd)[1:] + kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = outi for ii in range(1,nb): - Ds[ii,:,:], us[ii,:,:] = self.config._reflect_geom(Ds[ii-1,:,:]) - # kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. - coefs[ii,:] = self.config._reflect_coefs(coefs=coefs, - indout=indouts[ii,:,:]) + Ds = self.config._reflect_geom(Ds[ii-1,:,:]) + kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. - self._dgeom['dreflect'] = None + self._dgeom['dreflect'] = {'nb':nb, 'Type':Type, 'Types':Types, + 'us':us, 'kouts':kouts, 'indouts':indouts} diff --git a/tofu/version.py b/tofu/version.py index 6e339ed9c..a29ee71f2 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.0-598-gac65ced' \ No newline at end of file +__version__='1.4.0-599-g4eb2ce5' \ No newline at end of file From 99a6895cb68e36329f5341676adc73e255fc390e Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Mon, 26 Aug 2019 22:06:26 +0200 Subject: [PATCH 07/19] [Issue 158] Rays.get_reflections_as_cam() and Rays.add_reflections() operational, TODO: plots --- tofu/geom/_core.py | 126 ++++++++++++++++++++++++++------------------- tofu/version.py | 2 +- 2 files changed, 74 insertions(+), 54 deletions(-) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 63c04ffbb..58e532e50 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -3190,7 +3190,7 @@ def _get_keys_dgeom(): lk = ['D','u','pinhole', 'nRays', 'kIn', 'kOut', 'PkIn', 'PkOut', 'vperp', 'indout', 'kRMin', 'PRMin', 'RMin', 'isImage', - 'Etendues', 'Surfaces', 'dX12'] + 'Etendues', 'Surfaces', 'dX12', 'dreflect'] return lk @staticmethod @@ -3463,26 +3463,35 @@ def _complete_dX12(self, dgeom): - def _prepare_inputs_kInOut(self): - if self._method=='ref': - # Prepare input + def _prepare_inputs_kInOut(self, D=None, u=None): + + # Prepare input + if D is None: D = np.ascontiguousarray(self.D) + else: + D = np.ascontiguousarray(D) + if u is None: u = np.ascontiguousarray(self.u) + else: + u = np.ascontiguousarray(u) + assert D.shape == u.shape - # Get reference - lS = self.lStruct_computeInOut + # Get reference + lS = self.lStruct_computeInOut - lSIn = [ss for ss in lS if ss._InOut=='in'] - if len(lSIn)==0: - msg = "self.config must have at least a StructIn subclass !" - assert len(lSIn)>0, msg - elif len(lSIn)>1: - S = lSIn[np.argmin([ss.dgeom['Surf'] for ss in lSIn])] - else: - S = lSIn[0] + lSIn = [ss for ss in lS if ss._InOut=='in'] + if len(lSIn)==0: + msg = "self.config must have at least a StructIn subclass !" + assert len(lSIn) > 0, msg + iref = np.argmin([ss.dgeom['Surf'] for ss in lSIn]) + S = lSIn[iref] + + VPoly = S.Poly_closed + VVIn = S.dgeom['VIn'] + largs = [D, u, VPoly, VVIn] + + if self._method=='ref': - VPoly = S.Poly_closed - VVIn = S.dgeom['VIn'] Lim = S.Lim nLim = S.noccur VType = self.config.Id.Type @@ -3494,7 +3503,6 @@ def _prepare_inputs_kInOut(self): lSVIn.append(ss.dgeom['VIn']) lSLim.append(ss.Lim) lSnLim.append(ss.noccur) - largs = [D, u, VPoly, VVIn] dkwd = dict(Lim=Lim, nLim=nLim, LSPoly=lSPoly, LSLim=lSLim, lSnLim=lSnLim, LSVIn=lSVIn, VType=VType, @@ -3502,23 +3510,7 @@ def _prepare_inputs_kInOut(self): EpsA=1.e-9, EpsB=1.e-9, EpsPlane=1.e-9, Test=True) elif self._method=='optimized': - # Prepare input - D = np.ascontiguousarray(self.D) - u = np.ascontiguousarray(self.u) - # Get reference - lS = self.lStruct_computeInOut - - lSIn = [ss for ss in lS if ss._InOut=='in'] - if len(lSIn)==0: - msg = "self.config must have at least a StructIn subclass !" - assert len(lSIn)>0, msg - elif len(lSIn)>1: - S = lSIn[np.argmin([ss.dgeom['Surf'] for ss in lSIn])] - else: - S = lSIn[0] - VPoly = S.Poly_closed - VVIn = S.dgeom['VIn'] if np.size(np.shape(S.Lim)) > 1 : Lim = np.asarray([S.Lim[0][0], S.Lim[0][1]]) else: @@ -3560,7 +3552,6 @@ def _prepare_inputs_kInOut(self): lSVInx = np.asarray(lSVInx) lSVIny = np.asarray(lSVIny) - largs = [D, u, VPoly, VVIn] dkwd = dict(ves_lims=Lim, nstruct_tot=num_tot_structs, nstruct_lim=num_lim_structs, @@ -3575,12 +3566,6 @@ def _prepare_inputs_kInOut(self): rmin=-1, forbid=True, eps_uz=1.e-6, eps_vz=1.e-9, eps_a=1.e-9, eps_b=1.e-9, eps_plane=1.e-9, test=True) - else: - # -------------------------------- - # Here I can prepare the inputs as requested by your routine - pass - # -------------------------------- - return largs, dkwd def _compute_kInOut(self): @@ -3803,7 +3788,7 @@ def _set_dmisc(self, color=None): # Reflections ########### - def get_reflections_as_cam(self, Type=None, Name=None): + def get_reflections_as_cam(self, Type=None, Name=None, nb=None): """ Return a camera made of reflected LOS Reflected LOS can be of 3 types: @@ -3815,15 +3800,45 @@ def get_reflections_as_cam(self, Type=None, Name=None): return as an independent camera (CamLOS1D) """ + # Check inputs + if nb is None: + nb = 1 + nb = int(nb) + assert nb > 0 Ds = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u us, Types = self.config._reflect_geom(self.u, self._dgeom['vperp'], Type=Type) if Name is None: Name = self.Id.Name + '_Reflect%s'%str(Type) clas = Rays if self.__class__.__name__ == Rays else CamLOS1D - return clas(dgeom=(Ds,us), config=self.config, + + # Run first iteration + Types = np.full((nb,self.nRays), 0, dtype=int) + Ds = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u + us, Types[0,:] = self.config._reflect_geom(self.u, + self._dgeom['vperp'], + Type=Type) + lcam = [clas(dgeom=(Ds,us), config=self.config, Exp=self.Id.Exp, Diag=self.Id.Diag, - Name=Name, shot=self.Id.shot) + Name=Name, shot=self.Id.shot)] + if nb == 1: + return lcam[0], Types[0,:] + + largs, dkwd = self._prepare_inputs_kInOut(D=Ds, u=us) + kouts, vperps, indouts = _GG.LOS_Calc_PInOut_VesStruct(*largs, **dkwd)[1:] + + # Run other iterations + for ii in range(1,nb): + Ds = Ds + (kouts[None,:]-1.e-12) * us + us, Types[ii,:] = self.config._reflect_geom(us, vperps, + Type=Type) + kouts, vperps, indouts = _GG.LOS_Calc_PInOut_VesStruct(Ds, us, + *largs[2:], + **dkwd)[1:] + lcam.append(clas(dgeom=(Ds,us), config=self.config, + Exp=self.Id.Exp, Diag=self.Id.Diag, + Name=Name, shot=self.Id.shot)) + return lcam, Types def add_reflections(self, Type=None, nb=None, coefs=None): @@ -3851,27 +3866,32 @@ def add_reflections(self, Type=None, nb=None, coefs=None): assert coefs.shape == (self.nRays,) assert np.all(coefs >= 0.) and np.all(coefs <= 1.) - # Get info from config - nRays = self.Rays + # Prepare output + nRays = self.nRays Types = np.full((nb,nRays), 0, dtype=int) us = np.full((nb,3,nRays), np.nan, dtype=float) kouts = np.full((nb,nRays), np.nan, dtype=float) indouts = np.full((nb,3,nRays), 0, dtype=int) - vperp = np.full((nb,3,nRays), np.nan, dtype=float) + vperps = np.full((nb,3,nRays), np.nan, dtype=float) - # Prepare inputs - largs, dkwd = self._prepare_inputs_kInOut() - - Ds[0,:,:] = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u + # Run first iteration + Ds = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u us[0,:,:], Types[0,:] = self.config._reflect_geom(self.u, self._dgeom['vperp'], Type=Type) + largs, dkwd = self._prepare_inputs_kInOut(D=Ds, u=us[0,:,:]) outi = _GG.LOS_Calc_PInOut_VesStruct(*largs, **dkwd)[1:] - kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = outi + kouts[0,:], vperps[0,:,:], indouts[0,:,:] = outi + # Run other iterations for ii in range(1,nb): - Ds = self.config._reflect_geom(Ds[ii-1,:,:]) - kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = self. + Ds = Ds + (kouts[ii-1:ii,:]-1.e-12) * us[ii-1,:,:] + us[ii,:,:], Types[ii,:] = self.config._reflect_geom(us[ii-1,:,:], + vperps[ii-1,:,:], + Type=Type) + outi = _GG.LOS_Calc_PInOut_VesStruct(Ds, us[ii,:,:], + *largs[2:], **dkwd)[1:] + kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = outi self._dgeom['dreflect'] = {'nb':nb, 'Type':Type, 'Types':Types, 'us':us, 'kouts':kouts, 'indouts':indouts} diff --git a/tofu/version.py b/tofu/version.py index a29ee71f2..1e8c74911 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.0-599-g4eb2ce5' \ No newline at end of file +__version__='1.4.0-600-gc1bd190' \ No newline at end of file From cc321d1e097cbca5f05ce77c2de3c9bdde09b8bd Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Tue, 27 Aug 2019 09:25:27 +0200 Subject: [PATCH 08/19] [Issue 158] Started updating plots => Rays.plot() => tf.geom._plot.Rays_plot() => Rays._get_plotL() => tf.geom._comp.LOS_CrossProj(), TBF --- tofu/geom/_comp.py | 63 +++++++++++++++++++++++++++++++--------------- tofu/geom/_core.py | 34 +++++++++++++++++-------- tofu/geom/_plot.py | 19 +++++++++----- tofu/version.py | 2 +- 4 files changed, 81 insertions(+), 37 deletions(-) diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index a5dc73ccd..dfc675ab3 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -482,7 +482,7 @@ def LOS_PRMin(Ds, dus, kPOut=None, Eps=1.e-12, Test=True): return kRMin -def LOS_CrossProj(VType, Ds, us, kPIns, kPOuts, kRMins, +def LOS_CrossProj(VType, Ds, us, kIns, kOuts, kRMins, Lplot='In', proj='All', multi=False): """ Compute the parameters to plot the poloidal projection of the LOS """ assert type(VType) is str and VType.lower() in ['tor','lin'] @@ -490,34 +490,57 @@ def LOS_CrossProj(VType, Ds, us, kPIns, kPOuts, kRMins, assert type(proj) is str proj = proj.lower() assert proj in ['cross','hor','all','3d'] - assert Ds.ndim==2 and Ds.shape==us.shape - nL = Ds.shape[1] - k0 = kPIns if Lplot.lower()=='in' else np.zeros((nL,)) + lc = [Ds.ndim in [2,3], Ds.shape == us.shape] + if not all(lc): + msg = "Ds and us must have the same shape and dim in [2,3]:\n" + msg += " - provided Ds.shape: %s\n"%str(Ds.shape) + msg += " - provided us.shape: %s"%str(us.shape) + raise Exception(msg) + lc = [kIns.shape == kOuts.shape, + kIns.shape[-1] == Ds.shape[-1], kIns.ndim == Ds.ndim-1] + if not all(lc): + msg = "kIns and kOuts must have the same shape and ndim = Ds.ndim-1:\n" + msg += " - Ds.shape : %s\n"%str(Ds.shape) + msg += " - kIns.shape : %s\n"%str(kIns.shape) + msg += " - kOutss.shape: %s"%str(kOuts.shape) + raise Exception(msg) + + # Prepare inputs + if Ds.ndim == 2: + Ds = np.array([Ds]) + us = np.array([us]) + kIns = np.array([kIns]) + kOuts = np.array([kOuts]) + + nref, _, nL = Ds.shape + + k0 = kIns if Lplot.lower() == 'in' else np.zeros((nL,)) if VType.lower()=='tor' and proj in ['cross','all']: - CrossProjAng = np.arccos(np.sqrt(us[0,:]**2+us[1,:]**2) + CrossProjAng = np.arccos(np.sqrt(us[:,0,:]**2+us[:,1,:]**2) /np.sqrt(np.sum(us**2,axis=0))) nkp = np.ceil(25.*(1 - (CrossProjAng/(np.pi/4)-1)**2) + 2) - ks = np.max([kRMins,kPIns],axis=0) if Lplot.lower()=='in' else kRMins + ks = np.max([kRMins,kIns],axis=0) if Lplot.lower()=='in' else kRMins pts0 = [] if multi: for ii in range(0,nL): - if np.isnan(kPOuts[ii]): - pts0.append( np.array([[np.nan,np.nan], - [np.nan,np.nan]]) ) - else: - k = np.linspace(k0[ii],kPOuts[ii],nkp[ii],endpoint=True) - k = np.unique(np.append(k,ks[ii])) - pp = Ds[:,ii:ii+1] + k[np.newaxis,:]*us[:,ii:ii+1] - pts0.append( np.array([np.hypot(pp[0,:],pp[1,:]),pp[2,:]]) ) + for jj in range(0,nref): + if np.isnan(kOuts[jj,ii]): + pts0.append( np.array([[np.nan,np.nan], + [np.nan,np.nan]]) ) + else: + k = np.linspace(k0[ii], kOuts[jj,ii], nkp[ii], endpoint=True) + k = np.unique(np.append(k,ks[ii])) + pp = Ds[:,ii:ii+1] + k[np.newaxis,:]*us[:,ii:ii+1] + pts0.append( np.array([np.hypot(pp[0,:],pp[1,:]),pp[2,:]]) ) else: for ii in range(0,nL): - if np.isnan(kPOuts[ii]): + if np.isnan(kOuts[ii]): pts0.append(np.array([[np.nan,np.nan,np.nan], [np.nan,np.nan,np.nan], [np.nan,np.nan,np.nan]])) else: - k = np.linspace(k0[ii],kPOuts[ii],nkp[ii],endpoint=True) + k = np.linspace(k0[ii],kOuts[ii],nkp[ii],endpoint=True) k = np.append(np.unique(np.append(k,ks[ii])),np.nan) pts0.append( Ds[:,ii:ii+1] + k[np.newaxis,:]*us[:,ii:ii+1] ) pts0 = np.concatenate(tuple(pts0),axis=1) @@ -527,21 +550,21 @@ def LOS_CrossProj(VType, Ds, us, kPIns, kPOuts, kRMins, pts = [] if multi: for ii in range(0,nL): - if np.isnan(kPOuts[ii]): + if np.isnan(kOuts[ii]): pts.append( np.array([[np.nan,np.nan], [np.nan,np.nan], [np.nan,np.nan]]) ) else: - k = np.array([k0[ii],kPOuts[ii]]) + k = np.array([k0[ii],kOuts[ii]]) pts.append( Ds[:,ii:ii+1] + k[np.newaxis,:]*us[:,ii:ii+1] ) else: for ii in range(0,nL): - if np.isnan(kPOuts[ii]): + if np.isnan(kOuts[ii]): pts.append(np.array([[np.nan,np.nan,np.nan], [np.nan,np.nan,np.nan], [np.nan,np.nan,np.nan]])) else: - k = np.array([k0[ii],kPOuts[ii],np.nan]) + k = np.array([k0[ii],kOuts[ii],np.nan]) pts.append( Ds[:,ii:ii+1] + k[np.newaxis,:]*us[:,ii:ii+1] ) pts = np.concatenate(tuple(pts),axis=1) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 58e532e50..24ab5b729 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -3920,12 +3920,12 @@ def _strip_dgeom(self, strip=0): if strip==1: lkeep = ['D','u','pinhole','nRays', 'kIn','kOut','vperp','indout', 'kRMin', - 'Etendues','Surfaces','isImage','dX12'] + 'Etendues','Surfaces','isImage','dX12', 'dreflect'] utils.ToFuObject._strip_dict(self._dgeom, lkeep=lkeep) elif self._dstrip['strip']<=1 and strip>=2: lkeep = ['D','u','pinhole','nRays', 'kIn','kOut','vperp','indout', - 'Etendues','Surfaces','isImage','dX12'] + 'Etendues','Surfaces','isImage','dX12', 'dreflect'] utils.ToFuObject._strip_dict(self._dgeom, lkeep=lkeep) def _strip_dconfig(self, strip=0, verb=True): @@ -4358,7 +4358,7 @@ def get_subset(self, indch=None, Name=None): obj = self.__class__(fromdict=dd) return obj - def _get_plotL(self, Lplot='Tot', proj='All', ind=None, multi=False): + def _get_plotL(self, reflections=True, Lplot='Tot', proj='All', ind=None, multi=False): """ Get the (R,Z) coordinates of the cross-section projections """ ind = self._check_indch(ind) if ind.size>0: @@ -4370,9 +4370,18 @@ def _get_plotL(self, Lplot='Tot', proj='All', ind=None, multi=False): kRMin = self._dgeom['kRMin'][ind] else: kRMin = None - pts = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, - kPIn, kPOut, kRMin, proj=proj, - Lplot=Lplot, multi=multi) + + # Add reflections ? + c0 = (reflections and self._dgeom.get('dreflect') is not None + and self._dgeom['dreflect'].get('us') is not None) + if c0: + pts = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, + kPIn, kPOut, kRMin, proj=proj, + Lplot=Lplot, multi=multi) + else: + pts = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, + kPIn, kPOut, kRMin, proj=proj, + Lplot=Lplot, multi=multi) else: pts = None return pts @@ -4923,7 +4932,8 @@ def calc_signal_from_Plasma2D(self, plasma2d, t=None, newcalc=False, - def plot(self, lax=None, proj='all', Lplot=_def.LOSLplot, element='L', + def plot(self, lax=None, proj='all', reflections=True, + Lplot=_def.LOSLplot, element='L', element_config='P', Leg='', dL=None, dPtD=_def.LOSMd, dPtI=_def.LOSMd, dPtO=_def.LOSMd, dPtR=_def.LOSMd, dPtP=_def.LOSMd, dLeg=_def.TorLegd, multi=False, ind=None, @@ -4937,15 +4947,18 @@ def plot(self, lax=None, proj='all', Lplot=_def.LOSLplot, element='L', Parameters ---------- - Lax : list / plt.Axes + lax : list / plt.Axes The axes for plotting (list of 2 axes if Proj='All') If None a new figure with new axes is created - Proj : str + proj : str Flag specifying the kind of projection: - 'Cross' : cross-section - 'Hor' : horizontal - 'All' : both cross-section and horizontal (on 2 axes) - '3d' : a (matplotlib) 3d plot + projections:bool + Flag indicating whether to plot also the reflected rays + Assuming some reflected rays are present (self.add_reflections()) element : str Flag specifying which elements to plot Each capital letter corresponds to an element: @@ -5004,7 +5017,8 @@ def plot(self, lax=None, proj='all', Lplot=_def.LOSLplot, element='L', """ - return _plot.Rays_plot(self, Lax=lax, Proj=proj, Lplot=Lplot, + return _plot.Rays_plot(self, Lax=lax, Proj=proj, + reflections=reflections, Lplot=Lplot, element=element, element_config=element_config, Leg=Leg, dL=dL, dPtD=dPtD, dPtI=dPtI, dPtO=dPtO, dPtR=dPtR, dPtP=dPtP, dLeg=dLeg, multi=multi, ind=ind, diff --git a/tofu/geom/_plot.py b/tofu/geom/_plot.py index e0452874e..3657c6395 100644 --- a/tofu/geom/_plot.py +++ b/tofu/geom/_plot.py @@ -862,8 +862,8 @@ def Get_FieldsFrom_LLOS(L,Fields): -def Rays_plot(GLos, Lax=None, Proj='all', Lplot=_def.LOSLplot, - element='LDIORP', element_config='P', +def Rays_plot(GLos, Lax=None, Proj='all', reflections=True, + Lplot=_def.LOSLplot, element='LDIORP', element_config='P', Leg=None, dL=None, dPtD=_def.LOSMd, dPtI=_def.LOSMd, dPtO=_def.LOSMd, dPtR=_def.LOSMd, dPtP=_def.LOSMd, dLeg=_def.TorLegd, multi=False, @@ -918,7 +918,9 @@ def Rays_plot(GLos, Lax=None, Proj='all', Lplot=_def.LOSLplot, if len(ind)>0 and not element=='': if Proj=='3d': - Lax[0] = _Rays_plot_3D(GLos, ax=Lax[0], Elt=element, Lplot=Lplot, + Lax[0] = _Rays_plot_3D(GLos, ax=Lax[0], + reflections=reflections, + Elt=element, Lplot=Lplot, Leg=Leg, dL=dL, dPtD=dPtD, dPtI=dPtI, dPtO=dPtO, dPtR=dPtR, dPtP=dPtP, dLeg=None, multi=multi, ind=ind, @@ -929,7 +931,9 @@ def Rays_plot(GLos, Lax=None, Proj='all', Lplot=_def.LOSLplot, fs=fs, wintit=wintit, Type=GLos.config.Type)) if Proj in ['cross','all']: - Lax[0] = _Rays_plot_Cross(GLos, ax=Lax[0], Elt=element, Lplot=Lplot, + Lax[0] = _Rays_plot_Cross(GLos, ax=Lax[0], + reflections=reflections, + Elt=element, Lplot=Lplot, Leg=Leg, dL=dL, dPtD=dPtD, dPtI=dPtI, dPtO=dPtO, dPtR=dPtR, dPtP=dPtP, dLeg=None, multi=multi, ind=ind, @@ -937,7 +941,9 @@ def Rays_plot(GLos, Lax=None, Proj='all', Lplot=_def.LOSLplot, Test=Test) if Proj in ['hor','all']: ii = 0 if Proj=='hor' else 1 - Lax[ii] = _Rays_plot_Hor(GLos, ax=Lax[ii], Elt=element, Lplot=Lplot, + Lax[ii] = _Rays_plot_Hor(GLos, ax=Lax[ii], + reflections=reflections, + Elt=element, Lplot=Lplot, Leg=Leg, dL=dL, dPtD=dPtD, dPtI=dPtI, dPtO=dPtO, dPtR=dPtR, dPtP=dPtP, dLeg=None, multi=multi, ind=ind, @@ -952,7 +958,8 @@ def Rays_plot(GLos, Lax=None, Proj='all', Lplot=_def.LOSLplot, -def _Rays_plot_Cross(L,Leg=None,Lplot='Tot', Elt='LDIORP',ax=None, +def _Rays_plot_Cross(L,Leg=None, reflections=True, + Lplot='Tot', Elt='LDIORP',ax=None, dL=_def.LOSLd, dPtD=_def.LOSMd, dPtI=_def.LOSMd, dPtO=_def.LOSMd, dPtR=_def.LOSMd, dPtP=_def.LOSMd, dLeg=_def.TorLegd, multi=False, ind=None, diff --git a/tofu/version.py b/tofu/version.py index 1e8c74911..a988e6bb6 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.0-600-gc1bd190' \ No newline at end of file +__version__='1.4.0-601-g99a6895' \ No newline at end of file From c17968a97c3bc734782b176522b0d01adf20222d Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Tue, 27 Aug 2019 10:57:17 +0200 Subject: [PATCH 09/19] [Issue 158] Revised Rays.add_reflections() to put reflections as last dimension, and updated Rays._get_plotL() and tf.geom._comp.LOS_PRMin() accordingly => test and update tf.geom._plot.Rays_plot() and Rays.plot() --- tofu/geom/_comp.py | 63 +++++++++++++++++++--------------- tofu/geom/_core.py | 85 +++++++++++++++++++++++++--------------------- tofu/version.py | 2 +- 3 files changed, 84 insertions(+), 66 deletions(-) diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index dfc675ab3..315d1dab2 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -449,41 +449,51 @@ def _get_phithetaproj_dist(poly_closed, ax, Dtheta, nDtheta, ############################################################################### """ -def LOS_PRMin(Ds, dus, kPOut=None, Eps=1.e-12, Test=True): +def LOS_PRMin(Ds, us, kPOut=None, Eps=1.e-12, squeeze=True, Test=True): """ Compute the point on the LOS where the major radius is minimum """ if Test: - assert Ds.ndim in [1,2] and 3 in Ds.shape and Ds.shape==dus.shape - assert kPOut is None or (Ds.ndim==1 and not hasattr(kPOut,'__iter__')) or (Ds.ndim==2 and kPOut.shape==(Ds.size/3,)) - - v = Ds.ndim==1 - if v: - Ds = Ds.reshape((3,1)) - dus = dus.reshape((3,1)) - if kPOut is not None: - kPOut = np.array([kPOut]) - - kRMin = np.nan*np.ones((Ds.shape[1],)) - uparN = np.sqrt(dus[0,:]**2 + dus[1,:]**2) + assert Ds.ndim in [1,2,3] and 3 in Ds.shape and Ds.shape == us.shape + if kOut is not None: + kOut = np.atleast_1d(kOut) + assert kOut.size == Ds.size/3 + + v = Ds.ndim == 1 + if Ds.ndim == 1: + Ds, us = Ds[:,None,None], us[:,None,None] + elif Ds.ndim == 2: + Ds, us = Ds[:,:,None], us[:,:,None] + if kOut is not None: + kOut = kOut[:,None] + nref, _, nlos = Ds.shape + + kRMin = np.full((nlos,nref), np.nan) + uparN = np.sqrt(us[0,:,:]**2 + us[1,:,:]**2) # Case with u vertical - ind = uparN>Eps + ind = uparN > Eps kRMin[~ind] = 0. # Else - kRMin[ind] = -(dus[0,ind]*Ds[0,ind]+dus[1,ind]*Ds[1,ind])/uparN[ind]**2 + kRMin[ind] = -(us[0,ind]*Ds[0,ind] + us[1,ind]*Ds[1,ind]) / uparN[ind]**2 # Check - kRMin[kRMin<=0.] = 0. - if kPOut is not None: - kRMin[kRMin>kPOut] = kPOut[kRMin>kPOut] - - if v: - kRMin = kRMin[0] + kRMin[kRMin <= 0.] = 0. + if kOut is not None: + kRMin[kRMin > kOut] = kOut[kRMin > kOut] + + # squeeze + if squeeze: + if nref == 1 and nlos == 11: + kRMin = kRMin[0,0] + elif nref == 1: + kRMin = kRMin[:,0] + elif nlos == 1: + kRMin = kRMin[0,:] return kRMin def LOS_CrossProj(VType, Ds, us, kIns, kOuts, kRMins, - Lplot='In', proj='All', multi=False): + proj='All', multi=False): """ Compute the parameters to plot the poloidal projection of the LOS """ assert type(VType) is str and VType.lower() in ['tor','lin'] assert Lplot.lower() in ['tot','in'] @@ -496,8 +506,8 @@ def LOS_CrossProj(VType, Ds, us, kIns, kOuts, kRMins, msg += " - provided Ds.shape: %s\n"%str(Ds.shape) msg += " - provided us.shape: %s"%str(us.shape) raise Exception(msg) - lc = [kIns.shape == kOuts.shape, - kIns.shape[-1] == Ds.shape[-1], kIns.ndim == Ds.ndim-1] + lc = [kIns.shape == (Ds.shape[-1],), + kOuts.shape[-1] == Ds.shape[-1], kOuts.ndim == Ds.ndim-1] if not all(lc): msg = "kIns and kOuts must have the same shape and ndim = Ds.ndim-1:\n" msg += " - Ds.shape : %s\n"%str(Ds.shape) @@ -509,18 +519,17 @@ def LOS_CrossProj(VType, Ds, us, kIns, kOuts, kRMins, if Ds.ndim == 2: Ds = np.array([Ds]) us = np.array([us]) - kIns = np.array([kIns]) kOuts = np.array([kOuts]) nref, _, nL = Ds.shape - k0 = kIns if Lplot.lower() == 'in' else np.zeros((nL,)) + #k0 = kIns if Lplot.lower() == 'in' else np.zeros((nL,)) if VType.lower()=='tor' and proj in ['cross','all']: CrossProjAng = np.arccos(np.sqrt(us[:,0,:]**2+us[:,1,:]**2) /np.sqrt(np.sum(us**2,axis=0))) nkp = np.ceil(25.*(1 - (CrossProjAng/(np.pi/4)-1)**2) + 2) - ks = np.max([kRMins,kIns],axis=0) if Lplot.lower()=='in' else kRMins + #ks = np.max([kRMins,kIns],axis=0) if Lplot.lower()=='in' else kRMins pts0 = [] if multi: for ii in range(0,nL): diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 24ab5b729..bd19024ad 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -3841,7 +3841,7 @@ def get_reflections_as_cam(self, Type=None, Name=None, nb=None): return lcam, Types - def add_reflections(self, Type=None, nb=None, coefs=None): + def add_reflections(self, Type=None, nb=None): """ Add relfected LOS to the camera Reflected LOS can be of 3 types: @@ -3859,42 +3859,37 @@ def add_reflections(self, Type=None, nb=None, coefs=None): nb = 1 nb = int(nb) assert nb > 0 - if coefs is not None: - coefs = np.atleast_1d(coefs).astype(float).ravel() - if coefs.size == 1: - coefs = np.repeat(coefs, self.nRays) - assert coefs.shape == (self.nRays,) - assert np.all(coefs >= 0.) and np.all(coefs <= 1.) # Prepare output nRays = self.nRays - Types = np.full((nb,nRays), 0, dtype=int) - us = np.full((nb,3,nRays), np.nan, dtype=float) - kouts = np.full((nb,nRays), np.nan, dtype=float) - indouts = np.full((nb,3,nRays), 0, dtype=int) - vperps = np.full((nb,3,nRays), np.nan, dtype=float) + Types = np.full((Rays,nb), 0, dtype=int) + Ds = np.full((3,nRays,nb), np.nan, dtype=float) + us = np.full((3,nRays,nb), np.nan, dtype=float) + kouts = np.full((nRays,nb), np.nan, dtype=float) + indouts = np.full((3,nRays,nb), 0, dtype=int) + vperps = np.full((3,nRays,nb), np.nan, dtype=float) # Run first iteration - Ds = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u - us[0,:,:], Types[0,:] = self.config._reflect_geom(self.u, + Ds[:,:,0] = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u + us[:,:,0], Types[:,0] = self.config._reflect_geom(self.u, self._dgeom['vperp'], Type=Type) - largs, dkwd = self._prepare_inputs_kInOut(D=Ds, u=us[0,:,:]) + largs, dkwd = self._prepare_inputs_kInOut(D=Ds[:,:,0], u=us[:,:,0]) outi = _GG.LOS_Calc_PInOut_VesStruct(*largs, **dkwd)[1:] - kouts[0,:], vperps[0,:,:], indouts[0,:,:] = outi + kouts[:,0], vperps[:,:,0], indouts[:,:,0] = outi # Run other iterations for ii in range(1,nb): - Ds = Ds + (kouts[ii-1:ii,:]-1.e-12) * us[ii-1,:,:] - us[ii,:,:], Types[ii,:] = self.config._reflect_geom(us[ii-1,:,:], - vperps[ii-1,:,:], + Ds[:,:,ii] = Ds[:,:,ii-1] + (kouts[:,ii-1:ii]-1.e-12) * us[:,:,ii-1] + us[:,:,ii], Types[:,ii] = self.config._reflect_geom(us[:,:,ii-1], + vperps[:,:,ii-1], Type=Type) - outi = _GG.LOS_Calc_PInOut_VesStruct(Ds, us[ii,:,:], + outi = _GG.LOS_Calc_PInOut_VesStruct(Ds[:,:,ii], us[:,:,ii], *largs[2:], **dkwd)[1:] - kouts[ii,:], vperps[ii,:,:], indouts[ii,:,:] = outi + kouts[:,ii], vperps[:,:,ii], indouts[:,:,ii] = outi self._dgeom['dreflect'] = {'nb':nb, 'Type':Type, 'Types':Types, - 'us':us, 'kouts':kouts, 'indouts':indouts} + 'Ds':Ds, 'us':us, 'kouts':kouts, 'indouts':indouts} @@ -4358,30 +4353,44 @@ def get_subset(self, indch=None, Name=None): obj = self.__class__(fromdict=dd) return obj - def _get_plotL(self, reflections=True, Lplot='Tot', proj='All', ind=None, multi=False): + def _get_plotL(self, reflections=True, Lplot='Tot', + proj='All', ind=None, multi=False): """ Get the (R,Z) coordinates of the cross-section projections """ ind = self._check_indch(ind) - if ind.size>0: - Ds, us = self.D[:,ind], self.u[:,ind] - if ind.size==1: - Ds, us = Ds.reshape((3,1)), us.reshape((3,1)) - kPIn, kPOut = self.kIn[ind], self.kOut[ind] - if self.config.Id.Type=='Tor': - kRMin = self._dgeom['kRMin'][ind] + if ind.size > 0: + if Lplot.lower() == 'tot': + Ds = self.D[:,ind] else: - kRMin = None + Ds = self.D[:,ind] + self.kIn[None,ind] * self.u[:,ind] + kOuts = np.atleast_1d(self.kOut[ind])[:,None] + if ind.size == 1: + Ds, us = Ds[:,None], us[:,None] + Ds, us = Ds[:,:,None], us[:,:,None] + kRMin = None # Add reflections ? c0 = (reflections and self._dgeom.get('dreflect') is not None and self._dgeom['dreflect'].get('us') is not None) if c0: - pts = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, - kPIn, kPOut, kRMin, proj=proj, - Lplot=Lplot, multi=multi) - else: - pts = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, - kPIn, kPOut, kRMin, proj=proj, - Lplot=Lplot, multi=multi) + Dsadd = self._dgeom['dreflect']['Ds'][:,ind,:] + usadd = self._dgeom['dreflect']['us'][:,ind,:] + kOutsadd = self._dgeom['dreflect']['kouts'][ind,:] + if ind.size == 1: + Dsadd, usadd = Dsadd[:,None,:], usadd[:,None,:] + kOutsadd = kOutsadd[:,None,:] + Ds = np.concatenate((Ds, Dsadd), axis=-1) + us = np.concatenate((us, usadd), axis=-1) + kOuts = np.concatenate((kOuts, kOutsadd), axis=-1) + if self.config.Id.Type == 'Tor': + kRMin = _comp.LOS_PRMin(Ds, us, kPOut=kOuts, + Eps=1.e-12, squeeze=False) + + elif self.config.Id.Type == 'Tor': + kRMin = self._dgeom['kRMin'][ind][:,None] + + pts = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, + kOuts, kRMin, proj=proj, + Lplot=Lplot, multi=multi) else: pts = None return pts diff --git a/tofu/version.py b/tofu/version.py index a988e6bb6..781536c89 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.0-601-g99a6895' \ No newline at end of file +__version__='1.4.0-602-gcc321d1' \ No newline at end of file From 9629ab2307aea1c286129f05afcee4f16660f354 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Tue, 27 Aug 2019 17:17:27 +0200 Subject: [PATCH 10/19] [Issue 158] tf.geom._plot.Rays_plot() being updated and tested, but bug found in tofu/geom/_sampling_tools.pyx => opened issue #160 --- tofu/geom/_comp.py | 160 +++++++++++++++++++++------------------------ tofu/geom/_core.py | 37 ++++++----- tofu/geom/_plot.py | 33 ++++++---- tofu/version.py | 2 +- 4 files changed, 114 insertions(+), 118 deletions(-) diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index 315d1dab2..b59463882 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -449,7 +449,7 @@ def _get_phithetaproj_dist(poly_closed, ax, Dtheta, nDtheta, ############################################################################### """ -def LOS_PRMin(Ds, us, kPOut=None, Eps=1.e-12, squeeze=True, Test=True): +def LOS_PRMin(Ds, us, kOut=None, Eps=1.e-12, squeeze=True, Test=True): """ Compute the point on the LOS where the major radius is minimum """ if Test: assert Ds.ndim in [1,2,3] and 3 in Ds.shape and Ds.shape == us.shape @@ -463,8 +463,9 @@ def LOS_PRMin(Ds, us, kPOut=None, Eps=1.e-12, squeeze=True, Test=True): elif Ds.ndim == 2: Ds, us = Ds[:,:,None], us[:,:,None] if kOut is not None: - kOut = kOut[:,None] - nref, _, nlos = Ds.shape + if kOut.ndim == 1: + kOut = kOut[:,None] + _, nlos, nref = Ds.shape kRMin = np.full((nlos,nref), np.nan) uparN = np.sqrt(us[0,:,:]**2 + us[1,:,:]**2) @@ -492,107 +493,94 @@ def LOS_PRMin(Ds, us, kPOut=None, Eps=1.e-12, squeeze=True, Test=True): return kRMin -def LOS_CrossProj(VType, Ds, us, kIns, kOuts, kRMins, - proj='All', multi=False): +def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, + num_threads=16, Test=True): """ Compute the parameters to plot the poloidal projection of the LOS """ assert type(VType) is str and VType.lower() in ['tor','lin'] - assert Lplot.lower() in ['tot','in'] - assert type(proj) is str - proj = proj.lower() - assert proj in ['cross','hor','all','3d'] - lc = [Ds.ndim in [2,3], Ds.shape == us.shape] + dproj = {'cross':('R','Z'), 'hor':('x,y'), 'all':('R','Z','x','y'), + '3d':('x','y','z')} + assert type(proj) in [str, tuple] + if type(proj) is tuple: + assert all([type(pp) is str for pp in proj]) + lcoords = proj + else: + proj = proj.lower() + assert proj in dproj.keys() + lcoords = dproj[proj] + + lc = [Ds.ndim == 3, Ds.shape == us.shape] if not all(lc): msg = "Ds and us must have the same shape and dim in [2,3]:\n" msg += " - provided Ds.shape: %s\n"%str(Ds.shape) msg += " - provided us.shape: %s"%str(us.shape) raise Exception(msg) - lc = [kIns.shape == (Ds.shape[-1],), - kOuts.shape[-1] == Ds.shape[-1], kOuts.ndim == Ds.ndim-1] + lc = [kOuts.size == Ds.size/3, kOuts.shape == Ds.shape[1:]] if not all(lc): - msg = "kIns and kOuts must have the same shape and ndim = Ds.ndim-1:\n" + msg = "kOuts must have the same shape and ndim = Ds.ndim-1:\n" msg += " - Ds.shape : %s\n"%str(Ds.shape) - msg += " - kIns.shape : %s\n"%str(kIns.shape) msg += " - kOutss.shape: %s"%str(kOuts.shape) raise Exception(msg) # Prepare inputs - if Ds.ndim == 2: - Ds = np.array([Ds]) - us = np.array([us]) - kOuts = np.array([kOuts]) - - nref, _, nL = Ds.shape - - #k0 = kIns if Lplot.lower() == 'in' else np.zeros((nL,)) + _, nlos, nseg = Ds.shape + + # Detailed sampling for 'tor' and ('cross' or 'all') + R, Z = None, None + if 'R' in lcoords or 'Z' in lcoords: + angcross = np.arccos(np.sqrt(us[0,...]**2 + us[1,...]**2) + /np.sqrt(np.sum(us**2, axis=0))) + resnk = np.ceil(25.*(1 - (angcross/(np.pi/4)-1)**2) + 5) + resnk = 1./resnk.ravel() + + # Use optimized get sample + DL = np.vstack((np.zeros((nlos*nseg,),dtype=float), kOuts.ravel())) + k, reseff, lind = _GG.LOS_get_sample(nlos*nseg, resnk, DL, + dmethod='rel', method='simps', + num_threads=num_threads, Test=Test) + + assert lind.size == nseg*nlos - 1 + ind = lind[::nseg] + nbrep = np.r_[lind[0], np.diff(lind), k.size - lind[-1]] + pts = (np.repeat(Ds.reshape((3,nlos*nseg)), nbrep, axis=1) + + k[None,:] * np.repeat(us.reshape((3,nlos*nseg)), nbrep, + axis=1)) + import ipdb # DB + ipdb.set_trace() # DB - if VType.lower()=='tor' and proj in ['cross','all']: - CrossProjAng = np.arccos(np.sqrt(us[:,0,:]**2+us[:,1,:]**2) - /np.sqrt(np.sum(us**2,axis=0))) - nkp = np.ceil(25.*(1 - (CrossProjAng/(np.pi/4)-1)**2) + 2) - #ks = np.max([kRMins,kIns],axis=0) if Lplot.lower()=='in' else kRMins - pts0 = [] if multi: - for ii in range(0,nL): - for jj in range(0,nref): - if np.isnan(kOuts[jj,ii]): - pts0.append( np.array([[np.nan,np.nan], - [np.nan,np.nan]]) ) - else: - k = np.linspace(k0[ii], kOuts[jj,ii], nkp[ii], endpoint=True) - k = np.unique(np.append(k,ks[ii])) - pp = Ds[:,ii:ii+1] + k[np.newaxis,:]*us[:,ii:ii+1] - pts0.append( np.array([np.hypot(pp[0,:],pp[1,:]),pp[2,:]]) ) - else: - for ii in range(0,nL): - if np.isnan(kOuts[ii]): - pts0.append(np.array([[np.nan,np.nan,np.nan], - [np.nan,np.nan,np.nan], - [np.nan,np.nan,np.nan]])) - else: - k = np.linspace(k0[ii],kOuts[ii],nkp[ii],endpoint=True) - k = np.append(np.unique(np.append(k,ks[ii])),np.nan) - pts0.append( Ds[:,ii:ii+1] + k[np.newaxis,:]*us[:,ii:ii+1] ) - pts0 = np.concatenate(tuple(pts0),axis=1) - pts0 = np.array([np.hypot(pts0[0,:],pts0[1,:]),pts0[2,:]]) - - if not (VType.lower()=='tor' and proj=='cross'): - pts = [] - if multi: - for ii in range(0,nL): - if np.isnan(kOuts[ii]): - pts.append( np.array([[np.nan,np.nan], - [np.nan,np.nan], - [np.nan,np.nan]]) ) - else: - k = np.array([k0[ii],kOuts[ii]]) - pts.append( Ds[:,ii:ii+1] + k[np.newaxis,:]*us[:,ii:ii+1] ) - else: - for ii in range(0,nL): - if np.isnan(kOuts[ii]): - pts.append(np.array([[np.nan,np.nan,np.nan], - [np.nan,np.nan,np.nan], - [np.nan,np.nan,np.nan]])) - else: - k = np.array([k0[ii],kOuts[ii],np.nan]) - pts.append( Ds[:,ii:ii+1] + k[np.newaxis,:]*us[:,ii:ii+1] ) - pts = np.concatenate(tuple(pts),axis=1) - - if proj=='hor': - pts = [pp[:2,:] for pp in pts] if multi else pts[:2,:] - elif proj=='cross': - if VType.lower()=='tor': - pts = pts0 + if 'R' in lcoords: + R = np.hypot(pts[0,:],pts[1,:]).split(ind) + if 'Z' in lcoords: + Z = pts[2,:].split(ind) else: - pts = [pp[1:,:] for pp in pts] if multi else pts[1:,:] - elif proj=='all': + if 'R' in lcoords: + R = np.insert(np.hypot(pts[0,:],pts[1,:]), ind, np.nan) + if 'Z' in lcoords: + Z = np.insert(pts[2,:], ind, np.nan) + + # Normal sampling => pts + # unnecessary only if 'tor' and 'cross' + x, y, z = None, None, None + if 'x' in lcoords or 'y' in lcoords or 'z' in lcoords: + Ds = np.concatenate((Ds, Ds[:,:,-1:] + kOuts[None,:,-1:]*us[:,:,-1:]), + axis=-1) if multi: - if VType.lower()=='tor': - pts = [(p0,pp[:2,:]) for (p0,pp) in zip(*[pts0,pts])] - else: - pts = (pts[1:,:],pts[:2,:]) + if 'x' in lcoords: + x = Ds[0,...] + if 'y' in lcoords: + y = Ds[1,...] + if 'z' in lcoords: + z = Ds[2,...] else: - pts = (pts0,pts[:2,:]) if VType.lower()=='tor' else (pts[1:,:],pts[:2,:]) - return pts + nancoords = np.full((3,nlos,1), np.nan) + Ds = np.concatenate((Ds,nancoords), axis=-1) + if 'x' in lcoords: + x = Ds[0,...].ravel() + if 'y' in lcoords: + y = Ds[1,...].ravel() + if 'z' in lcoords: + z = Ds[2,...].ravel() + return R, Z, x, y, z diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index bd19024ad..c5185aae8 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -3646,14 +3646,16 @@ def compute_dgeom(self, extra=True, plotdebug=True): def _compute_dgeom_kRMin(self): # Get RMin if Type is Tor if self.config.Id.Type=='Tor': - kRMin = _comp.LOS_PRMin(self.D, self.u, kPOut=self.kOut, Eps=1.e-12) + kRMin = np.atleast_1d(_comp.LOS_PRMin(self.D, self.u, + kOut=self.kOut, + Eps=1.e-12, squeeze=True)) else: kRMin = None self._dgeom.update({'kRMin':kRMin}) def _compute_dgeom_extra1(self): if self._dgeom['kRMin'] is not None: - PRMin = self.D + self._dgeom['kRMin'][np.newaxis,:]*self.u + PRMin = self.D + self._dgeom['kRMin'][None,:]*self.u RMin = np.hypot(PRMin[0,:],PRMin[1,:]) else: PRMin, RMin = None, None @@ -3862,7 +3864,7 @@ def add_reflections(self, Type=None, nb=None): # Prepare output nRays = self.nRays - Types = np.full((Rays,nb), 0, dtype=int) + Types = np.full((nRays,nb), 0, dtype=int) Ds = np.full((3,nRays,nb), np.nan, dtype=float) us = np.full((3,nRays,nb), np.nan, dtype=float) kouts = np.full((nRays,nb), np.nan, dtype=float) @@ -3880,13 +3882,14 @@ def add_reflections(self, Type=None, nb=None): # Run other iterations for ii in range(1,nb): - Ds[:,:,ii] = Ds[:,:,ii-1] + (kouts[:,ii-1:ii]-1.e-12) * us[:,:,ii-1] - us[:,:,ii], Types[:,ii] = self.config._reflect_geom(us[:,:,ii-1], - vperps[:,:,ii-1], - Type=Type) - outi = _GG.LOS_Calc_PInOut_VesStruct(Ds[:,:,ii], us[:,:,ii], + Dsi = Ds[:,:,ii-1] + (kouts[None,:,ii-1]-1.e-12) * us[:,:,ii-1] + usi, Types[:,ii] = self.config._reflect_geom(us[:,:,ii-1], + vperps[:,:,ii-1], + Type=Type) + outi = _GG.LOS_Calc_PInOut_VesStruct(Dsi, usi, *largs[2:], **dkwd)[1:] kouts[:,ii], vperps[:,:,ii], indouts[:,:,ii] = outi + Ds[:,:,ii], us[:,:,ii] = Dsi, usi self._dgeom['dreflect'] = {'nb':nb, 'Type':Type, 'Types':Types, 'Ds':Ds, 'us':us, 'kouts':kouts, 'indouts':indouts} @@ -4358,11 +4361,12 @@ def _get_plotL(self, reflections=True, Lplot='Tot', """ Get the (R,Z) coordinates of the cross-section projections """ ind = self._check_indch(ind) if ind.size > 0: + us = self.u[:,ind] if Lplot.lower() == 'tot': Ds = self.D[:,ind] else: - Ds = self.D[:,ind] + self.kIn[None,ind] * self.u[:,ind] - kOuts = np.atleast_1d(self.kOut[ind])[:,None] + Ds = self.D[:,ind] + self.kIn[None,ind] * us + kOuts = np.atleast_1d(self.kOut[ind])[:,None] if ind.size == 1: Ds, us = Ds[:,None], us[:,None] Ds, us = Ds[:,:,None], us[:,:,None] @@ -4377,23 +4381,22 @@ def _get_plotL(self, reflections=True, Lplot='Tot', kOutsadd = self._dgeom['dreflect']['kouts'][ind,:] if ind.size == 1: Dsadd, usadd = Dsadd[:,None,:], usadd[:,None,:] - kOutsadd = kOutsadd[:,None,:] + kOutsadd = kOutsadd[None,:] Ds = np.concatenate((Ds, Dsadd), axis=-1) us = np.concatenate((us, usadd), axis=-1) kOuts = np.concatenate((kOuts, kOutsadd), axis=-1) if self.config.Id.Type == 'Tor': - kRMin = _comp.LOS_PRMin(Ds, us, kPOut=kOuts, + kRMin = _comp.LOS_PRMin(Ds, us, kOut=kOuts, Eps=1.e-12, squeeze=False) elif self.config.Id.Type == 'Tor': kRMin = self._dgeom['kRMin'][ind][:,None] - pts = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, - kOuts, kRMin, proj=proj, - Lplot=Lplot, multi=multi) + R, Z, x, y, z = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, + kOuts, proj=proj, multi=multi) else: - pts = None - return pts + R, Z, x, y, z = None, None, None, None, None + return R, Z, x, y, z def get_sample(self, res=None, resMode='abs', DL=None, method='sum', ind=None, pts=False, compact=True, num_threads=_NUM_THREADS, Test=True): diff --git a/tofu/geom/_plot.py b/tofu/geom/_plot.py index 3657c6395..ef0442a93 100644 --- a/tofu/geom/_plot.py +++ b/tofu/geom/_plot.py @@ -972,12 +972,13 @@ def _Rays_plot_Cross(L,Leg=None, reflections=True, Type=L.Ves.Type) if 'L' in Elt: - pts = L._get_plotL(Lplot=Lplot, proj='Cross', ind=ind, multi=multi) + R, Z, _, _, _ = L._get_plotL(Lplot=Lplot, proj='Cross', + reflections=reflections, ind=ind, multi=multi) if multi: - for ii in range(0,len(pts)): - ax.plot(pts[ii][0,:], pts[ii][1,:], label=Leg[ii], **dL) + for ii in range(0,len(R)): + ax.plot(R[ii], Z[ii], label=Leg[ii], **dL) else: - ax.plot(pts[0,:], pts[1,:], label=Leg, **dL) + ax.plot(R, Z, label=Leg, **dL) for kk in dPts.keys(): if kk in Elt: @@ -1008,7 +1009,8 @@ def _Rays_plot_Cross(L,Leg=None, reflections=True, -def _Rays_plot_Hor(L, Leg=None, Lplot='Tot', Elt='LDIORP',ax=None, +def _Rays_plot_Hor(L, Leg=None, reflections=True, + Lplot='Tot', Elt='LDIORP',ax=None, dL=_def.LOSLd, dPtD=_def.LOSMd, dPtI=_def.LOSMd, dPtO=_def.LOSMd, dPtR=_def.LOSMd, dPtP=_def.LOSMd, dLeg=_def.TorLegd, multi=False, ind=None, @@ -1021,12 +1023,13 @@ def _Rays_plot_Hor(L, Leg=None, Lplot='Tot', Elt='LDIORP',ax=None, ax = _def.Plot_LOSProj_DefAxes('Hor', fs=fs, wintit=wintit, Type=L.Ves.Type) if 'L' in Elt: - pts = L._get_plotL(Lplot=Lplot, proj='hor', ind=ind, multi=multi) + _, _, x, y, _ = L._get_plotL(Lplot=Lplot, proj='hor', + reflections=reflections, ind=ind, multi=multi) if multi: - for ii in range(0,len(pts)): - ax.plot(pts[ii][0,:], pts[ii][1,:], label=Leg[ii], **dL) + for ii in range(0,len(x)): + ax.plot(x[ii], y[ii], label=Leg[ii], **dL) else: - ax.plot(pts[0,:], pts[1,:], label=Leg, **dL) + ax.plot(x, y, label=Leg, **dL) for kk in dPts.keys(): if kk in Elt: @@ -1053,7 +1056,8 @@ def _Rays_plot_Hor(L, Leg=None, Lplot='Tot', Elt='LDIORP',ax=None, -def _Rays_plot_3D(L,Leg=None,Lplot='Tot',Elt='LDIORr',ax=None, +def _Rays_plot_3D(L, Leg=None, reflections=True, + Lplot='Tot',Elt='LDIORr',ax=None, dL=_def.LOSLd, dPtD=_def.LOSMd, dPtI=_def.LOSMd, dPtO=_def.LOSMd, dPtR=_def.LOSMd, dPtP=_def.LOSMd, dLeg=_def.TorLegd, multi=False, ind=None, @@ -1066,13 +1070,14 @@ def _Rays_plot_3D(L,Leg=None,Lplot='Tot',Elt='LDIORr',ax=None, ax = _def.Plot_3D_plt_Tor_DefAxes(fs=fs, wintit=wintit) if 'L' in Elt: - pts = L._get_plotL(Lplot=Lplot, proj='3d', ind=ind, multi=multi) + _, _, x, y, z = L._get_plotL(Lplot=Lplot, proj='3d', + reflections=reflections, ind=ind, multi=multi) if multi: - for ii in range(0,len(pts)): - ax.plot(pts[ii][0,:], pts[ii][1,:], pts[ii][2,:], + for ii in range(0,len(x)): + ax.plot(x[ii], y[ii], z[ii], label=Leg[ii], **dL) else: - ax.plot(pts[0,:], pts[1,:], pts[2,:], label=Leg, **dL) + ax.plot(x, y, z, label=Leg, **dL) for kk in dPts.keys(): if kk in Elt: diff --git a/tofu/version.py b/tofu/version.py index 781536c89..2b314b04b 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.0-602-gcc321d1' \ No newline at end of file +__version__='1.4.0-604-g88b4711' \ No newline at end of file From 10aa10f763c7c3959c9fcb7844af68669cd8e34e Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 28 Aug 2019 07:42:59 +0200 Subject: [PATCH 11/19] [Issue 158] Debugging reflections --- tofu/geom/_comp.py | 41 ++++++++++++++++++++++++++--------------- tofu/geom/_plot.py | 6 ++++-- tofu/version.py | 2 +- 3 files changed, 31 insertions(+), 18 deletions(-) diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index b59463882..9b30a305d 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -494,7 +494,7 @@ def LOS_PRMin(Ds, us, kOut=None, Eps=1.e-12, squeeze=True, Test=True): def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, - num_threads=16, Test=True): + num_threads=16, return_pts=False, Test=True): """ Compute the parameters to plot the poloidal projection of the LOS """ assert type(VType) is str and VType.lower() in ['tor','lin'] dproj = {'cross':('R','Z'), 'hor':('x,y'), 'all':('R','Z','x','y'), @@ -507,6 +507,8 @@ def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, proj = proj.lower() assert proj in dproj.keys() lcoords = dproj[proj] + if return_pts: + assert proj in ['cross','hor', '3d'] lc = [Ds.ndim == 3, Ds.shape == us.shape] if not all(lc): @@ -534,8 +536,9 @@ def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, # Use optimized get sample DL = np.vstack((np.zeros((nlos*nseg,),dtype=float), kOuts.ravel())) + # TBF: should be 'rel', but 'abs' for debugging, until Issue160 fixed k, reseff, lind = _GG.LOS_get_sample(nlos*nseg, resnk, DL, - dmethod='rel', method='simps', + dmethod='abs', method='simps', num_threads=num_threads, Test=Test) assert lind.size == nseg*nlos - 1 @@ -544,19 +547,23 @@ def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, pts = (np.repeat(Ds.reshape((3,nlos*nseg)), nbrep, axis=1) + k[None,:] * np.repeat(us.reshape((3,nlos*nseg)), nbrep, axis=1)) - import ipdb # DB - ipdb.set_trace() # DB - - if multi: - if 'R' in lcoords: - R = np.hypot(pts[0,:],pts[1,:]).split(ind) - if 'Z' in lcoords: - Z = pts[2,:].split(ind) + if return_pts: + pts = np.array([np.hypot(pts[0,:],pts[1,:]), pts[2,:]]) + if multi: + pts = pts.split(ind) + else: + pts = np.insert(pts, ind, np.nan) else: - if 'R' in lcoords: - R = np.insert(np.hypot(pts[0,:],pts[1,:]), ind, np.nan) - if 'Z' in lcoords: - Z = np.insert(pts[2,:], ind, np.nan) + if multi: + if 'R' in lcoords: + R = np.hypot(pts[0,:],pts[1,:]).split(ind) + if 'Z' in lcoords: + Z = pts[2,:].split(ind) + else: + if 'R' in lcoords: + R = np.insert(np.hypot(pts[0,:],pts[1,:]), ind, np.nan) + if 'Z' in lcoords: + Z = np.insert(pts[2,:], ind, np.nan) # Normal sampling => pts # unnecessary only if 'tor' and 'cross' @@ -580,7 +587,11 @@ def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, y = Ds[1,...].ravel() if 'z' in lcoords: z = Ds[2,...].ravel() - return R, Z, x, y, z + + if return_pts: + return pts + else: + return R, Z, x, y, z diff --git a/tofu/geom/_plot.py b/tofu/geom/_plot.py index ef0442a93..2f4ee91fc 100644 --- a/tofu/geom/_plot.py +++ b/tofu/geom/_plot.py @@ -1489,8 +1489,10 @@ def _Cam12D_plottouch(cam, key=None, ind=None, quant='lengths', nchMax=_nchMax, dax['cross'][0], dax['hor'][0] = out if cam._isLOS(): - lCross = cam._get_plotL(Lplot=Lplot, proj='cross', multi=True) - lHor = cam._get_plotL(Lplot=Lplot, proj='hor', multi=True) + # lCross = cam._get_plotL(Lplot=Lplot, proj='cross', multi=True) + lR, lZ = cam._get_plotL(Lplot=Lplot, proj='cross', multi=True)[:2] + # lHor = cam._get_plotL(Lplot=Lplot, proj='hor', multi=True) + lx, ly = cam._get_plotL(Lplot=Lplot, proj='hor', multi=True)[2:4] if Bck and nD == 2: crossbck = [lCross[indbck[0]],nan2,lCross[indbck[1]],nan2, lCross[indbck[2]],nan2,lCross[indbck[3]]] diff --git a/tofu/version.py b/tofu/version.py index 2b314b04b..cb83b26f1 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.0-604-g88b4711' \ No newline at end of file +__version__='1.4.0-617-g944cb49' \ No newline at end of file From 76c5f9ae8031a31ea8e1b68af797522de1567164 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 28 Aug 2019 09:20:24 +0200 Subject: [PATCH 12/19] [Issue 158] Rays.plot() and Rays .plot_touch() operational with reflections for all Cams. TODO: calc_signal() --- tofu/geom/_comp.py | 59 +++++++++++++++++++++++++++++----------------- tofu/geom/_core.py | 14 ++++++----- tofu/geom/_plot.py | 8 +++---- tofu/version.py | 2 +- 4 files changed, 51 insertions(+), 32 deletions(-) diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index 9b30a305d..ae7a868af 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -536,29 +536,29 @@ def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, # Use optimized get sample DL = np.vstack((np.zeros((nlos*nseg,),dtype=float), kOuts.ravel())) - # TBF: should be 'rel', but 'abs' for debugging, until Issue160 fixed k, reseff, lind = _GG.LOS_get_sample(nlos*nseg, resnk, DL, - dmethod='abs', method='simps', + dmethod='rel', method='simps', num_threads=num_threads, Test=Test) assert lind.size == nseg*nlos - 1 - ind = lind[::nseg] + ind = lind[nseg-1::nseg] nbrep = np.r_[lind[0], np.diff(lind), k.size - lind[-1]] pts = (np.repeat(Ds.reshape((3,nlos*nseg)), nbrep, axis=1) + k[None,:] * np.repeat(us.reshape((3,nlos*nseg)), nbrep, axis=1)) + if return_pts: pts = np.array([np.hypot(pts[0,:],pts[1,:]), pts[2,:]]) if multi: - pts = pts.split(ind) + pts = np.split(pts, ind, axis=1) else: - pts = np.insert(pts, ind, np.nan) + pts = np.insert(pts, ind, np.nan, axis=1) else: if multi: if 'R' in lcoords: - R = np.hypot(pts[0,:],pts[1,:]).split(ind) + R = np.split(np.hypot(pts[0,:],pts[1,:]), ind) if 'Z' in lcoords: - Z = pts[2,:].split(ind) + Z = np.split(pts[2,:], ind) else: if 'R' in lcoords: R = np.insert(np.hypot(pts[0,:],pts[1,:]), ind, np.nan) @@ -569,24 +569,41 @@ def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, # unnecessary only if 'tor' and 'cross' x, y, z = None, None, None if 'x' in lcoords or 'y' in lcoords or 'z' in lcoords: - Ds = np.concatenate((Ds, Ds[:,:,-1:] + kOuts[None,:,-1:]*us[:,:,-1:]), + pts = np.concatenate((Ds, Ds[:,:,-1:] + kOuts[None,:,-1:]*us[:,:,-1:]), axis=-1) + + pts = pts.reshape((3,nlos*(nseg+1))) if multi: - if 'x' in lcoords: - x = Ds[0,...] - if 'y' in lcoords: - y = Ds[1,...] - if 'z' in lcoords: - z = Ds[2,...] + ind = np.arange(1,nlos)*(nseg+1) else: nancoords = np.full((3,nlos,1), np.nan) - Ds = np.concatenate((Ds,nancoords), axis=-1) - if 'x' in lcoords: - x = Ds[0,...].ravel() - if 'y' in lcoords: - y = Ds[1,...].ravel() - if 'z' in lcoords: - z = Ds[2,...].ravel() + pts = np.concatenate((pts,nancoords), axis=-1) + + if return_pts: + assert proj in ['hor','3d'] + if multi: + if proj == 'hor': + pts = np.split(pts[:2,:], ind, axis=1) + else: + pts = np.split(pts, ind, axis=1) + elif proj == 'hor': + pts = pts[:2,:] + + else: + if multi: + if 'x' in lcoords: + x = np.split(pts[0,:], ind) + if 'y' in lcoords: + y = np.split(pts[1,:], ind) + if 'z' in lcoords: + z = np.split(pts[2,:], ind) + else: + if 'x' in lcoords: + x = pts[0,:] + if 'y' in lcoords: + y = pts[1,:] + if 'z' in lcoords: + z = pts[2,:] if return_pts: return pts diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index c5185aae8..e3e0976f6 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -4357,16 +4357,17 @@ def get_subset(self, indch=None, Name=None): return obj def _get_plotL(self, reflections=True, Lplot='Tot', - proj='All', ind=None, multi=False): + proj='All', ind=None, return_pts=False, multi=False): """ Get the (R,Z) coordinates of the cross-section projections """ ind = self._check_indch(ind) if ind.size > 0: us = self.u[:,ind] + kOuts = np.atleast_1d(self.kOut[ind])[:,None] if Lplot.lower() == 'tot': Ds = self.D[:,ind] else: Ds = self.D[:,ind] + self.kIn[None,ind] * us - kOuts = np.atleast_1d(self.kOut[ind])[:,None] + kOuts = kOuts - np.atleast_1d(self.kIn[ind])[:,None] if ind.size == 1: Ds, us = Ds[:,None], us[:,None] Ds, us = Ds[:,:,None], us[:,:,None] @@ -4392,11 +4393,12 @@ def _get_plotL(self, reflections=True, Lplot='Tot', elif self.config.Id.Type == 'Tor': kRMin = self._dgeom['kRMin'][ind][:,None] - R, Z, x, y, z = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, - kOuts, proj=proj, multi=multi) + out = _comp.LOS_CrossProj(self.config.Id.Type, Ds, us, + kOuts, proj=proj, + return_pts=return_pts, multi=multi) else: - R, Z, x, y, z = None, None, None, None, None - return R, Z, x, y, z + out = None + return out def get_sample(self, res=None, resMode='abs', DL=None, method='sum', ind=None, pts=False, compact=True, num_threads=_NUM_THREADS, Test=True): diff --git a/tofu/geom/_plot.py b/tofu/geom/_plot.py index 2f4ee91fc..bba7a02e3 100644 --- a/tofu/geom/_plot.py +++ b/tofu/geom/_plot.py @@ -1489,10 +1489,10 @@ def _Cam12D_plottouch(cam, key=None, ind=None, quant='lengths', nchMax=_nchMax, dax['cross'][0], dax['hor'][0] = out if cam._isLOS(): - # lCross = cam._get_plotL(Lplot=Lplot, proj='cross', multi=True) - lR, lZ = cam._get_plotL(Lplot=Lplot, proj='cross', multi=True)[:2] - # lHor = cam._get_plotL(Lplot=Lplot, proj='hor', multi=True) - lx, ly = cam._get_plotL(Lplot=Lplot, proj='hor', multi=True)[2:4] + lCross = cam._get_plotL(Lplot=Lplot, proj='cross', + return_pts=True, multi=True) + lHor = cam._get_plotL(Lplot=Lplot, proj='hor', + return_pts=True, multi=True) if Bck and nD == 2: crossbck = [lCross[indbck[0]],nan2,lCross[indbck[1]],nan2, lCross[indbck[2]],nan2,lCross[indbck[3]]] diff --git a/tofu/version.py b/tofu/version.py index cb83b26f1..cbfe684eb 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.0-617-g944cb49' \ No newline at end of file +__version__='1.4.0-618-g10aa10f' \ No newline at end of file From 4f1df50a6010b615896de8533cb6d8695657ec4e Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 28 Aug 2019 11:07:35 +0200 Subject: [PATCH 13/19] [Issue 158] Added dreflect to Struct for Struct an dsegment-wise handling of reflections. TODO: test / check consistency of indout and debug unit tests for Rays._get_plotlL() --- tofu/geom/_core.py | 117 +++++++++++++++++++++++++++++++++++++++++++-- tofu/version.py | 2 +- 2 files changed, 115 insertions(+), 4 deletions(-) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index e3e0976f6..92599038d 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -122,6 +122,7 @@ class Struct(utils.ToFuObject): 'dgeom':{'Type':'Tor', 'Lim':[], 'arrayorder':'C'}, 'dsino':{}, 'dphys':{}, + 'dreflect':{'Type':'specular'}, 'dmisc':{'color':'k'}} _dplot = {'cross':{'Elt':'P', 'dP':{'color':'k','lw':2}, @@ -141,7 +142,7 @@ class Struct(utils.ToFuObject): 'linewidth':0., 'antialiased':False}, 'Lim':None, 'Nstep':50}} - + _DREFLECT_DTYPES = {'specular':0, 'diffusive':1, 'ccube':2} # Does not exist beofre Python 3.6 !!! def __init_subclass__(cls, color='k', **kwdargs): @@ -184,6 +185,7 @@ def _reset(self): self._dgeom = dict.fromkeys(self._get_keys_dgeom()) self._dsino = dict.fromkeys(self._get_keys_dsino()) self._dphys = dict.fromkeys(self._get_keys_dphys()) + self._dreflect = dict.fromkeys(self._get_keys_dreflect()) self._dmisc = dict.fromkeys(self._get_keys_dmisc()) #self._dplot = copy.deepcopy(self.__class__._ddef['dplot']) @@ -237,6 +239,11 @@ def _get_largs_dphys(): largs = ['lSymbols'] return largs + @staticmethod + def _get_largs_dreflect(): + largs = ['Types', 'coefs'] + return largs + @staticmethod def _get_largs_dmisc(): largs = ['color'] @@ -402,6 +409,25 @@ def _checkformat_inputs_dphys(lSymbols=None): lSymbols = np.asarray(lSymbols,dtype=str) return lSymbols + def _checkformat_inputs_dreflect(self, Types=None, coefs=None): + if Types is None: + Types = self._ddef['dreflect']['Type'] + + assert type(Types) in [str, np.ndarray] + if type(Types) is str: + assert Types in self._DREFLECT_DTYPES.keys() + Types = np.full((self.nseg,), self._DREFLECT_DTYPES[Types], dtype=int) + else: + Types = Types.astype(int).ravel() + assert Types.shape == () + Typesu = np.unique(Types) + lc = np.array([Typesu == vv + for vv in self._DREFLECT_DTYPES.values()]) + assert np.all(np.any(Types, axis=0)) + + assert coefs is None + return Types, coefs + @classmethod def _checkformat_inputs_dmisc(cls, color=None): if color is None: @@ -432,6 +458,11 @@ def _get_keys_dphys(): lk = ['lSymbols'] return lk + @staticmethod + def _get_keys_dreflect(): + lk = ['Types', 'coefs'] + return lk + @staticmethod def _get_keys_dmisc(): lk = ['color'] @@ -450,10 +481,13 @@ def _init(self, Poly=None, Type=_Type, kwdgeom = self._extract_kwdargs(allkwds, largs) largs = self._get_largs_dphys() kwdphys = self._extract_kwdargs(allkwds, largs) + largs = self._get_largs_dreflect() + kwdreflect = self._extract_kwdargs(allkwds, largs) largs = self._get_largs_dmisc() kwdmisc = self._extract_kwdargs(allkwds, largs) self._set_dgeom(**kwdgeom) self.set_dphys(**kwdphys) + self.set_dreflect(**kwdreflect) self._set_dmisc(**kwdmisc) self._dstrip['strip'] = 0 @@ -489,6 +523,11 @@ def set_dphys(self, lSymbols=None): lSymbols = self._checkformat_inputs_dphys(lSymbols) self._dphys['lSymbols'] = lSymbols + def set_dreflect(self, Types=None, coefs=None): + Types, coefs = self._checkformat_inputs_dreflect(Types=Types, coefs=coefs) + self._dreflect['Types'] = Types + self._dreflect['coefs'] = coefs + def _set_color(self, color=None): color = self._checkformat_inputs_dmisc(color=color) self._dmisc['color'] = color @@ -512,6 +551,9 @@ def _strip_dsino(self, lkeep=['RefPt','nP']): def _strip_dphys(self, lkeep=['lSymbols']): utils.ToFuObject._strip_dict(self._dphys, lkeep=lkeep) + def _strip_dreflect(self, lkeep=['lSymbols']): + utils.ToFuObject._strip_dict(self._dreflect, lkeep=lkeep) + def _strip_dmisc(self, lkeep=['color']): utils.ToFuObject._strip_dict(self._dmisc, lkeep=lkeep) @@ -543,6 +585,14 @@ def _rebuild_dphys(self, lkeep=['lSymbols']): lkeep=lkeep, dname='dphys') self.set_dphys(lSymbols=self.dphys['lSymbols']) + def _rebuild_dreflect(self, lkeep=['lSymbols']): + reset = utils.ToFuObject._test_Rebuild(self._dreflect, lkeep=lkeep) + if reset: + utils.ToFuObject._check_Fields4Rebuild(self._dreflect, + lkeep=lkeep, dname='dreflect') + self.set_dreflect(Types=self.dreflect['Types'], + coefs=self.dreflect['coefs']) + def _rebuild_dmisc(self, lkeep=['color']): reset = utils.ToFuObject._test_Rebuild(self._dmisc, lkeep=lkeep) if reset: @@ -560,7 +610,7 @@ def _strip_init(cls): nMax = max(cls._dstrip['allowed']) doc = """ 1: Remove dsino expendables - 2: Remove also dgeom, dphys and dmisc expendables""" + 2: Remove also dgeom, dphys, dreflect and dmisc expendables""" doc = utils.ToFuObjectBase.strip.__doc__.format(doc,nMax) if sys.version[0]=='2': cls.strip.__func__.__doc__ = doc @@ -576,22 +626,26 @@ def _strip(self, strip=0): self._rebuild_dgeom() self._rebuild_dsino() self._rebuild_dphys() + self._rebuild_dreflect() self._rebuild_dmisc() elif strip==1: self._strip_dsino() self._rebuild_dgeom() self._rebuild_dphys() + self._rebuild_dreflect() self._rebuild_dmisc() else: self._strip_dsino() self._strip_dgeom() self._strip_dphys() + self._strip_dreflect() self._strip_dmisc() def _to_dict(self): dout = {'dgeom':{'dict':self.dgeom, 'lexcept':None}, 'dsino':{'dict':self.dsino, 'lexcept':None}, 'dphys':{'dict':self.dphys, 'lexcept':None}, + 'dreflect':{'dict':self.dreflect, 'lexcept':None}, 'dmisc':{'dict':self.dmisc, 'lexcept':None}, 'dplot':{'dict':self._dplot, 'lexcept':None}} return dout @@ -600,6 +654,7 @@ def _from_dict(self, fd): self._dgeom.update(**fd['dgeom']) self._dsino.update(**fd['dsino']) self._dphys.update(**fd['dphys']) + self._dreflect.update(**fd['dreflect']) self._dmisc.update(**fd['dmisc']) if 'dplot' in fd.keys(): self._dplot.update(**fd['dplot']) @@ -624,6 +679,10 @@ def Poly_closed(self): """ Returned the closed polygon """ return np.hstack((self._dgeom['Poly'],self._dgeom['Poly'][:,0:1])) @property + def nseg(self): + """ Retunr the number of segmnents constituting the closed polygon """ + return self._dgeom['Poly'].shape[1] + @property def pos(self): return self._dgeom['pos'] @property @@ -649,6 +708,9 @@ def dsino(self): def dphys(self): return self._dphys @property + def dreflect(self): + return self._dreflect + @property def dmisc(self): return self._dmisc @@ -952,11 +1014,60 @@ def _get_phithetaproj_dist(self, refpt=None, ntheta=None, nphi=None, return dist, nDphi, Dphi, nDtheta, Dtheta + def get_reflections(self, u, vperp, indout2=None, lamb=None): + """ Return the reflected unit vectors from input unit vectors and vperp + + The reflected unit vector depends on the incoming LOS (u), + the local normal unit vector (vperp), and the polygon segment hit + (indout2) + Future releases: dependence on lambda + + Also return per-LOS reflection Types (0:specular, 1:diffusive, 2:ccube) + + """ + + # Get per-LOS reflection Types and associated indices + Types = self._dreflect['Types'][indout2] + indspec = Types == 0 + inddiff = Types == 1 + indcorn = Types == 2 + + # Get reflected unit vectors + u2 = np.full(u.shape, np.nan) + if np.any(np.logical_or(indspec,inddiff)): + vpar = np.array([vperp[1,:]*u[2,:] - vperp[2,:]*u[1,:], + vperp[2,:]*u[0,:] - vperp[0,:]*u[2,:], + vperp[0,:]*u[1,:] - vperp[1,:]*u[0,:]]) + vpar = np.array([vpar[1,:]*vperp[2,:] - vpar[2,:]*vperp[1,:], + vpar[2,:]*vperp[0,:] - vpar[0,:]*vperp[2,:], + vpar[0,:]*vperp[1,:] - vpar[1,:]*vperp[0,:]]) + vpar = vpar / np.sqrt(np.sum(vpar**2, axis=0))[None,:] + + if np.any(indspec): + # Compute u2 for specular + sca = np.sum(u[:,indspec]*vperp[:,indspec],axis=0,keepdims=True) + sca2 = np.sum(u[:,indspec]*vpar[:,indspec],axis=0,keepdims=True) + assert np.all(sca<=0.) and np.all(sca>=-1.) + assert np.all(sca2>=0.) and np.all(sca<=1.) + u2[:,indspec] = - sca*vperp[:,indspec] + sca2*vpar[:,indspec] + + if np.any(inddiff): + # Compute u2 for diffusive + sca = 2.*(np.random.random((1,inddiff.sum()))-0.5) + u2[:,inddiff] = (np.sqrt(1.-sca**2) * vperp[:,inddiff] + + sca * vpar[:,inddiff]) + + if np.any(indcorn): + u2[:,indcorn] = -u[:,indcorn] + return u2, Types + + + def plot(self, lax=None, proj='all', element='PIBsBvV', dP=None, dI=_def.TorId, dBs=_def.TorBsd, dBv=_def.TorBvd, dVect=_def.TorVind, dIHor=_def.TorITord, dBsHor=_def.TorBsTord, dBvHor=_def.TorBvTord, Lim=None, Nstep=_def.TorNTheta, - dLeg=_def.TorLegd, indices=False, + dLeg=_def.TorLegd, indices=True, draw=True, fs=None, wintit=None, Test=True): """ Plot the polygon defining the vessel, in chosen projection diff --git a/tofu/version.py b/tofu/version.py index cbfe684eb..ef58c5bf6 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.0-618-g10aa10f' \ No newline at end of file +__version__='1.4.0-620-g6f448be' \ No newline at end of file From fcae638b43c3865098009a688f0cd83f72f56578 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 28 Aug 2019 12:31:13 +0200 Subject: [PATCH 14/19] [Issue 158] Debugged unit tests and Struct dreflect to be tested with Rays.add_reflections --- tofu/geom/_comp.py | 3 ++- tofu/geom/_core.py | 4 ++-- tofu/version.py | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index ae7a868af..1ee06cabe 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -572,12 +572,13 @@ def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, pts = np.concatenate((Ds, Ds[:,:,-1:] + kOuts[None,:,-1:]*us[:,:,-1:]), axis=-1) - pts = pts.reshape((3,nlos*(nseg+1))) if multi: ind = np.arange(1,nlos)*(nseg+1) + pts = pts.reshape((3,nlos*(nseg+1))) else: nancoords = np.full((3,nlos,1), np.nan) pts = np.concatenate((pts,nancoords), axis=-1) + pts = pts.reshape((3,nlos*(nseg+2))) if return_pts: assert proj in ['hor','3d'] diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 92599038d..cfe759094 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -551,7 +551,7 @@ def _strip_dsino(self, lkeep=['RefPt','nP']): def _strip_dphys(self, lkeep=['lSymbols']): utils.ToFuObject._strip_dict(self._dphys, lkeep=lkeep) - def _strip_dreflect(self, lkeep=['lSymbols']): + def _strip_dreflect(self, lkeep=['Types','coefs']): utils.ToFuObject._strip_dict(self._dreflect, lkeep=lkeep) def _strip_dmisc(self, lkeep=['color']): @@ -585,7 +585,7 @@ def _rebuild_dphys(self, lkeep=['lSymbols']): lkeep=lkeep, dname='dphys') self.set_dphys(lSymbols=self.dphys['lSymbols']) - def _rebuild_dreflect(self, lkeep=['lSymbols']): + def _rebuild_dreflect(self, lkeep=['Types','coefs']): reset = utils.ToFuObject._test_Rebuild(self._dreflect, lkeep=lkeep) if reset: utils.ToFuObject._check_Fields4Rebuild(self._dreflect, diff --git a/tofu/version.py b/tofu/version.py index ef58c5bf6..45d35422e 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.0-620-g6f448be' \ No newline at end of file +__version__='1.4.0-620-g4f1df50' \ No newline at end of file From 89af1099d1a975b1b1830d1f5dacc500567ff22c Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 28 Aug 2019 12:31:48 +0200 Subject: [PATCH 15/19] [Issue 158] Overloaded operator __repr__ for amll tofu objects --- tofu/data/_core.py | 3 ++- tofu/data/_plot.py | 30 ++++++++++++++++++++---------- tofu/imas2tofu/_core.py | 5 +++++ tofu/utils.py | 16 ++++++++++------ tofu/version.py | 2 +- 5 files changed, 38 insertions(+), 18 deletions(-) diff --git a/tofu/data/_core.py b/tofu/data/_core.py index 6f600c83c..e587bf19d 100644 --- a/tofu/data/_core.py +++ b/tofu/data/_core.py @@ -2956,7 +2956,7 @@ def _get_keyingroup(self, str_, group=None, msgstr=None, raise_=False): if msg is not None: msg += "\n\nRequested %s could not be identified !\n"%msgstr msg += "Please provide a valid (unique) key/name/quant/dim:\n\n" - msg += '\n\n'.join(self.get_summary(verb=False, return_='msg')) + msg += self.get_summary(verb=False, return_='msg') if raise_: raise Exception(msg) return key, msg @@ -3002,6 +3002,7 @@ def get_summary(self, sep=' ', line='-', just='l', sep=sep, line=line, table_sep=table_sep, verb=verb, return_=return_) + #--------------------- # Methods for adding ref / quantities #--------------------- diff --git a/tofu/data/_plot.py b/tofu/data/_plot.py index 9a1df4264..67a885884 100644 --- a/tofu/data/_plot.py +++ b/tofu/data/_plot.py @@ -562,8 +562,10 @@ def _DataCam12D_plot(lData, key=None, nchMax=_nchMax, ntMax=_ntMax, if c1 and 'LOS' in lData[0]._dgeom['lCam'][0].Id.Cls: lCross, lHor, llab = [], [], [] for cc in lData[0]._dgeom['lCam']: - lCross += cc._get_plotL(Lplot=Lplot, proj='cross', multi=True) - lHor += cc._get_plotL(Lplot=Lplot, proj='hor', multi=True) + lCross += cc._get_plotL(Lplot=Lplot, proj='cross', + return_pts=True, multi=True) + lHor += cc._get_plotL(Lplot=Lplot, proj='hor', + return_pts=True, multi=True) if bck and nD == 2: crossbck = [lCross[indbck[0]],nan2,lCross[indbck[1]],nan2, lCross[indbck[2]],nan2,lCross[indbck[3]]] @@ -1312,8 +1314,10 @@ def _DataCam12D_plot_spectral(lData, key=None, if c1 and 'LOS' in lData[0]._dgeom['lCam'][0].Id.Cls: lCross, lHor, llab = [], [], [] for cc in lData[0]._dgeom['lCam']: - lCross += cc._get_plotL(Lplot=Lplot, proj='cross', multi=True) - lHor += cc._get_plotL(Lplot=Lplot, proj='hor', multi=True) + lCross += cc._get_plotL(Lplot=Lplot, proj='cross', + return_pts=True, multi=True) + lHor += cc._get_plotL(Lplot=Lplot, proj='hor', + return_pts=True, multi=True) if bck and nD == 2: crossbck = [lCross[indbck[0]],nan2,lCross[indbck[1]],nan2, lCross[indbck[2]],nan2,lCross[indbck[3]]] @@ -2040,8 +2044,10 @@ def _DataCam12D_plot_combine(lData, key=None, nchMax=_nchMax, ntMax=_ntMax, llHor[ii] = [None for jj in range(0,len(lData[ii]._dgeom['lCam']))] for jj in range(0,len(lData[ii]._dgeom['lCam'])): cc = lData[ii]._dgeom['lCam'][jj] - llCross[ii][jj] = cc._get_plotL(Lplot=Lplot, proj='cross', multi=True) - llHor[ii][jj] = cc._get_plotL(Lplot=Lplot, proj='hor', multi=True) + llCross[ii][jj] = cc._get_plotL(Lplot=Lplot, proj='cross', + return_pts=True, multi=True) + llHor[ii][jj] = cc._get_plotL(Lplot=Lplot, proj='hor', + return_pts=True, multi=True) if c2 and lis2D[ii] and bck: indbck = np.r_[lindr[ii][0,0], lindr[ii][0,-1], @@ -2635,8 +2641,10 @@ def _Data1D_plot_spectrogram(Data, tf, f, lpsd, lang, if c1 and 'LOS' in Data._dgeom['lCam'][0].Id.Cls: lCross, lHor, llab = [], [], [] for cc in Data._dgeom['lCam']: - lCross += cc._get_plotL(Lplot=Lplot, proj='cross', multi=True) - lHor += cc._get_plotL(Lplot=Lplot, proj='hor', multi=True) + lCross += cc._get_plotL(Lplot=Lplot, proj='cross', + return_pts=True, multi=True) + lHor += cc._get_plotL(Lplot=Lplot, proj='hor', + return_pts=True, multi=True) if bck and cc._is2D(): crossbck = [lCross[indbck[0]],nan2,lCross[indbck[1]],nan2, lCross[indbck[2]],nan2,lCross[indbck[3]]] @@ -3291,8 +3299,10 @@ def _Data_plot_svd(Data, chronos, s, topos, modes=None, if c1 and 'LOS' in Data._dgeom['lCam'][0].Id.Cls: lCross, lHor, llab = [], [], [] for cc in Data._dgeom['lCam']: - lCross += cc._get_plotL(Lplot=Lplot, proj='cross', multi=True) - lHor += cc._get_plotL(Lplot=Lplot, proj='hor', multi=True) + lCross += cc._get_plotL(Lplot=Lplot, proj='cross', + return_pts=True, multi=True) + lHor += cc._get_plotL(Lplot=Lplot, proj='hor', + return_pts=True, multi=True) if bck and cc._is2D(): crossbck = [lCross[indbck[0]],nan2,lCross[indbck[1]],nan2, lCross[indbck[2]],nan2,lCross[indbck[3]]] diff --git a/tofu/imas2tofu/_core.py b/tofu/imas2tofu/_core.py index 42ab6dce1..87f83258a 100644 --- a/tofu/imas2tofu/_core.py +++ b/tofu/imas2tofu/_core.py @@ -1516,6 +1516,11 @@ def get_summary(self, sep=' ', line='-', just='l', raise Exception("Valid return_ values are: %s"%str(lok)) return out + def __repr__(self): + if hasattr(self, 'get_summary'): + return self.get_summary(return_='msg', verb=False) + else: + return object.__repr__(self) #--------------------- # Methods for returning data diff --git a/tofu/utils.py b/tofu/utils.py index 067ad97ab..920c6fa78 100644 --- a/tofu/utils.py +++ b/tofu/utils.py @@ -1417,19 +1417,23 @@ def _get_summary(cls, lar, lcol, elif return_ == 'array': out = lar elif return_ == 'msg': - out = lmsg + out = table_sep.join(lmsg) else: lok = [False, True, 'msg', 'array'] raise Exception("Valid return_ values are: %s"%str(lok)) return out - def get_summary(self, sep=' ', line='-', just='l', - table_sep=None, verb=True, return_=False): - """ Summary description of the object content as a np.array of str """ - pass - + # def get_summary(self, sep=' ', line='-', just='l', + # table_sep=None, verb=True, return_=False): + # """ Summary description of the object content as a np.array of str """ + # pass + def __repr__(self): + if hasattr(self, 'get_summary'): + return self.get_summary(return_='msg', verb=False) + else: + return object.__repr__(self) ############################# diff --git a/tofu/version.py b/tofu/version.py index 45d35422e..27ca5b6a6 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.0-620-g4f1df50' \ No newline at end of file +__version__='1.4.0-621-gfcae638' \ No newline at end of file From 611bf7a33b9a5c0a486e10186944c67a9bb47fc6 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 28 Aug 2019 15:06:41 +0200 Subject: [PATCH 16/19] [Issue 158] __repr__ fixed for MultiIDSLoader --- tofu/imas2tofu/_core.py | 8 +++++--- tofu/utils.py | 17 ++++++----------- tofu/version.py | 2 +- 3 files changed, 12 insertions(+), 15 deletions(-) diff --git a/tofu/imas2tofu/_core.py b/tofu/imas2tofu/_core.py index 87f83258a..d2aeff233 100644 --- a/tofu/imas2tofu/_core.py +++ b/tofu/imas2tofu/_core.py @@ -1455,8 +1455,10 @@ def get_ids(self, ids=None, occ=None): #--------------------- def get_summary(self, sep=' ', line='-', just='l', - verb=True, return_=False): + table_sep=None, verb=True, return_=False): """ Summary description of the object content as a np.array of str """ + if table_sep is None: + table_sep = '\n\n' # ----------------------- # idd @@ -1502,7 +1504,7 @@ def get_summary(self, sep=' ', line='-', just='l', else: msg1 = '' if verb: - msg = '\n\n'.join([msg0,msg1]) + msg = table_sep.join([msg0,msg1]) print(msg) if return_ != False: if return_ == True: @@ -1510,7 +1512,7 @@ def get_summary(self, sep=' ', line='-', just='l', elif return_ == 'array': out = (a0, a1) elif return_ == 'msg': - out = (msg0, msg1) + out = table_sep.join([msg0,msg1]) else: lok = [False, True, 'msg', 'array'] raise Exception("Valid return_ values are: %s"%str(lok)) diff --git a/tofu/utils.py b/tofu/utils.py index 920c6fa78..782cd219f 100644 --- a/tofu/utils.py +++ b/tofu/utils.py @@ -1402,7 +1402,6 @@ def _get_summary(cls, lar, lcol, if table_sep is None: table_sep = '\n\n' - # Out if verb or return_ in [True,'msg']: lmsg = [cls._getcharray(ar, col, sep=sep, line=line, just=just) @@ -1417,23 +1416,19 @@ def _get_summary(cls, lar, lcol, elif return_ == 'array': out = lar elif return_ == 'msg': - out = table_sep.join(lmsg) + out = lmsg else: lok = [False, True, 'msg', 'array'] raise Exception("Valid return_ values are: %s"%str(lok)) return out - # def get_summary(self, sep=' ', line='-', just='l', - # table_sep=None, verb=True, return_=False): - # """ Summary description of the object content as a np.array of str """ - # pass + def get_summary(self, sep=' ', line='-', just='l', + table_sep=None, verb=True, return_=False): + """ Summary description of the object content as a np.array of str """ + pass + - def __repr__(self): - if hasattr(self, 'get_summary'): - return self.get_summary(return_='msg', verb=False) - else: - return object.__repr__(self) ############################# diff --git a/tofu/version.py b/tofu/version.py index 27ca5b6a6..323951183 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.0-621-gfcae638' \ No newline at end of file +__version__='1.4.0-622-g89af109' \ No newline at end of file From 88b7c2d88abefe90b2d50717375acb73a5c45763 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 28 Aug 2019 18:54:30 +0200 Subject: [PATCH 17/19] [Issue 158] indout now refers to config.lStruct instead of lStruct_compute_kInOut, debugged and tested, and reflection type is inferred from Struct in config --- tofu/geom/_comp.py | 1 + tofu/geom/_core.py | 233 ++++++++++++++++++++++++--------------------- tofu/utils.py | 16 +++- tofu/version.py | 2 +- 4 files changed, 140 insertions(+), 112 deletions(-) diff --git a/tofu/geom/_comp.py b/tofu/geom/_comp.py index 1ee06cabe..84c3fa804 100644 --- a/tofu/geom/_comp.py +++ b/tofu/geom/_comp.py @@ -536,6 +536,7 @@ def LOS_CrossProj(VType, Ds, us, kOuts, proj='All', multi=False, # Use optimized get sample DL = np.vstack((np.zeros((nlos*nseg,),dtype=float), kOuts.ravel())) + k, reseff, lind = _GG.LOS_get_sample(nlos*nseg, resnk, DL, dmethod='rel', method='simps', num_threads=num_threads, Test=Test) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index cfe759094..bcecc2ba7 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -416,10 +416,10 @@ def _checkformat_inputs_dreflect(self, Types=None, coefs=None): assert type(Types) in [str, np.ndarray] if type(Types) is str: assert Types in self._DREFLECT_DTYPES.keys() - Types = np.full((self.nseg,), self._DREFLECT_DTYPES[Types], dtype=int) + Types = np.full((self.nseg+2,), self._DREFLECT_DTYPES[Types], dtype=int) else: Types = Types.astype(int).ravel() - assert Types.shape == () + assert Types.shape == (self.nseg+2,) Typesu = np.unique(Types) lc = np.array([Typesu == vv for vv in self._DREFLECT_DTYPES.values()]) @@ -1013,21 +1013,8 @@ def _get_phithetaproj_dist(self, refpt=None, ntheta=None, nphi=None, return dist, nDphi, Dphi, nDtheta, Dtheta - - def get_reflections(self, u, vperp, indout2=None, lamb=None): - """ Return the reflected unit vectors from input unit vectors and vperp - - The reflected unit vector depends on the incoming LOS (u), - the local normal unit vector (vperp), and the polygon segment hit - (indout2) - Future releases: dependence on lambda - - Also return per-LOS reflection Types (0:specular, 1:diffusive, 2:ccube) - - """ - - # Get per-LOS reflection Types and associated indices - Types = self._dreflect['Types'][indout2] + @staticmethod + def _get_reflections_ufromTypes(u, vperp, Types): indspec = Types == 0 inddiff = Types == 1 indcorn = Types == 2 @@ -1059,7 +1046,27 @@ def get_reflections(self, u, vperp, indout2=None, lamb=None): if np.any(indcorn): u2[:,indcorn] = -u[:,indcorn] - return u2, Types + return u2 + + def get_reflections(self, indout2, u=None, vperp=None): + """ Return the reflected unit vectors from input unit vectors and vperp + + The reflected unit vector depends on the incoming LOS (u), + the local normal unit vector (vperp), and the polygon segment hit + (indout2) + Future releases: dependence on lambda + + Also return per-LOS reflection Types (0:specular, 1:diffusive, 2:ccube) + + """ + + # Get per-LOS reflection Types and associated indices + Types = self._dreflect['Types'][indout2] + u2 = None + if u is not None: + assert vperp is not None + u2 = self._get_reflections_ufromTypes(u, vperp, Types) + return Types, u2 @@ -2500,6 +2507,34 @@ def get_summary(self, sep=' ', line='-', just='l', sep=sep, line=line, table_sep=table_sep, verb=verb, return_=return_) + def get_reflections(self, indout, u=None, vperp=None): + + # Get global Types array + lS = self.lStruct + + # Version only usable when indout returns npts+1 and npts+2 instead of + # -1 and -2 + # ls = [ss._dreflect['Types'].size for ss in lS] + # Types = np.empty((len(lS), np.max(ls)), dtype=int) + # for ii,ss in enumerate(lS): + # Types[ii,:ls[ii]] = ss._dreflect['Types'] + # # Deduce Types + # Types = Types[indout[0,:], indout[2,:]] + + iu = np.unique(indout[0,:]) + Types = np.empty((indout.shape[1],), dtype=int) + for ii in iu: + ind = indout[0,:] == ii + Types[ind] = lS[ii]._dreflect['Types'][indout[2,ind]] + + # Deduce u2 + u2 = None + if u is not None: + assert vperp is not None + u2 = Struct._get_reflections_ufromTypes(u, vperp, Types) + return Types, u2 + + def _get_phithetaproj_dist(self, refpt=None, ntheta=None, nphi=None, theta=None, phi=None): # Prepare repf @@ -2675,11 +2710,11 @@ def _reflect_Types(self, indout=None, Type=None, nRays=None): assert Type in ['specular', 'diffusive', 'ccube'] Types = np.full((nRays,), _DREFLECT[Type], dtype=int) else: - Types = None + Types = self.get_reflections(indout)[0] return Types - def _reflect_geom(self, u, vperp, indout=None, Type=None): + def _reflect_geom(self, u=None, vperp=None, indout=None, Type=None): assert u.shape == vperp.shape and u.shape[0] == 3 if indout is not None: assert indout.shape == (3,u.shape[1]) @@ -2687,37 +2722,8 @@ def _reflect_geom(self, u, vperp, indout=None, Type=None): # Get Types of relection for each Ray Types = self._reflect_Types(indout=indout, Type=Type, nRays=u.shape[1]) - # Get indices of each Type - indspec = Types == 0 - inddiff = Types == 1 - indcorn = Types == 2 - - u2 = np.full(u.shape, np.nan) - if np.any(np.logical_or(indspec,inddiff)): - vpar = np.array([vperp[1,:]*u[2,:] - vperp[2,:]*u[1,:], - vperp[2,:]*u[0,:] - vperp[0,:]*u[2,:], - vperp[0,:]*u[1,:] - vperp[1,:]*u[0,:]]) - vpar = np.array([vpar[1,:]*vperp[2,:] - vpar[2,:]*vperp[1,:], - vpar[2,:]*vperp[0,:] - vpar[0,:]*vperp[2,:], - vpar[0,:]*vperp[1,:] - vpar[1,:]*vperp[0,:]]) - vpar = vpar / np.sqrt(np.sum(vpar**2, axis=0))[None,:] - - if np.any(indspec): - # Compute u2 for specular - sca = np.sum(u[:,indspec]*vperp[:,indspec],axis=0,keepdims=True) - sca2 = np.sum(u[:,indspec]*vpar[:,indspec],axis=0,keepdims=True) - assert np.all(sca<=0.) and np.all(sca>=-1.) - assert np.all(sca2>=0.) and np.all(sca<=1.) - u2[:,indspec] = - sca*vperp[:,indspec] + sca2*vpar[:,indspec] - - if np.any(inddiff): - # Compute u2 for diffusive - sca = 2.*(np.random.random((1,inddiff.sum()))-0.5) - u2[:,inddiff] = (np.sqrt(1.-sca**2) * vperp[:,inddiff] - + sca * vpar[:,inddiff]) - - if np.any(indcorn): - u2[:,indcorn] = -u[:,indcorn] + # Deduce u2 + u2 = Struct._get_reflections_ufromTypes(u, vperp, Types) return u2, Types @@ -3299,7 +3305,7 @@ def _checkformat_inputs_dmisc(cls, color=None): @staticmethod def _get_keys_dgeom(): lk = ['D','u','pinhole', 'nRays', - 'kIn', 'kOut', 'PkIn', 'PkOut', 'vperp', 'indout', + 'kIn', 'kOut', 'PkIn', 'PkOut', 'vperp', 'indout', 'indStruct', 'kRMin', 'PRMin', 'RMin', 'isImage', 'Etendues', 'Surfaces', 'dX12', 'dreflect'] return lk @@ -3574,9 +3580,9 @@ def _complete_dX12(self, dgeom): - def _prepare_inputs_kInOut(self, D=None, u=None): + def _prepare_inputs_kInOut(self, D=None, u=None, indStruct=None): - # Prepare input + # Prepare input: D, u if D is None: D = np.ascontiguousarray(self.D) else: @@ -3587,8 +3593,10 @@ def _prepare_inputs_kInOut(self, D=None, u=None): u = np.ascontiguousarray(u) assert D.shape == u.shape - # Get reference - lS = self.lStruct_computeInOut + # Get reference: lS + if indStruct is None: + indStruct = self.indStruct_computeInOut + lS = [ss for ii,ss in enumerate(self.config.lStruct) if ii in indStruct] lSIn = [ss for ss in lS if ss._InOut=='in'] if len(lSIn)==0: @@ -3677,12 +3685,17 @@ def _prepare_inputs_kInOut(self, D=None, u=None): rmin=-1, forbid=True, eps_uz=1.e-6, eps_vz=1.e-9, eps_a=1.e-9, eps_b=1.e-9, eps_plane=1.e-9, test=True) - return largs, dkwd + return indStruct, largs, dkwd - def _compute_kInOut(self): + def _compute_kInOut(self, largs=None, dkwd=None, indStruct=None): # Prepare inputs - largs, dkwd = self._prepare_inputs_kInOut() + if largs is None: + indStruct, largs, dkwd =\ + self._prepare_inputs_kInOut(indStruct=indStruct) + else: + assert dkwd is not None + assert indStruct is not None if self._method=='ref': # call the dedicated function @@ -3696,7 +3709,10 @@ def _compute_kInOut(self): kIn, kOut, vperp, indout = out else: pass - return kIn, kOut, vperp, indout + + # Make sure indices refer to lStruct + indout[0,:] = indStruct[indout[0,:]] + return kIn, kOut, vperp, indout, indStruct def compute_dgeom(self, extra=True, plotdebug=True): @@ -3711,7 +3727,7 @@ def compute_dgeom(self, extra=True, plotdebug=True): self._dgeom = self._complete_dX12(self._dgeom) # Perform computation of kIn and kOut - kIn, kOut, vperp, indout = self._compute_kInOut() + kIn, kOut, vperp, indout, indStruct = self._compute_kInOut() # Clean up (in case of nans) ind = np.isnan(kIn) @@ -3746,7 +3762,8 @@ def compute_dgeom(self, extra=True, plotdebug=True): kIn[ind] = 0. # Update dgeom - dd = {'kIn':kIn, 'kOut':kOut, 'vperp':vperp, 'indout':indout} + dd = {'kIn':kIn, 'kOut':kOut, 'vperp':vperp, + 'indout':indout, 'indStruct':indStruct} self._dgeom.update(dd) # Run extra computations @@ -3918,9 +3935,6 @@ def get_reflections_as_cam(self, Type=None, Name=None, nb=None): nb = 1 nb = int(nb) assert nb > 0 - Ds = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u - us, Types = self.config._reflect_geom(self.u, self._dgeom['vperp'], - Type=Type) if Name is None: Name = self.Id.Name + '_Reflect%s'%str(Type) clas = Rays if self.__class__.__name__ == Rays else CamLOS1D @@ -3928,8 +3942,9 @@ def get_reflections_as_cam(self, Type=None, Name=None, nb=None): # Run first iteration Types = np.full((nb,self.nRays), 0, dtype=int) Ds = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u - us, Types[0,:] = self.config._reflect_geom(self.u, - self._dgeom['vperp'], + us, Types[0,:] = self.config._reflect_geom(u=self.u, + vperp=self._dgeom['vperp'], + indout=self._dgeom['indout'], Type=Type) lcam = [clas(dgeom=(Ds,us), config=self.config, Exp=self.Id.Exp, Diag=self.Id.Diag, @@ -3937,17 +3952,19 @@ def get_reflections_as_cam(self, Type=None, Name=None, nb=None): if nb == 1: return lcam[0], Types[0,:] - largs, dkwd = self._prepare_inputs_kInOut(D=Ds, u=us) - kouts, vperps, indouts = _GG.LOS_Calc_PInOut_VesStruct(*largs, **dkwd)[1:] + indStruct, largs, dkwd = self._prepare_inputs_kInOut(D=Ds, u=us, + indStruct=self._dgeom['indStruct']) + outi = self._compute_kInOut(largs=largs, dkwd=dkwd, indStruct=indStruct) + kouts, vperps, indouts = outi[1:-1] # Run other iterations for ii in range(1,nb): Ds = Ds + (kouts[None,:]-1.e-12) * us - us, Types[ii,:] = self.config._reflect_geom(us, vperps, - Type=Type) - kouts, vperps, indouts = _GG.LOS_Calc_PInOut_VesStruct(Ds, us, - *largs[2:], - **dkwd)[1:] + us, Types[ii,:] = self.config._reflect_geom(u=us, vperp=vperps, + indout=indouts, Type=Type) + outi = self._compute_kInOut(largs=[Dsi,usi,largs[2],largs[3]], + dkwd=dkwd, indStruct=indStruct) + kouts, vperps, indouts = outi[1:-1] lcam.append(clas(dgeom=(Ds,us), config=self.config, Exp=self.Id.Exp, Diag=self.Id.Diag, Name=Name, shot=self.Id.shot)) @@ -3984,22 +4001,25 @@ def add_reflections(self, Type=None, nb=None): # Run first iteration Ds[:,:,0] = self.D + (self._dgeom['kOut'][None,:]-1.e-12) * self.u - us[:,:,0], Types[:,0] = self.config._reflect_geom(self.u, - self._dgeom['vperp'], + us[:,:,0], Types[:,0] = self.config._reflect_geom(u=self.u, + vperp=self._dgeom['vperp'], + indout=self._dgeom['indout'], Type=Type) - largs, dkwd = self._prepare_inputs_kInOut(D=Ds[:,:,0], u=us[:,:,0]) - outi = _GG.LOS_Calc_PInOut_VesStruct(*largs, **dkwd)[1:] - kouts[:,0], vperps[:,:,0], indouts[:,:,0] = outi + indStruct, largs, dkwd = self._prepare_inputs_kInOut(D=Ds[:,:,0], u=us[:,:,0], + indStruct=self._dgeom['indStruct']) + outi = self._compute_kInOut(largs=largs, dkwd=dkwd, indStruct=indStruct) + kouts[:,0], vperps[:,:,0], indouts[:,:,0] = outi[1:-1] # Run other iterations for ii in range(1,nb): Dsi = Ds[:,:,ii-1] + (kouts[None,:,ii-1]-1.e-12) * us[:,:,ii-1] - usi, Types[:,ii] = self.config._reflect_geom(us[:,:,ii-1], - vperps[:,:,ii-1], + usi, Types[:,ii] = self.config._reflect_geom(u=us[:,:,ii-1], + vperp=vperps[:,:,ii-1], + indout=indouts[:,:,ii-1], Type=Type) - outi = _GG.LOS_Calc_PInOut_VesStruct(Dsi, usi, - *largs[2:], **dkwd)[1:] - kouts[:,ii], vperps[:,:,ii], indouts[:,:,ii] = outi + outi = self._compute_kInOut(largs=[Dsi,usi,largs[2],largs[3]], + dkwd=dkwd, indStruct=indStruct) + kouts[:,ii], vperps[:,:,ii], indouts[:,:,ii] = outi[1:-1] Ds[:,:,ii], us[:,:,ii] = Dsi, usi self._dgeom['dreflect'] = {'nb':nb, 'Type':Type, 'Types':Types, @@ -4028,12 +4048,12 @@ def _strip_dgeom(self, strip=0): # strip if strip==1: lkeep = ['D','u','pinhole','nRays', - 'kIn','kOut','vperp','indout', 'kRMin', + 'kIn','kOut','vperp','indout', 'indStruct', 'kRMin', 'Etendues','Surfaces','isImage','dX12', 'dreflect'] utils.ToFuObject._strip_dict(self._dgeom, lkeep=lkeep) elif self._dstrip['strip']<=1 and strip>=2: lkeep = ['D','u','pinhole','nRays', - 'kIn','kOut','vperp','indout', + 'kIn','kOut','vperp','indout','indStruct', 'Etendues','Surfaces','isImage','dX12', 'dreflect'] utils.ToFuObject._strip_dict(self._dgeom, lkeep=lkeep) @@ -4214,17 +4234,17 @@ def config(self): return self._dconfig['Config'] @property - def lStruct_computeInOut(self): + def indStruct_computeInOut(self): compute = self.config.get_compute() lS = self.config.lStruct - lSI, lSO = [], [] - for ii in range(0,self._dconfig['Config']._dStruct['nObj']): + iI, iO = [], [] + for ii in range(0,len(lS)): if compute[ii]: - if lS[ii]._InOut=='in': - lSI.append(lS[ii]) - elif lS[ii]._InOut=='out': - lSO.append(lS[ii]) - return lSI+lSO + if lS[ii]._InOut == 'in': + iI.append(ii) + elif lS[ii]._InOut == 'out': + iO.append(ii) + return np.r_[iI + iO] @property def Etendues(self): @@ -4331,7 +4351,7 @@ def select(self, key=None, val=None, touch=None, log='any', out=int): Tuple indicating you want the rays that are touching some specific elements of self.config: - touch[0] : str / int or list of such str : a 'Cls_Name' string indicating the element - int : the index of the element in self.lStruct_computeInOut + int : the index of the element in self.config.lStruct - touch[1] : int / list of int Indices of the desired segments on the polygon (i.e.: of the cross-section polygon of the above element) @@ -4406,10 +4426,10 @@ def _check_touch(tt): raise Exception(msg) if cS: - lS = self.lStruct_computeInOut k0, k1 = touch[ii].split('_') - ind = [jj for jj in range(0,len(lS)) - if lS[jj].Id.Cls==k0 and lS[jj].Id.Name==k1] + lS = self.config.lStruct + ind = [jj for jj,ss in enumerate(lS) + if ss.Id.Cls == k0 and ss.Id.Name == k1] assert len(ind)==1 touch[ii] = [ind[0]] elif c0: @@ -4727,8 +4747,8 @@ def calc_kInkOut_IsoFlux(self, lPoly, lVIn=None, Lim=None, 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] = kin, kout + # PIn, POut, kin, kout, VperpIn, vperp, IIn, indout = out[] + 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]], @@ -5234,9 +5254,10 @@ def get_touch_dict(self, ind=None, out=bool): raise Exception(msg) dElt = {} + lS = self.config.lStruct ind = self._check_indch(ind, out=bool) - for ss in self.lStruct_computeInOut: - kn = "%s_%s"%(ss.__class__.__name__, ss.Id.Name) + for ii in self.indStruct_computeInOut: + kn = "%s_%s"%(lS[ii].__class__.__name__, lS[ii].Id.Name) indtouch = self.select(touch=kn, out=bool) if np.any(indtouch): indok = indtouch & ind @@ -5246,7 +5267,7 @@ def get_touch_dict(self, ind=None, out=bool): indok = indok.nonzero()[0] indout = indout.nonzero()[0] dElt[kn] = {'indok':indok, 'indout':indout, - 'col':ss.get_color()} + 'col':lS[ii].get_color()} return dElt def get_touch_colors(self, ind=None, dElt=None, @@ -5334,9 +5355,9 @@ def get_summary(self, sep=' ', line='-', just='l', ar1.append( [str(vv) for vv in v] ) # call base method - self._get_summary([ar0, ar1], [col0, col1], - sep=sep, line=line, table_sep=table_sep, - verb=verb, return_=return_) + return self._get_summary([ar0, ar1], [col0, col1], + sep=sep, line=line, table_sep=table_sep, + verb=verb, return_=return_) def __add__(self, other): if not other.__class__.__name__ == self.__class__.__name__: diff --git a/tofu/utils.py b/tofu/utils.py index 782cd219f..7679c6864 100644 --- a/tofu/utils.py +++ b/tofu/utils.py @@ -1402,6 +1402,7 @@ def _get_summary(cls, lar, lcol, if table_sep is None: table_sep = '\n\n' + # Out if verb or return_ in [True,'msg']: lmsg = [cls._getcharray(ar, col, sep=sep, line=line, just=just) @@ -1416,18 +1417,23 @@ def _get_summary(cls, lar, lcol, elif return_ == 'array': out = lar elif return_ == 'msg': - out = lmsg + out = table_sep.join(lmsg) else: lok = [False, True, 'msg', 'array'] raise Exception("Valid return_ values are: %s"%str(lok)) return out - def get_summary(self, sep=' ', line='-', just='l', - table_sep=None, verb=True, return_=False): - """ Summary description of the object content as a np.array of str """ - pass + # def get_summary(self, sep=' ', line='-', just='l', + # table_sep=None, verb=True, return_=False): + # """ Summary description of the object content as a np.array of str """ + # pass + def __repr__(self): + if hasattr(self, 'get_summary'): + return self.get_summary(return_='msg', verb=False) + else: + return object.__repr__(self) diff --git a/tofu/version.py b/tofu/version.py index 323951183..bd82af6cc 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.0-622-g89af109' \ No newline at end of file +__version__='1.4.0-623-g611bf7a' \ No newline at end of file From 9358d1126bdd8d5b6842e505f61effef1d2817ee Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 28 Aug 2019 19:26:43 +0200 Subject: [PATCH 18/19] [Issue 158] Added basics for reflections for calc_signal() and calc_signal_from_Plasma2D() --- tofu/geom/_core.py | 26 ++++++++++++++++++++++++++ tofu/version.py | 2 +- 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index bcecc2ba7..bfd2d4c3e 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -4871,6 +4871,7 @@ def _calc_signal_postformat(self, sig, Brightness=True, dataname=None, t=None, def calc_signal(self, func, t=None, ani=None, fkwdargs={}, Brightness=True, res=None, DL=None, resMode='abs', method='sum', minimize='calls', num_threads=16, + reflections=True, coefs=None, ind=None, out=object, plot=True, dataname=None, fs=None, dmargin=None, wintit=None, invert=True, units=None, draw=True, connect=True, newcalc=True): @@ -4934,6 +4935,18 @@ def calc_signal(self, func, t=None, ani=None, fkwdargs={}, Brightness=True, dmethod=resMode, method=method, ani=ani, t=t, fkwdargs=fkwdargs, minimize=minimize, num_threads=num_threads, Test=True) + c0 = (reflections and self._dgeom['dreflect'] is not None + and self._dgeom['dreflect'].get('nb',0) > 0) + if c0: + if coefs is None: + coefs = 1. + for ii in range(self._dgeom['dreflect']['nb']): + Dsi = np.ascontiguousarray(self._dgeom['dreflect']['Ds'][:,:,ii]) + usi = np.ascontiguousarray(self._dgeom['dreflect']['us'][:,:,ii]) + s += coefs*_GG.LOS_calc_signal(func, Dsi, usi, res, DL, + dmethod=resMode, method=method, ani=ani, + t=t, fkwdargs=fkwdargs, minimize=minimize, + num_threads=num_threads, Test=True) # Integrate if s.ndim == 2: @@ -4987,6 +5000,7 @@ def calc_signal_from_Plasma2D(self, plasma2d, t=None, newcalc=False, res=None, DL=None, resMode='abs', method='sum', minimize='calls', num_threads=16, + reflections=True, coefs=None, ind=None, out=object, plot=True, dataname=None, fs=None, dmargin=None, wintit=None, invert=True, units=None, draw=True, connect=True): @@ -5036,6 +5050,18 @@ def calc_signal_from_Plasma2D(self, plasma2d, t=None, newcalc=False, dmethod=resMode, method=method, ani=ani, t=t, fkwdargs={}, minimize=minimize, Test=True, num_threads=num_threads) + c0 = (reflections and self._dgeom['dreflect'] is not None + and self._dgeom['dreflect'].get('nb',0) > 0) + if c0: + if coefs is None: + coefs = 1. + for ii in range(self._dgeom['dreflect']['nb']): + Dsi = np.ascontiguousarray(self._dgeom['dreflect']['Ds'][:,:,ii]) + usi = np.ascontiguousarray(self._dgeom['dreflect']['us'][:,:,ii]) + sig += coefs*_GG.LOS_calc_signal(funcbis, Dsi, usi, res, DL, + dmethod=resMode, method=method, ani=ani, + t=t, fkwdargs=fkwdargs, minimize=minimize, + num_threads=num_threads, Test=True) else: # Get ptsRZ along LOS // Which to choose ??? pts, reseff, indpts = self.get_sample(res, resMode=resMode, DL=DL, method=method, ind=ind, diff --git a/tofu/version.py b/tofu/version.py index bd82af6cc..aa4f615d4 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.0-623-g611bf7a' \ No newline at end of file +__version__='1.4.0-624-g88b7c2d' \ No newline at end of file From c2c38d86882a70e530f8a1550bee40903109d380 Mon Sep 17 00:00:00 2001 From: VEZINET Didier Date: Wed, 28 Aug 2019 19:34:42 +0200 Subject: [PATCH 19/19] [Issue 158] Debugged unit tests for Data.select_ch() (due to new indout indexing) --- tofu/tests/tests02_data/tests03_core.py | 7 ++++--- tofu/version.py | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/tofu/tests/tests02_data/tests03_core.py b/tofu/tests/tests02_data/tests03_core.py index 0a09faef2..2dbe55793 100644 --- a/tofu/tests/tests02_data/tests03_core.py +++ b/tofu/tests/tests02_data/tests03_core.py @@ -254,12 +254,13 @@ def test03_select_t(self): def test04_select_ch(self): for oo in self.lobj: if oo.dgeom['lCam'] is not None: - name = [k for k in oo.config.dStruct['lorder'] + name = [(ii,k) for ii,k in + enumerate(oo.config.dStruct['lorder']) if 'Ves' in k or 'PlasmaDomain' in k] assert len(name) == 1 - ind = oo.select_ch(touch=name[0], out=bool) + ind = oo.select_ch(touch=name[0][1], out=bool) assert ind.sum() > 0, (ind.sum(), ind) - assert np.allclose(ind, oo.select_ch(touch=0, + assert np.allclose(ind, oo.select_ch(touch=name[0][0], out=bool)) if len(oo.dchans().keys()) > 0: ind = oo.select_ch(key='Name', val=['c0','c10'], diff --git a/tofu/version.py b/tofu/version.py index aa4f615d4..ab9e7aaa3 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.0-624-g88b7c2d' \ No newline at end of file +__version__='1.4.0-625-g9358d11' \ No newline at end of file