diff --git a/examples/tutorials/tuto_plot_custom_emissivity.py b/examples/tutorials/tuto_plot_custom_emissivity.py index dc8676037..9d05cdb2c 100644 --- a/examples/tutorials/tuto_plot_custom_emissivity.py +++ b/examples/tutorials/tuto_plot_custom_emissivity.py @@ -31,13 +31,15 @@ # Now, we define an emissivity function that depends on r and z coordinates. # We can plot its profile in the (0, X, Z) plane. + def emissivity(pts, t=None, vect=None): """Custom emissivity as a function of geometry. :param pts: ndarray of shape (3, n_points) (each column is a xyz coordinate) :param t: optional, time parameter to add a time dependency to the emissivity function - :param vect: + :param vect: optional, ndarray of shape (3, n_points), if anisotropic + emissivity, unit direction vectors (X,Y,Z) :return: - emissivity -- 2D array holding the emissivity for each point in the input grid @@ -88,6 +90,7 @@ def project_to_2D(xyz): method="sum", newcalc=True, plot=False, + ani=False, t=time_vector) sig.plot(ntMax=1) diff --git a/tofu/geom/_GG.pyx b/tofu/geom/_GG.pyx index c645e7e5d..4ffa3a155 100644 --- a/tofu/geom/_GG.pyx +++ b/tofu/geom/_GG.pyx @@ -9,10 +9,7 @@ from warnings import warn import numpy as np import scipy.integrate as scpintg from matplotlib.path import Path -if sys.version[0]=='3': - from inspect import signature as insp -elif sys.version[0]=='2': - from inspect import getargspec as insp + # -- Cython libraries imports -------------------------------------------------- from cpython cimport bool from cpython.array cimport array, clone @@ -851,7 +848,7 @@ def discretize_vpoly(double[:,::1] VPoly, double dL, reference Radius coordinates, not shifted even if DIn <> 0. If DIn == 0, then Rref_arr = PtsCross[0, ...] VPolybis : - + """ cdef int NP=VPoly.shape[1] cdef int[1] sz_vb, sz_ot @@ -2771,74 +2768,6 @@ def LOS_get_sample(int nlos, dL, double[:,::1] los_lims, str dmethod='abs', ###################################################################### # Signal calculation ###################################################################### -cdef get_insp(ff): - out = insp(ff) - if sys.version[0]=='3': - pars = out.parameters.values() - na = np.sum([(pp.kind==pp.POSITIONAL_OR_KEYWORD - and pp.default is pp.empty) for pp in pars]) - kw = [pp.name for pp in pars if (pp.kind==pp.POSITIONAL_OR_KEYWORD - and pp.default is not pp.empty)] - else: - nat, nak = len(out.args), len(out.defaults) - na = nat-nak - kw = [out.args[ii] for ii in range(nat-1,na-1,-1)][::-1] - return na, kw - - - -def check_ff(ff, t=None, Ani=None, bool Vuniq=False): - cdef bool ani - cdef int nt = -1 - # ... - stre = "Input emissivity function (ff)" - assert hasattr(ff,'__call__'), stre+" must be a callable (function)!" - na, kw = get_insp(ff) - assert na==1, stre+" must take only one positional argument: ff(Pts)!" - assert 't' in kw, stre+" must have kwarg 't=None' for time vector!" - C = type(t) in [int,float,np.int64,np.float64] or hasattr(t,'__iter__') - assert t is None or C, "Arg t must be None, a scalar or an iterable!" - Pts = np.array([[1,2],[3,4],[5,6]]) - NP = Pts.shape[1] - try: - out = ff(Pts, t=t) - except Exception: - Str = stre+" must take one positional arg: a (3,N) np.ndarray" - assert False, Str - if hasattr(t,'__iter__'): - nt = len(t) - Str = ("ff(Pts,t=t), where Pts is a (3,N) np.array and " - +"t a len()=nt iterable, must return a (nt,N) np.ndarray!") - assert type(out) is np.ndarray and out.shape==(nt,NP), Str - else: - Str = ("When fed a (3,N) np.array only, or if t is a scalar," - +" ff must return a (N,) np.ndarray!") - assert type(out) is np.ndarray and out.shape==(NP,), Str - - ani = ('Vect' in kw) if Ani is None else Ani - if ani: - Str = "If Ani=True, ff must take a keyword argument 'Vect=None'!" - assert 'Vect' in kw, Str - Vect = np.array([1,2,3]) if Vuniq else np.ones(Pts.shape) - try: - out = ff(Pts, Vect=Vect, t=t) - except Exception: - Str = "If Ani=True, ff must handle multiple points Pts (3,N) with " - if Vuniq: - Str += "a unique common vector (Vect as a len()=3 iterable)" - else: - Str += "multiple vectors (Vect as a (3,N) np.ndarray)" - assert False, Str - if hasattr(t,'__iter__'): - Str = ("If Ani=True, ff must return a (nt,N) np.ndarray when " - +"Pts is (3,N), Vect is provided and t is (nt,)") - assert type(out) is np.ndarray and out.shape==(nt,NP), Str - else: - Str = ("If Ani=True, ff must return a (nt,N) np.ndarray when " - +"Pts is (3,N), Vect is provided and t is (nt,)") - assert type(out) is np.ndarray and out.shape==(NP,), Str - return ani - def integrate1d(y, double dx, t=None, str method='sum'): """ Generic integration method ['sum','simps','romb'] @@ -2931,7 +2860,6 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, cdef double[1] loc_eff_res cdef double[::1] reseff_mv cdef double[::1] res_mv - cdef double[::1,:] sig_mv cdef double[:,::1] val_mv cdef double[:,::1] pts_mv cdef double[:,::1] usbis_mv @@ -2999,7 +2927,6 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, n_imode = _st.get_nb_imode(imode) # Initialization result sig = np.empty((nt, nlos), dtype=float, order='F') - sig_mv = sig # If the resolution is the same for every LOS, we create a tab if res_is_list : res_arr = np.asarray(res) @@ -3096,11 +3023,12 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, pts, usbis = _st.call_get_sample_single_ani(lims[0, ii], lims[1, ii], res_mv[ii], - n_dmode, n_imode, + n_dmode, + n_imode, &loc_eff_res[0], &nb_rows[0], - ray_orig[:,ii:ii+1], - ray_vdir[:,ii:ii+1]) + ray_orig[:, ii:ii+1], + ray_vdir[:, ii:ii+1]) # loop over time for calling and integrating for jj in range(nt): val = func(pts, t=ltime[jj], vect=-usbis, **fkwdargs) @@ -3200,7 +3128,7 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, ray_orig[:, ii:ii+1], ray_vdir[:, ii:ii+1]) val_2d = func(pts, t=t, vect=-usbis, **fkwdargs) - sig_mv[:, ii] = np.sum(val_2d, axis=-1)*loc_eff_res[0] + sig[:, ii] = np.sum(val_2d, axis=-1)*loc_eff_res[0] elif n_imode == 1: # simpson integration mode for ii in range(nlos): pts, usbis = _st.call_get_sample_single_ani(lims[0, ii], @@ -3233,39 +3161,42 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, # -- not anisotropic ----------------------------------------------- if n_imode == 0: # "sum" integration mode for ii in range(nlos): - pts = _st.call_get_sample_single(lims[0,ii], lims[1,ii], + pts = _st.call_get_sample_single(lims[0, ii], + lims[1, ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], &nb_rows[0], - ray_orig[:,ii:ii+1], - ray_vdir[:,ii:ii+1]) + ray_orig[:, ii:ii+1], + ray_vdir[:, ii:ii+1]) val_2d = func(pts, t=t, **fkwdargs) - sig[:, ii] = np.sum(val_2d,axis=-1)*loc_eff_res[0] + sig[:, ii] = np.sum(val_2d, axis=-1)*loc_eff_res[0] elif n_imode == 1: # "simpson" integration mode for ii in range(nlos): - pts = _st.call_get_sample_single(lims[0,ii], lims[1,ii], + pts = _st.call_get_sample_single(lims[0, ii], + lims[1, ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], &nb_rows[0], - ray_orig[:,ii:ii+1], - ray_vdir[:,ii:ii+1]) + ray_orig[:, ii:ii+1], + ray_vdir[:, ii:ii+1]) val = func(pts, t=t, **fkwdargs) sig[:, ii] = scpintg.simps(val, x=None, axis=-1, dx=loc_eff_res[0]) elif n_imode == 2: # "romberg" integration mode for ii in range(nlos): - pts = _st.call_get_sample_single(lims[0,ii], lims[1,ii], + pts = _st.call_get_sample_single(lims[0, ii], + lims[1, ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], &nb_rows[0], - ray_orig[:,ii:ii+1], - ray_vdir[:,ii:ii+1]) + ray_orig[:, ii:ii+1], + ray_vdir[:, ii:ii+1]) val = func(pts, t=t, **fkwdargs) sig[:, ii] = scpintg.romb(val, show=False, axis=1, - dx=loc_eff_res[0]) + dx=loc_eff_res[0]) return sig @@ -3379,7 +3310,7 @@ cdef inline void NEW_LOS_sino_Tor(double orig0, double orig1, double orig2, kPMin = res[0] if is_LOS_Mode and kPMin > kOut: kPMin = kOut - + # Computing the point's coordinates......................................... cdef double PMin0 = orig0 + kPMin * dirv0 cdef double PMin1 = orig1 + kPMin * dirv1 diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index af9c03e15..1bd4a0f98 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -5760,6 +5760,100 @@ def calc_min_rho_from_Plasma2D(self, plasma, t=None, log='min', for vv in np.split(val, lind, axis=-1)]) return vals, pts, t + def get_inspector(self, ff): + out = inspect.signature(ff) + pars = out.parameters.values() + na = np.sum([(pp.kind == pp.POSITIONAL_OR_KEYWORD + and pp.default is pp.empty) for pp in pars]) + kw = [pp.name for pp in pars if (pp.kind == pp.POSITIONAL_OR_KEYWORD + and pp.default is not pp.empty)] + return na, kw + + def check_ff(self, ff, t=None, ani=None): + # Initialization of function wrapper + wrapped_ff = ff + + # Define unique error message giving all info in a concise way + # Optionnally add error-specific line afterwards + msg = ("User-defined emissivity function ff must:\n" + + "\t- be a callable (function)\n" + + "\t- take only one positional arg " + + "and at least one keyword arg:\n" + + "\t\t - ff(pts, t=None), where:\n" + + "\t\t\t - pts is a (3, npts) of (x, y, z) coordinates\n" + + "\t\t\t - t can be None / scalar / iterable of len(t) = nt\n" + + "\t- Always return a 2d (nt, npts) np.ndarray, where:\n" + + "\t\t - nt = len(t) if t is an iterable\n" + + "\t\t - nt = 1 if t is None or scalar\n" + + "\t\t - npts is the number of pts (pts.shape[1])\n\n" + + "\t- Optionally, ff can take an extra keyword arg:\n" + + "\t\t - ff(pts, vect=None, t=None), where:\n" + + "\t\t\t - vect is a (3, npts) np.ndarray\n" + + "\t\t\t - vect contains the (x, y, z) coordinates " + + "of the units vectors of the photon emission directions" + + "for each pts. Present only for anisotropic emissivity, " + + "unless specifically indicated otherwise " + + "(with ani=False in LOS_calc_signal).\n" + + "\t\t\tDoes not affect the outpout shape (still (nt, npts))") + + # .. Checking basic definition of function .......................... + if not hasattr(ff, '__call__'): + msg += "\n\n => ff must be a callable (function)!" + raise Exception(msg) + + npos_args, kw = self.get_inspector(ff) + if npos_args != 1: + msg += "\n\n => ff must take only 1 positional arg: ff(pts)!" + raise Exception(msg) + + if 't' not in kw: + msg += "\n\n => ff must have kwarg 't=None' for time vector!" + raise Exception(msg) + + # .. Checking time vector ......................................... + ltypeok = [int, float, np.int64, np.float64] + is_t_type_valid = (type(t) in ltypeok or hasattr(t, '__iter__')) + if not (t is None or is_t_type_valid): + msg += "\n\n => t must be None, scalar or iterable !" + raise Exception(msg) + nt = len(t) if hasattr(t, '__iter__') else 1 + + # .. Test anisotropic case ....................................... + if ani is None: + is_ani = ('vect' in kw) + else: + assert isinstance(ani, bool) + is_ani = ani + + # .. Testing outputs ............................................... + test_pts = np.array([[1, 2], [3, 4], [5, 6]]) + npts = test_pts.shape[1] + if is_ani: + vect = np.ones(test_pts.shape) + try: + out = ff(test_pts, vect=vect, t=t) + except Exception: + msg += "\n\n => ff must take ff(pts, vect=vect, t=t) !" + raise Exception(msg) + else: + try: + out = ff(test_pts, t=t) + except Exception: + msg += "\n\n => ff must take a ff(pts, t=t) !" + raise Exception(msg) + + if not (isinstance(out, np.ndarray) and (out.shape == (nt, npts) + or out.shape == (npts,))): + msg += "\n\n => wrong output (always 2d np.ndarray) !" + raise Exception(msg) + + if nt == 1 and out.shape == (npts,): + def wrapped_ff(*args, **kwargs): + res_ff = ff(*args, **kwargs) + return np.reshape(res_ff, (1, -1)) + + return is_ani, wrapped_ff + def _calc_signal_preformat(self, ind=None, DL=None, t=None, out=object, Brightness=True): msg = "Arg out must be in [object,np.ndarray]" @@ -5886,7 +5980,7 @@ def calc_signal( self, func, t=None, - ani=False, + ani=None, fkwdargs={}, Brightness=True, res=None, @@ -5922,11 +6016,6 @@ def calc_signal( => the method returns W/m2 (resp. W/m2/sr) The line is sampled using :meth:`~tofu.geom.LOS.get_sample`, - The integral can be computed using three different methods: - - 'sum': A numpy.sum() on the local values (x segments lengths) - - 'simps': using :meth:`scipy.integrate.simps` - - 'romb': using :meth:`scipy.integrate.romb` - Except func, arguments common to :meth:`~tofu.geom.LOS.get_sample` Parameters @@ -5941,6 +6030,15 @@ def calc_signal( - vect: None / (3,N) np.ndarray, unit direction vectors (X,Y,Z) Should return at least: - val : (N,) np.ndarray, local emissivity values + method : string, the integral can be computed using 3 different methods + - 'sum': A numpy.sum() on the local values (x segments) DEFAULT + - 'simps': using :meth:`scipy.integrate.simps` + - 'romb': using :meth:`scipy.integrate.romb` + minimize : string, method to minimize for computation optimization + - "calls": minimal number of calls to `func` (default) + - "memory": slowest method, to use only if "out of memory" error + - "hybrid": mix of before-mentioned methods. + Returns ------- @@ -5966,6 +6064,7 @@ def calc_signal( # Launch # NB : find a way to exclude cases with DL[0,:]>=DL[1,:] !! # Exclude Rays not seeing the plasma if newcalc: + ani, func = self.check_ff(func, t=t, ani=ani) s = _GG.LOS_calc_signal( func, Ds, diff --git a/tofu/tests/tests09_tutorials/tests01_runall.py b/tofu/tests/tests09_tutorials/tests01_runall.py index 969554fa6..e20341724 100644 --- a/tofu/tests/tests09_tutorials/tests01_runall.py +++ b/tofu/tests/tests09_tutorials/tests01_runall.py @@ -148,13 +148,22 @@ def _test_tuto(cls, tuto=None, pathtuto=_HERE, root=_TFROOT): src = os.path.join(pathtuto, tuto + '.py') target = os.path.join(root, tuto + '.py') shutil.copyfile(src, target) + error = None try: cmd = 'python ' + target - out = subprocess.run(cmd, shell=True, check=True, + out = subprocess.run(cmd, shell=True, check=False, stdout=subprocess.PIPE, stderr=subprocess.PIPE) + ls = out.stderr.decode().split('\n') + if any(['Error' in ss and not 'ModuleNotFoundError' in ss + for ss in ls]): + error = '\n'.join(ls) except Exception as err: - raise err + error = err + + if error is not None: + msg = str(error) + raise Exception(msg) plt.close('all') # Remove temporary files and saved files diff --git a/tofu/utils.py b/tofu/utils.py index f727a83ca..475b6a338 100644 --- a/tofu/utils.py +++ b/tofu/utils.py @@ -2484,7 +2484,8 @@ def func(li, val=val, n1=n1, n2=n2): else: assert type(val) is np.ndarray - # val = np.atleast_1d(val.squeeze()) + if val.ndim > len(lrids) and val.ndim > ninds: + val = np.atleast_1d(np.squeeze(val)) ndim = val.ndim c0 = ndim >= len(lrids) and len(lrids) >= ninds and ndim >= ninds if not c0: diff --git a/tofu/version.py b/tofu/version.py index ac22e7403..31088c3e2 100644 --- a/tofu/version.py +++ b/tofu/version.py @@ -1,2 +1,2 @@ # Do not edit, pipeline versioning governed by git tags! -__version__ = '1.4.1-232-g56df4eb' +__version__ = '1.4.1-264-g1f49781'