diff --git a/examples/tutorials/plot_basic_tutorial.py b/examples/tutorials/tuto_plot_basic.py similarity index 100% rename from examples/tutorials/plot_basic_tutorial.py rename to examples/tutorials/tuto_plot_basic.py diff --git a/examples/tutorials/plot_create_geometry.py b/examples/tutorials/tuto_plot_create_geometry.py similarity index 98% rename from examples/tutorials/plot_create_geometry.py rename to examples/tutorials/tuto_plot_create_geometry.py index 12306022b..dcd3a586e 100644 --- a/examples/tutorials/plot_create_geometry.py +++ b/examples/tutorials/tuto_plot_create_geometry.py @@ -9,6 +9,12 @@ import matplotlib.pyplot as plt import tofu.geom as tfg +import os +import tofu as tf +print(os.getcwd()) +print(tf.__version__) +print(tf.__path__) + ############################################################################### # Creating an empty Vessel @@ -86,7 +92,7 @@ # Plot the polygon by default in two projections (cross-section and horizontal) # and return the list of axes -Lax = ves.plot(element="P") +lax = ves.plot(element="P") ############################################################################### diff --git a/examples/tutorials/plot_custom_emissivity.py b/examples/tutorials/tuto_plot_custom_emissivity.py similarity index 96% rename from examples/tutorials/plot_custom_emissivity.py rename to examples/tutorials/tuto_plot_custom_emissivity.py index 349aa3d45..dc8676037 100644 --- a/examples/tutorials/plot_custom_emissivity.py +++ b/examples/tutorials/tuto_plot_custom_emissivity.py @@ -84,8 +84,11 @@ def project_to_2D(xyz): sig, units = cam2d.calc_signal(emissivity, res=0.01, reflections=False, + minimize="hybrid", + method="sum", newcalc=True, plot=False, t=time_vector) + sig.plot(ntMax=1) -plt.show() +plt.show(block=False) diff --git a/tofu/geom/_GG.pyx b/tofu/geom/_GG.pyx index 0462264b8..c645e7e5d 100644 --- a/tofu/geom/_GG.pyx +++ b/tofu/geom/_GG.pyx @@ -2920,12 +2920,12 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, cdef int sz1_us, sz2_us cdef int sz1_dls, sz2_dls cdef int n_imode, n_dmode + cdef int nlos + cdef int nt=0, ii, jj cdef bint res_is_list cdef bint C0, C1 cdef list ltime cdef double loc_r - cdef unsigned int nlos - cdef unsigned int nt=0, ii, jj cdef long[1] nb_rows cdef long[::1] indbis cdef double[1] loc_eff_res @@ -2946,7 +2946,7 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, cdef long* ind_arr = NULL cdef double* reseff_arr = NULL cdef double** coeff_ptr = NULL - # .. ray_orig shape needed for testing and in algo ............................... + # .. ray_orig shape needed for testing and in algo ......................... sz1_ds = ray_orig.shape[0] nlos = ray_orig.shape[1] res_is_list = hasattr(res, '__iter__') @@ -2959,7 +2959,8 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, assert sz1_ds == 3, "Dim 0 of arg ray_orig should be 3" assert sz1_us == 3, "Dim 0 of arg ray_vdir should be 3" assert sz1_dls == 2, "Dim 0 of arg lims should be 2" - error_message = "Args ray_orig, ray_vdir, lims should have same dimension 1" + error_message = ("Args ray_orig, ray_vdir, and lims " + + "should have same dimension 1") assert nlos == sz2_us == sz2_dls, error_message C0 = not res_is_list and res > 0. C1 = res_is_list and len(res)==nlos and np.all(res>0.) @@ -2997,7 +2998,7 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, n_dmode = _st.get_nb_dmode(dmode) n_imode = _st.get_nb_imode(imode) # Initialization result - sig = np.empty((nt,nlos),dtype=float,order='F') + 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 : @@ -3008,7 +3009,8 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, # -------------------------------------------------------------------------- # Minimize function calls: sample (vect), call (once) and integrate if minim == 'calls': - if not method=='sum': + if n_imode != 0: + # Integration mode is Simpson or Romberg # Discretize all LOS k, reseff, ind = LOS_get_sample(nlos, res_arr, lims, dmethod=dmode, method=imode, @@ -3016,7 +3018,7 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, nbrep = np.r_[ind[0], np.diff(ind), k.size - ind[nlos-2]] # get pts and values usbis = np.repeat(ray_vdir, nbrep, axis=1) - pts = np.repeat(ray_orig, nbrep, axis=1) + k[None,:]*usbis + pts = np.repeat(ray_orig, nbrep, axis=1) + k[None, :]*usbis # memory view: reseff_mv = reseff indbis = np.concatenate(([0],ind,[k.size])) @@ -3025,10 +3027,10 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, coeff_ptr[0] = NULL reseff_arr = malloc(nlos*sizeof(double)) ind_arr = malloc(nlos*sizeof(long)) - # we sample lines of sight + # .. we sample lines of sight ...................................... _st.los_get_sample_core_var_res(nlos, - &lims[0,0], - &lims[1,0], + &lims[0, 0], + &lims[1, 0], n_dmode, n_imode, &res_arr[0], &coeff_ptr[0], @@ -3051,31 +3053,27 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, coeff_ptr[0], ind_arr, num_threads) + # .................................................................. if ani: val_2d = func(pts, t=t, vect=-usbis, **fkwdargs) else: val_2d = func(pts, t=t, **fkwdargs) # Integrate - if method=='sum': - # # .. original version ..................................... - # indbis = np.concatenate(([0],ind,[k.size])) - # for ii in range(1,nlos): - # sig[:,ii] = np.sum(val_2d[:,indbis[ii]:indbis[ii+1]], - # axis=-1)*reseff_mv[ii] - # # .......................................................... - # Calling integration function - _st.integrate_sum_nlos(nlos, nt, - val_2d, - sig_mv, - ind_arr, - reseff_arr, - num_threads) + if n_imode == 0: # "sum" integration mode + # .. integrating function .......................................... + for ii in range(nlos): + if ii > 0: + sig[:,ii] = np.sum(val_2d[:, ind_arr[ii-1]:ind_arr[ii]], + axis=-1) * reseff_arr[ii] + else: + sig[:,0] = np.sum(val_2d[:, 0:ind_arr[0]], + axis=-1) * reseff_arr[0] # Cleaning up... free(coeff_ptr[0]) free(coeff_ptr) free(reseff_arr) free(ind_arr) - elif method=='simps': + elif n_imode == 1: # "simpson" integration mode for ii in range(nlos): jj = indbis[ii] jjp1 = indbis[ii+1] @@ -3083,7 +3081,7 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, loc_r = reseff_mv[ii] sig[:,ii] = scpintg.simps(val_mv, x=None, dx=loc_r, axis=-1) - else: + else: # Romberg integration mode for ii in range(nlos): sig[:,ii] = scpintg.romb(val_2d[:,indbis[ii]:indbis[ii+1]], dx=reseff_mv[ii], axis=1, show=False) @@ -3093,9 +3091,10 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, elif minim == 'memory': # loop over LOS and parallelize if ani: - if n_imode == 0: + if n_imode == 0: # sum integration mode for ii in range(nlos): - pts, usbis = _st.call_get_sample_single_ani(lims[0,0], lims[1,0], + pts, usbis = _st.call_get_sample_single_ani(lims[0, ii], + lims[1, ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], @@ -3106,9 +3105,10 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, for jj in range(nt): val = func(pts, t=ltime[jj], vect=-usbis, **fkwdargs) sig[jj, ii] = np.sum(val)*loc_eff_res[0] - elif n_imode == 1: + elif n_imode == 1: # simpson integration mode for ii in range(nlos): - pts, usbis = _st.call_get_sample_single_ani(lims[0,0], lims[1,0], + pts, usbis = _st.call_get_sample_single_ani(lims[0, ii], + lims[1, ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], @@ -3120,9 +3120,10 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, val = func(pts, t=ltime[jj], vect=-usbis, **fkwdargs) sig[jj, ii] = scpintg.simps(val, x=None, dx=loc_eff_res[0]) - elif n_imode == 2: + elif n_imode == 2: # romberg integration mode for ii in range(nlos): - pts, usbis = _st.call_get_sample_single_ani(lims[0,0], lims[1,0], + pts, usbis = _st.call_get_sample_single_ani(lims[0, ii], + lims[1, ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], @@ -3135,10 +3136,11 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, sig[jj, ii] = scpintg.romb(val, show=False, dx=loc_eff_res[0]) else: - # -- not anisotropic ------------------------------------------------------ - if n_imode == 0: + # -- not anisotropic ----------------------------------------------- + if n_imode == 0: # "sum" integration mode for ii in range(nlos): - pts = _st.call_get_sample_single(lims[0,0], lims[1,0], + pts = _st.call_get_sample_single(lims[0, ii], + lims[1, ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], @@ -3149,9 +3151,10 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, for jj in range(nt): val = func(pts, t=ltime[jj], **fkwdargs) sig[jj, ii] = np.sum(val)*loc_eff_res[0] - elif n_imode == 1: + elif n_imode == 1: # "simpson" integration mode for ii in range(nlos): - pts = _st.call_get_sample_single(lims[0,0], lims[1,0], + pts = _st.call_get_sample_single(lims[0,ii], + lims[1,ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], @@ -3163,9 +3166,10 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, val = func(pts, t=ltime[jj], **fkwdargs) sig[jj, ii] = scpintg.simps(val, x=None, dx=loc_eff_res[0]) - elif n_imode == 2: + elif n_imode == 2: # "romberg" integration mode for ii in range(nlos): - pts = _st.call_get_sample_single(lims[0,0], lims[1,0], + pts = _st.call_get_sample_single(lims[0, ii], + lims[1, ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], @@ -3184,55 +3188,52 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, else: # loop over LOS if ani: - if n_imode == 0: - print("about to use new method") + if n_imode == 0: # sum integration mode for ii in range(nlos): - pts, usbis = _st.call_get_sample_single_ani(lims[0,0], lims[1,0], + 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]) val_2d = func(pts, t=t, vect=-usbis, **fkwdargs) - # this is almost always the quickest solution... but can - # probably be better. We'll investigate some time - # how to make it faster, and for the time being we leave - # the numpy alternative commented - _st.integrate_c_sum_mat(val_2d, &sig_mv[0,ii], nt, - nb_rows[0], - loc_eff_res[0], num_threads) - # sig_mv[:, ii] = np.sum(val, axis=-1)*loc_eff_res[0] - elif n_imode == 1: + sig_mv[:, 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,0], lims[1,0], + pts, usbis = _st.call_get_sample_single_ani(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, vect=-usbis, **fkwdargs) # integration sig[:, ii] = scpintg.simps(val, x=None, axis=-1, dx=loc_eff_res[0]) - elif n_imode == 2: + elif n_imode == 2: # romberg integration mode for ii in range(nlos): - pts, usbis = _st.call_get_sample_single_ani(lims[0,0], lims[1,0], + 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]) val = func(pts, t=t, vect=-usbis, **fkwdargs) sig[:, ii] = scpintg.romb(val, show=False, axis=1, dx=loc_eff_res[0]) else: - # -- not anisotropic ------------------------------------------------------ - if n_imode == 0: + # -- not anisotropic ----------------------------------------------- + if n_imode == 0: # "sum" integration mode for ii in range(nlos): - pts = _st.call_get_sample_single(lims[0,0], lims[1,0], + pts = _st.call_get_sample_single(lims[0,ii], lims[1,ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], @@ -3240,17 +3241,10 @@ 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, **fkwdargs) - # this is almost always the quickest solution... but can - # probably be better. We'll investigate some time - # how to make it faster, and for the time being we leave - # the numpy alternative commented - _st.integrate_c_sum_mat(val_2d, &sig_mv[0,ii], - nt, nb_rows[0], - loc_eff_res[0], num_threads) - # sig[:, ii] = np.sum(val,axis=-1)*loc_eff_res[0] - elif n_imode == 1: + 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,0], lims[1,0], + pts = _st.call_get_sample_single(lims[0,ii], lims[1,ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], @@ -3260,9 +3254,9 @@ def LOS_calc_signal(func, double[:,::1] ray_orig, double[:,::1] ray_vdir, res, val = func(pts, t=t, **fkwdargs) sig[:, ii] = scpintg.simps(val, x=None, axis=-1, dx=loc_eff_res[0]) - elif n_imode == 2: + elif n_imode == 2: # "romberg" integration mode for ii in range(nlos): - pts = _st.call_get_sample_single(lims[0,0], lims[1,0], + pts = _st.call_get_sample_single(lims[0,ii], lims[1,ii], res_mv[ii], n_dmode, n_imode, &loc_eff_res[0], diff --git a/tofu/geom/_basic_geom_tools.pxd b/tofu/geom/_basic_geom_tools.pxd index 7e91b3d26..99c5bb054 100644 --- a/tofu/geom/_basic_geom_tools.pxd +++ b/tofu/geom/_basic_geom_tools.pxd @@ -18,15 +18,6 @@ from cpython.array cimport array, clone cdef double _VSMALL cdef double _SMALL -cdef extern from "_fast_sum.c": - void sum_rows_blocks(double *orig, double *out, int n_rows, int n_cols) nogil - -cdef extern from "_fast_sum.c": - void sum_par_mat(double *orig, double *out, int n_rows, int n_cols) nogil - -cdef extern from "_fast_sum.c": - double sum_par_one_row(double *orig, int n_rows) nogil - # ============================================================================== # == Redifinition of functions # ============================================================================== @@ -107,11 +98,4 @@ cdef void compute_diff_div(const double[:, ::1] vec1, # ============================================================================== # == Matrix sum (np.sum) # ============================================================================== -cdef void sum_by_rows(double *orig, double *out, - int n_rows, int n_cols) nogil - -cdef void sum_naive_rows(double* orig, double* out, - int n_rows, int n_cols) nogil -cdef double sum_naive(double* orig, int n_cols) nogil - cdef long sum_naive_int(long* orig, int n_cols) nogil diff --git a/tofu/geom/_basic_geom_tools.pyx b/tofu/geom/_basic_geom_tools.pyx index dd87af10e..edcb4f502 100644 --- a/tofu/geom/_basic_geom_tools.pyx +++ b/tofu/geom/_basic_geom_tools.pyx @@ -370,23 +370,6 @@ cdef inline void sum_by_rows(double *orig, double *out, # ........... -cdef inline void sum_naive_rows(double* orig, double* out, - int n_rows, int n_cols) nogil: - cdef int ii, jj - for ii in prange(n_rows): - out[ii] = 0 - for jj in range(n_cols): - out[ii] += orig[ii*n_cols + jj] - - return - - -cdef inline double sum_naive(double* orig, int n_cols) nogil: - cdef int ii - cdef double out = 0. - for ii in prange(n_cols): - out += orig[ii] - return out cdef inline long sum_naive_int(long* orig, int n_cols) nogil: cdef int ii diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 103d3dc77..e429bf6c9 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -5886,7 +5886,7 @@ def calc_signal( self, func, t=None, - ani=None, + ani=False, fkwdargs={}, Brightness=True, res=None, @@ -5981,6 +5981,7 @@ def calc_signal( num_threads=num_threads, Test=True, ) + c0 = ( reflections and self._dgeom["dreflect"] is not None @@ -6013,11 +6014,11 @@ def calc_signal( ) # Integrate + # Creating the arrays with null everywhere.......... if s.ndim == 2: sig = np.full((s.shape[0], self.nRays), np.nan) else: sig = np.full((1, self.nRays), np.nan) - if t is None or len(t) == 1: sig[0, indok] = s else: @@ -6033,6 +6034,7 @@ def calc_signal( compact=True, pts=True, ) + if ani: nbrep = np.r_[ indpts[0], np.diff(indpts), pts.shape[1] - indpts[-1] @@ -6045,13 +6047,11 @@ def calc_signal( # This is the slowest step (~3.8 s with res=0.02 # and interferometer) val = func(pts, t=t, vect=vect) - # Integrate if val.ndim == 2: sig = np.full((val.shape[0], self.nRays), np.nan) else: sig = np.full((1, self.nRays), np.nan) - indpts = np.r_[0, indpts, pts.shape[1]] for ii in range(0, self.nRays): sig[:, ii] = ( @@ -6061,7 +6061,6 @@ def calc_signal( ) * reseff[ii] ) - # Format output return self._calc_signal_postformat( sig, diff --git a/tofu/geom/_fast_sum.c b/tofu/geom/_fast_sum.c deleted file mode 100644 index 3d4c5a7af..000000000 --- a/tofu/geom/_fast_sum.c +++ /dev/null @@ -1,56 +0,0 @@ -// C - algorithm for fast sum of a matrix (similar to np.sum(), sum row by row) -#include // for size_t - -#define MAX 8 -void sum_rows_blocks(double *orig, double *out, int n_rows, int n_cols){ - int n_blocks=n_rows/MAX; - int b; - int j; - int i; - for (b=0; blos_coeffs[0]) - pts = ray_orig + ksbis[None,:] * usbis + pts = ray_orig + ksbis[None, :] * usbis + # freeing memory used if los_coeffs != NULL: if los_coeffs[0] != NULL: free(los_coeffs[0]) @@ -1458,7 +1457,7 @@ cdef inline void los_get_sample_core_var_res(int nlos, return -# # -- utility for calc signal --------------------------------------------------- +# -- utility for calc signal --------------------------------------------------- cdef inline void los_get_sample_pts(int nlos, double* ptx, double* pty, @@ -1482,9 +1481,9 @@ cdef inline void los_get_sample_pts(int nlos, loc_vy = ray_vdir[1,0] loc_vz = ray_vdir[2,0] for ii in range(los_ind[0]): - ptx[ii] = loc_ox + coeff_ptr[ii] + loc_vx - pty[ii] = loc_oy + coeff_ptr[ii] + loc_vy - ptz[ii] = loc_oz + coeff_ptr[ii] + loc_vz + ptx[ii] = loc_ox + coeff_ptr[ii] * loc_vx + pty[ii] = loc_oy + coeff_ptr[ii] * loc_vy + ptz[ii] = loc_oz + coeff_ptr[ii] * loc_vz usx[ii] = loc_vx usy[ii] = loc_vy usz[ii] = loc_vz @@ -1497,71 +1496,10 @@ cdef inline void los_get_sample_pts(int nlos, loc_vy = ray_vdir[1,jj] loc_vz = ray_vdir[2,jj] for ii in range(los_ind[jj-1], los_ind[jj]): - ptx[ii] = loc_ox + coeff_ptr[ii] + loc_vx - pty[ii] = loc_oy + coeff_ptr[ii] + loc_vy - ptz[ii] = loc_oz + coeff_ptr[ii] + loc_vz + ptx[ii] = loc_ox + coeff_ptr[ii] * loc_vx + pty[ii] = loc_oy + coeff_ptr[ii] * loc_vy + ptz[ii] = loc_oz + coeff_ptr[ii] * loc_vz usx[ii] = loc_vx usy[ii] = loc_vy usz[ii] = loc_vz return - - -# -- calling sampling and intergrating with sum -------------------------------- -cdef inline void integrate_sum_nlos(int nlos, int nt, - double[:,::1] val_2d, - double[::1,:] sig_mv, - long* ind_arr, - double* reseff_arr, - int num_threads) nogil: - cdef int ii, jj - cdef int jjp1 - jj = 0 - jjp1 = ind_arr[0] - integrate_c_sum_mat(val_2d[:,jj:jjp1], - &sig_mv[0,0], - nt, jjp1 - jj, - reseff_arr[0], num_threads) - with nogil, parallel(num_threads=num_threads): - for ii in prange(1,nlos): - # sig[:,ii] = np.sum(val_2d[:,indbis[ii]:indbis[ii+1]], - # axis=-1)*reseff_mv[ii] - jj = ind_arr[ii-1] - jjp1 = ind_arr[ii] - integrate_c_sum_mat(val_2d[:,jj:jjp1], - &sig_mv[0,ii], - nt, jjp1 - jj, - reseff_arr[ii], num_threads) - return - -cdef inline void integrate_c_sum_mat(double[:,::1] val_mv, - double* sig, - int nrows, int ncols, - double loc_eff_res, - int num_threads) nogil: - cdef double* vsum - cdef int jj - # ... - vsum = malloc(nrows*sizeof(double)) - _bgt.sum_rows_blocks(&val_mv[0,0], &vsum[0], - nrows, ncols) - # _bgt.sum_by_rows(val_mv, &vsum[0], - # nrows, ncols) - # _bgt.sum_naive_rows(val_mv, &vsum[0], - # nrows, ncols) - # _bgt.sum_par_mat(val_mv, &vsum[0], - # nrows, ncols) - for jj in range(nrows): - sig[jj] = vsum[jj] * loc_eff_res - free(vsum) - return - - -cdef inline double integrate_c_sum_vec(double* val_mv, - int nsz, - double loc_eff_res, - int num_threads) nogil: - cdef double vsum - cdef int jj - # ... - vsum = _bgt.sum_par_one_row(val_mv, nsz) - return vsum * loc_eff_res diff --git a/tofu/tests/tests01_geom/tests03_core.py b/tofu/tests/tests01_geom/tests03_core.py index 0e0c2ebdb..8fdd22356 100644 --- a/tofu/tests/tests01_geom/tests03_core.py +++ b/tofu/tests/tests01_geom/tests03_core.py @@ -852,28 +852,37 @@ def ffT(Pts, t=None, vect=None): ind = None#[0,10,20,30,40] minimize = ["memory", "calls", "hybrid"] - for aa in [True, False]: - for rm in ["abs", "rel"]: + for typ in self.dobj.keys(): + c = 'CamLOS1D' + obj = self.dobj[typ][c] + for aa in [True, False]: + rm = 'rel' + # for rm in ["abs", "rel"]: + sigref, ii = None, 0 for dm in ["simps", "romb", "sum"]: for mmz in minimize: - for typ in self.dobj.keys(): - for c in self.dobj[typ].keys(): - obj = self.dobj[typ][c] - ff = ffT if obj.config.Id.Type=='Tor' else ffL - t = np.arange(0,10,10) - connect = (hasattr(plt.get_current_fig_manager(),'toolbar') - and getattr(plt.get_current_fig_manager(),'toolbar') - is not None) - out = obj.calc_signal(ff, t=t, ani=aa, - fkwdargs={}, - res=0.01, DL=None, - resMode=rm, - method=dm, minimize=mmz, - ind=ind, - plot=False, out=np.ndarray, - fs=(12,6), connect=connect) - sig, units = out - assert not np.all(np.isnan(sig)), str(ii) + ff = ffT if obj.config.Id.Type == 'Tor' else ffL + t = np.arange(0, 10, 10) + connect = (hasattr(plt.get_current_fig_manager(), + 'toolbar') + and getattr(plt.get_current_fig_manager(), + 'toolbar') + is not None) + out = obj.calc_signal(ff, t=t, ani=aa, + fkwdargs={}, + res=0.01, DL=None, + resMode=rm, + method=dm, minimize=mmz, + ind=ind, + plot=False, out=np.ndarray, + fs=(12, 6), connect=connect) + sig, units = out + assert not np.all(np.isnan(sig)), str(ii) + if sigref is not None: + assert np.allclose(sig, sigref) + if obj.nRays <= 100 and ii == 0: + sigref = sig + ii += 1 plt.close('all') def test11_plot(self): @@ -892,9 +901,10 @@ def test11_plot(self): lax = obj.plot(proj='hor', element='LDIO', Leg='KD', draw=False) except Exception as err: - msg = str(err) - msg += typ+' '+c - print(msg) + pass + # msg = str(err) + # msg += typ+' '+c + # print(msg) plt.close('all') def test12_plot_sino(self): @@ -927,6 +937,7 @@ def test14_saveload(self, verb=False): # Just to check the loaded version works fine obj2.strip(0, verb=verb) os.remove(pfe) + def test15_get_sample_same_res_unit(self): dmeths = ['rel', 'abs'] qmeths = ['simps', 'romb', 'sum'] diff --git a/tofu/tests/tests09_tutorials/__init__.py b/tofu/tests/tests09_tutorials/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tofu/tests/tests09_tutorials/tests01_runall.py b/tofu/tests/tests09_tutorials/tests01_runall.py new file mode 100644 index 000000000..969554fa6 --- /dev/null +++ b/tofu/tests/tests09_tutorials/tests01_runall.py @@ -0,0 +1,175 @@ + +# External modules +import os +import shutil +import types +import subprocess +import numpy as np +import warnings as warn +import matplotlib.pyplot as plt + +# Nose-specific +from nose import with_setup # optional + +# Importing package tofu.geom +import tofu as tf +from tofu import __version__ + + +_HERE = os.path.abspath(os.path.dirname(__file__)) +_TFROOT = _HERE[:-_HERE[::-1].index('/tofu'[::-1])-len('/tofu')] +# _TFROOT = tf.__path__[0][:-5] +# _PATHTUTO = os.path.join(_TFROOT, 'examples', 'tutorials') +keyVers = 'Vers' + + +####################################################### +# +# Setup and Teardown +# +####################################################### + +def setup_module(module, here=_HERE): + print("") # this is to get a newline after the dots + lf = os.listdir(here) + lf = [f for f in lf + if all([s in f for s in ['.npz']])] + lF = [] + for f in lf: + ff = f.split('_') + v = [fff[len(keyVers):] for fff in ff + if fff[:len(keyVers)] == keyVers] + msg = f + "\n " + str(ff) + "\n " + str(v) + assert len(v) == 1, msg + v = v[0] + if '.npz' in v: + v = v[:v.index('.npz')] + # print(v, __version__) + if v != __version__: + lF.append(f) + if len(lF) > 0: + print("Removing the following previous test files:") + for f in lF: + os.remove(os.path.join(here, f)) + # print("setup_module before anything in this file") + + +def teardown_module(module, here=_HERE): + # os.remove(VesTor.Id.SavePath + VesTor.Id.SaveName + '.npz') + # os.remove(VesLin.Id.SavePath + VesLin.Id.SaveName + '.npz') + # print("teardown_module after everything in this file") + # print("") # this is to get a newline + lf = os.listdir(here) + lf = [f for f in lf + if all([s in f for s in ['.npz']])] + lF = [] + for f in lf: + ff = f.split('_') + v = [fff[len(keyVers):] for fff in ff + if fff[:len(keyVers)] == keyVers] + msg = f + "\n "+str(ff) + "\n " + str(v) + assert len(v) == 1, msg + v = v[0] + if '.npz' in v: + v = v[:v.index('.npz')] + # print(v, __version__) + if v == __version__: + lF.append(f) + if len(lF) > 0: + print("Removing the following test files:") + for f in lF: + os.remove(os.path.join(here, f)) + + +# def my_setup_function(): +# print ("my_setup_function") + +# def my_teardown_function(): +# print ("my_teardown_function") + +# @with_setup(my_setup_function, my_teardown_function) +# def test_numbers_3_4(): +# print 'test_numbers_3_4 <============================ actual test code' +# assert multiply(3,4) == 12 + +# @with_setup(my_setup_function, my_teardown_function) +# def test_strings_a_3(): +# print 'test_strings_a_3 <============================ actual test code' +# assert multiply('a',3) == 'aaa' + + +####################################################### +# +# Struct subclasses +# +####################################################### + +def get_list_toturials(path=_HERE): + ltut = os.listdir(path) + ltut = [tt[:-3] for tt in ltut + if (all([ss in tt for ss in ['tuto', '.py']]) + and not any([ss in tt for ss in ['__init__', '.swp']]))] + return ltut + + +class Test00_tuto(object): + + @classmethod + def setup_class(cls, here=_HERE, pathtuto=_HERE, root=_TFROOT): + # print("") + # print("---- "+cls.__name__) + # ii = 0 + # for tuto in get_list_toturials(path=pathtuto): + # # Create local temporary copy to make sure the + # local version of tofu is imported + # method = 'test{0:02.0f}_{1}'.format(ii, tuto) + # fmethod = lambda tuto=tuto: cls.Test_tuto._test_tuto(tuto) + # setattr(cls, method, types.MethodType(fmethod, cls)) + # # setattr(cls, method, lambda cls: cls._test_tuto(tuto)) + # ii += 1 + pass + + @classmethod + def teardown_class(cls): + # print("teardown_class() after any methods in this class") + pass + + def setup(self): + # print("TestUM:setup() before each test method") + # self.setup_class() + pass + + def teardown(self): + # print("TestUM:teardown() after each test method") + pass + + @classmethod + 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) + try: + cmd = 'python ' + target + out = subprocess.run(cmd, shell=True, check=True, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE) + except Exception as err: + raise err + plt.close('all') + + # Remove temporary files and saved files + os.remove(target) + stdout = out.stdout.decode().split('\n') + lf = [stdout[ii].strip() for ii in range(len(stdout)) + if 'saved' in stdout[ii-1].lower()] + for ii in range(len(lf)): + os.remove(lf[ii]) + + def test01_plot_basic_tutorial(self): + self._test_tuto('tuto_plot_basic') + + def test02_plot_create_geometry(self): + self._test_tuto('tuto_plot_create_geometry') + + def test03_plot_custom_emissivity(self): + self._test_tuto('tuto_plot_custom_emissivity') diff --git a/tofu/tests/tests09_tutorials/tuto_plot_basic.py b/tofu/tests/tests09_tutorials/tuto_plot_basic.py new file mode 100644 index 000000000..50f99b86a --- /dev/null +++ b/tofu/tests/tests09_tutorials/tuto_plot_basic.py @@ -0,0 +1,78 @@ +""" +Getting started: 5 minutes tutorial for `tofu` +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This is a tutorial that aims to get a new user a little familiar with tofu's + structure. +""" + +# The following imports matplotlib, preferably using a +# backend that allows the plots to be interactive (Qt5Agg). +import matplotlib +try: + matplotlib.use('Qt5Agg') +except ImportError: + matplotlib.use(matplotlib.rcParams['backend']) + +############################################################################### +# We start by loading `tofu`. You might see some warnings at this stage since +# optional modules for `tofu` could +# be missing on the machine you are working on. This can be ignored safely. + +import numpy as np +import tofu as tf + +############################################################################### +# We can now create our first configuration. +# In `tofu` speak, a configuration is the geometry +# of the device and its structures. +# `tofu` provides pre-defined ones for your to try, +# so we're going to do just that: + +configB2 = tf.geom.utils.create_config("B2") + +############################################################################### +# The configuration can easily be visualized using the `.plot()` method: + +configB2.plot() + +############################################################################### +# Since `tofu` is all about tomography, +# let's create a 1D camera and plot its output. + +cam1d = tf.geom.utils.create_CamLOS1D( + config=configB2, + P=[3.4, 0, 0], + N12=100, + F=0.1, + D12=0.1, + angs=[np.pi, 0, 0], + Name="", + Exp="", + Diag="", +) +# interactive plot +cam1d.plot_touch() + +############################################################################### +# The principle is similar for 2D cameras. + +cam2d = tf.geom.utils.create_CamLOS2D( + config=configB2, + P=[3.4, 0, 0], + N12=100, + F=0.1, + D12=0.1, + angs=[np.pi, 0, 0], + Name="", + Exp="", + Diag="", +) +cam2d.plot_touch() + +############################################################################### +# What comes next is up to you! +# You could now play with the function parameters (change the cameras +# direction, refinement, aperture), +# with the plots (many are interactive) or create you own tomographic +# configuration. diff --git a/tofu/tests/tests09_tutorials/tuto_plot_create_geometry.py b/tofu/tests/tests09_tutorials/tuto_plot_create_geometry.py new file mode 100644 index 000000000..dcd3a586e --- /dev/null +++ b/tofu/tests/tests09_tutorials/tuto_plot_create_geometry.py @@ -0,0 +1,166 @@ +""" +Creating your first Geometry +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This is a tutorial that shows you how to create a simple geometry. +""" + +import numpy as np +import matplotlib.pyplot as plt +import tofu.geom as tfg + +import os +import tofu as tf +print(os.getcwd()) +print(tf.__version__) +print(tf.__path__) + + +############################################################################### +# Creating an empty Vessel +# ------------------------ +# +# If a vessel object does not exist yet, you have to create one (otherwise you +# can just load an existing one). +# A vessel object is basically defined by a 2D simple polygon +# (i.e.: non self-intersecting), that is then expanded linearly or toroidally +# depending on the desired configuration. +# This polygon limits the volume available for the plasma, where the emissivity +# can be non-zero. It is typically defined by the inner wall in a tokamak. +# +# Let's define the polygon limiting the vessel as a circle with a divertor-like +# shape at the bottom: + + +# Define the center, radius and lower limit +R0, Z0, rad, ZL = 2.0, 0.0, 1.0, -0.85 +# Define the key points in the divertor region below ZL +Div_R, Div_Z = [R0 - 0.2, R0, R0 + 0.2], [-1.2, -0.9, -1.2] +# Find the angles corresponding to ZL and span the rest +thet1 = np.arcsin((ZL - Z0) / rad) +thet2 = np.pi - thet1 +thet = np.linspace(thet1, thet2, 100) +# Assemble the polygon +poly_R = np.append(R0 + rad * np.cos(thet), Div_R) +poly_Z = np.append(Z0 + rad * np.sin(thet), Div_Z) +# Plot for checking +plt.figure(facecolor="w", figsize=(6, 6)) +plt.plot(poly_R, poly_Z) +plt.axis("equal") + +############################################################################### +# Notice that the polygon does not have to be closed, ToFu will anyway check +# that and close it automatically if necessary. +# Now let's feed this 2D polygon to the appropriate ToFu class and specify that +# it should be a toroidal type (if linear type 'Lin' is chosen, the length +# should be specified by the 'Lim' keyword argument). +# **tofu** also asks for a name to be associated to this instance, and an +# experiment ('Exp') and a shot number (useful when the same experiment changes +# geometry in time). + +# Create a toroidal Ves instance with name 'MyFirstVessel', associated to +# the experiment 'Misc' (for 'Miscellaneous') and shot number 0 +ves = tfg.Ves( + Name="MyFirstVessel", + Poly=[poly_R, poly_Z], + Type="Tor", + Exp="Misc", + shot=0 +) + +############################################################################### +# Now the vessel instance is created. It provides you several key attributes +# and methods (see :class:`~tofu.geom` for details). +# Among them the Id attribute is itself a class instance that contains all +# useful information about this vessel instance for identification, saving... +# In particular, that's where the name, the default saving path, the Type, the +# experiment, the shot number... are all stored. +# A default name for saving was also created that automatically includes not +# only the name you gave but also the module from which this instance was +# created (tofu.geom or tfg), the type of object, the experiment, the shot +# number... +# This recommended default pattern is useful for quick identification of saved +# object, it is advised not to modify it. + +print(ves.Id.SaveName) + +############################################################################### +# Now, we can simply visualise the created vessel by using the dedicated method +# (keyword argument 'Elt' specifies the elements of the instance we want to +# plot, typically one letter corresponds to one element, here we just want the +# polygon): + +# Plot the polygon by default in two projections (cross-section and horizontal) +# and return the list of axes +lax = ves.plot(element="P") + + +############################################################################### +# The created vessel instance, plotted in cross-section and horizontal +# projections +# Since the vessel is an important object (it defines where the plasma lives), +# all the other ToFu objects rely on it. It is thus important that you save +# it so that it can be used by other ToFu objects when necessary. + +ves.save(path="./") + +############################################################################### +# This method will save the instance as a numpy compressed file (.npz), using +# the path and file name found in ves.Id.SavePath and ves.Id.SaveName. While +# it is highly recommended to stick to the default value for the SaveName, +# but you can easily modify the saving path if you want by specifying it using +# keyword argument Path. + +############################################################################### +# Adding structural elements +# --------------------------- +# +# Like for a vessel, a structural element is mostly defined by a 2D polygon. +# If a vessel instance is provided, the type of the structural element +# (toroidal or linear) is automatically the same as the type of the vessel, +# otherwise the type must be specified. + +# A configuration, short for geometrical configuration is a set of vessel, +# and structural elements. + +# Define two polygons, one that does not enclose the vessel and one that does +thet = np.linspace(0.0, 2.0 * np.pi, 100) +poly1 = [[2.5, 3.5, 3.5, 2.5], [0.0, 0.0, 0.5, 0.5]] # a rectangle +poly2 = [R0 + 0.5 * np.cos(thet), -1.0 + 0.5 * np.sin(thet)] # a circle +poly3 = [[0.8, 1.3, 1.3, 0.8], [-0.5, -0.5, 0.5, 0.5]] # another rectangle +# Create the structural elements with the appropriate ToFu class, specifying +# the experiment and a shot number for keeping track of changes +s1 = tfg.PFC(Name="S1", + Poly=poly1, + Exp="Misc", + shot=0) +# now we create a structure that is not continuous along phi +# but is only defined within certain limits +s2 = tfg.PFC( + Name="S2", + Poly=poly2, + Exp="Misc", + shot=0, + Lim=[[0.0, np.pi], [np.pi / 2.0, np.pi * 3.0 / 2.0]], +) +# and another one, now defined as repetitions centered a position `pos` +# and with a certain `extent` +# here we wanted a structure uniformly repeated 5 times along phi +s3 = tfg.PFC( + Name="S3", + Poly=poly3, + # we dont take the last element of the list as it is 2.*pi = 0 + pos=np.linspace(0.0, 2.0 * np.pi, 6)[:-1], + # 5 repetitions + 5 empty spaces = 10 subdivision of 2. pi + extent=np.pi * 2.0 / 10.0, + Exp="Misc", + shot=0, +) +# Creating a configuration with vessel and structures +config = tfg.Config(Name="test", + Exp="Misc", + lStruct=[ves, s1, s2, s3]) +config.set_colors_random() # to see different colors +config.plot() +config.save() +# sphinx_gallery_thumbnail_number = 3 diff --git a/tofu/tests/tests09_tutorials/tuto_plot_custom_emissivity.py b/tofu/tests/tests09_tutorials/tuto_plot_custom_emissivity.py new file mode 100644 index 000000000..1da54ba20 --- /dev/null +++ b/tofu/tests/tests09_tutorials/tuto_plot_custom_emissivity.py @@ -0,0 +1,95 @@ +""" +Computing a camera image with custom emissivity +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This tutorial defines an emissivity that varies in space and computes +the signal received by a camera using this emissivity. +""" + +############################################################################### +# We start by loading a built-in `tofu` configuration and define a 2D camera. + +import matplotlib.pyplot as plt +import numpy as np +import tofu as tf + +configB2 = tf.geom.utils.create_config("B2") + +cam2d = tf.geom.utils.create_CamLOS2D( + config=configB2, + P=[3.4, 0, 0], + N12=100, + F=0.1, + D12=0.1, + angs=[np.pi, np.pi/6, 0], + Name="", + Exp="", + Diag="", +) + +############################################################################### +# 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) + :param t: optional, time parameter to add a time dependency to the + emissivity function + :param vect: + :return: + - emissivity -- 2D array holding the emissivity for each point in the + input grid + """ + r, z = np.hypot(pts[0, :], pts[1, :]), pts[2, :] + e = np.exp(-(r - 2.4) ** 2 / 0.2 ** 2 - z ** 2 / 0.4 ** 2) + if t is not None: + e = np.cos(np.atleast_1d(t))[:, None] * e[None, :] + else: + # as stated in documentation of calc_signal, e.ndim must be 2 + e = np.reshape(e, (1, -1)) + return e + + +y = np.linspace(2, 3, num=90) +z = np.linspace(-0.5, 0.5, num=100) +Y, Z = np.meshgrid(y, z) +X = np.zeros_like(Y) +pts = np.c_[X.ravel(), Y.ravel(), Z.ravel()].T +emissivity_vals = emissivity(pts) +emissivity_vals = emissivity_vals.reshape(X.shape) + + +def project_to_2D(xyz): + """Projection to (0, X, Z) plane.""" + return xyz[0], xyz[2] + + +fig, ax = plt.subplots() +ax.pcolormesh(Y, Z, emissivity_vals) +ax.set_xlabel('y') +ax.set_ylabel('z') +configB2.plot(lax=ax, proj='cross') +cam_center, = ax.plot(*project_to_2D(cam2d._dgeom['pinhole']), '*', ms=20) +ax.legend(handles=[cam_center], labels=['camera pinhole'], loc='upper right') + +############################################################################### +# Finally, we compute an image using the 2D camera and this emissivity. +# If we provide a time vector, the field will vary in a cosinusoidal fashion +# (see above definition) across time. + +time_vector = np.linspace(0, 2 * np.pi, num=100) + +sig, units = cam2d.calc_signal(emissivity, + res=0.01, + reflections=False, + minimize="hybrid", + method="sum", + newcalc=True, + plot=False, + t=time_vector) + +sig.plot(ntMax=1) +plt.show(block=False) diff --git a/tofu/version.py b/tofu/version.py index 1fc557b2f..ac22e7403 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-217-g754c8bb' +__version__ = '1.4.1-232-g56df4eb'