Feature #252: clearer error and check for emissivity function

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)

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

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,

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