add a gaussian_filter and find_peaks for scalarfield and vectorfield

pyphare/pyphare/core/gridlayout.py

@@ -44,6 +44,10 @@
         "Pyz": "primal",
         "Pzz": "primal",
         "tags": "dual",
+        "value": "primal",
+        "x": "primal",
+        "y": "primal",
+        "z": "primal",
     },
     "y": {
         "Bx": "dual",
@@ -79,6 +83,10 @@
         "Pyz": "primal",
         "Pzz": "primal",
         "tags": "dual",
+        "value": "primal",
+        "x": "primal",
+        "y": "primal",
+        "z": "primal",
     },
     "z": {
         "Bx": "dual",
@@ -114,6 +122,10 @@
         "Pyz": "primal",
         "Pzz": "primal",
         "tags": "dual",
+        "value": "primal",
+        "x": "primal",
+        "y": "primal",
+        "z": "primal",
     },
 }
 yee_centering_lower = {

pyphare/pyphare/core/operators.py

@@ -20,12 +20,12 @@ def _compute_dot_product(patch_datas, **kwargs):
 
 
 def _compute_sqrt(patch_datas, **kwargs):
-    ref_name = next(iter(patch_datas.keys()))
+    # ref_name = next(iter(patch_datas.keys())) TODO this is always "value"
 
     dset = np.sqrt(patch_datas["value"][:])
 
     return (
-        {"name": "value", "data": dset, "centering": patch_datas[ref_name].centerings},
+        {"name": "value", "data": dset, "centering": patch_datas["value"].centerings},
     )
 
 
@@ -139,3 +139,4 @@ def grad(hier, **kwargs):
     h = compute_hier_from(_compute_grad, hier, nb_ghosts=nb_ghosts)
 
     return VectorField(h)
+

pyphare/pyphare/core/ufuncs.py

@@ -0,0 +1,170 @@
+import numpy as np
+
+from pyphare.pharesee.hierarchy import ScalarField, VectorField
+from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from
+
+from scipy.ndimage import gaussian_filter
+
+
+
+def _compute_gaussian_filter_on_scalarfield(patch_datas, **kwargs):
+    from scipy.ndimage import gaussian_filter
+
+    ndim = patch_datas["value"].box.ndim
+    nb_ghosts = kwargs["nb_ghosts"]
+    sigma = kwargs["sigma"]
+    ds = np.asarray(patch_datas["value"][:])
+
+    ds_ = np.full(list(ds.shape), np.nan)
+
+    gf_ = gaussian_filter(ds, sigma=sigma)
+    select = tuple([slice(nb_ghosts, -nb_ghosts) for _ in range(ndim)])
+    ds_[select] = np.asarray(gf_[select])
+
+    return (
+        {"name": "value", "data": ds_, "centering": patch_datas["value"].centerings},
+    )
+
+
+def _compute_gaussian_filter_on_vectorfield(patch_datas, **kwargs):
+    from scipy.ndimage import gaussian_filter
+
+    ref_name = next(iter(patch_datas.keys()))
+
+    ndim = patch_datas[ref_name].box.ndim
+    nb_ghosts = kwargs["nb_ghosts"]
+    sigma = kwargs["sigma"]
+    ds_x = np.asarray(patch_datas["x"][:])
+    ds_y = np.asarray(patch_datas["y"][:])
+    ds_z = np.asarray(patch_datas["z"][:])
+
+    dsx_ = np.full(list(ds_x.shape), np.nan)
+    dsy_ = np.full(list(ds_y.shape), np.nan)
+    dsz_ = np.full(list(ds_z.shape), np.nan)
+
+    gfx_ = gaussian_filter(ds_x, sigma=sigma)
+    gfy_ = gaussian_filter(ds_y, sigma=sigma)
+    gfz_ = gaussian_filter(ds_z, sigma=sigma)
+
+    select = tuple([slice(nb_ghosts, -nb_ghosts) for _ in range(ndim)])
+
+    dsx_[select] = np.asarray(gfx_[select])
+    dsy_[select] = np.asarray(gfy_[select])
+    dsz_[select] = np.asarray(gfz_[select])
+
+    return (
+        {"name": "x", "data": dsx_, "centering": patch_datas["x"].centerings},
+        {"name": "y", "data": dsy_, "centering": patch_datas["y"].centerings},
+        {"name": "z", "data": dsz_, "centering": patch_datas["z"].centerings},
+    )
+
+
+# def gFilt(hier, **kwargs):
+#     sigma = kwargs.get("sigma", 2)
+# 
+#     # time0 = list(hier.times())[0]
+#     # level0 = 0
+#     # p0 = 0
+#     # pd0 = hier.levels(time0)[level0].patches[p0].patch_datas
+#     # key0 = list(pd0.keys())[0]
+#     # nb_ghosts = np.max(pd0[key0].ghosts_nbr)
+# 
+#     nb_ghosts = np.max(list(hier.level(0).patches[0].patch_datas.values())[0].ghosts_nbr)
+# 
+#     if nb_ghosts < sigma :
+#         print("nb_ghosts ({0}) < sigma ({1}) : your gaussian filter might be dirty".format(nb_ghosts, sigma))
+# 
+#     if hier.ndim == 1:
+#         if isinstance(hier, ScalarField) :
+#             h = compute_hier_from(_compute_gaussian_filter_on_scalarfield, hier, nb_ghosts=nb_ghosts, sigma=sigma)
+#             return ScalarField(h)
+#         elif isinstance(hier, VectorField) :
+#             h = compute_hier_from(_compute_gaussian_filter_on_vectorfield, hier, nb_ghosts=nb_ghosts, sigma=sigma)
+#             return VectorField(h)
+#         else:
+#             return NotImplemented
+#     else:
+#         return NotImplemented
+
+
+def gF(hier, **kwargs):
+    sigma = kwargs.get("sigma", 4)
+    if sigma == 1:
+        raise ValueError("sigma value has to be at least 2")
+    h_ = hier.__deepcopy__(memo={})
+    ndim = hier.ndim
+    n_pad = 4*sigma+1
+    # The gaussian filter is calculated on the box extended by
+    # n_pad. Hence, the number of points is large enough so that the value
+    # at the last point of the real box is as equal as possible to the one
+    # at the first point of the next box...
+
+    for time in h_.times():
+        interp_ = hier.interpol(time)
+        for lvl in h_.levels(time).values():
+            for patch in lvl.patches:
+                names = list(patch.patch_datas.keys())
+                box = patch.box
+
+                for name in names:
+                    pdata = patch.patch_datas[name]
+                    nb_ghosts = pdata.ghosts_nbr
+                    if not n_pad > nb_ghosts:
+                        raise ValueError('sigma value is too small')
+
+                    r_ = []
+                    for i in range(ndim):
+                        s_ = np.arange(box.lower[i]-n_pad, box.upper[i]+2+n_pad)*pdata.dl[i]
+                        r_.append(s_)
+
+                    func, _ = interp_[name]
+
+                    if ndim == 1:
+                        data = func(r_[0])
+                    elif ndim == 2:
+                        data = func(r_[0], r_[1])
+                    elif ndim == 3:
+                        data = func(r_[0], r_[1], r_[2])
+                    else:
+                        raise ValueError('unvalid dimension')
+
+                    gf_ = gaussian_filter(data, sigma=sigma)
+                    select = tuple([slice(n_pad-nb_ghosts[i], -n_pad+nb_ghosts[i]) for i in range(ndim)])
+                    pdata.dataset = np.asarray(gf_[select])
+
+    return h_
+
+
+def peakIds(hier, **kwargs):
+    from scipy.signal import find_peaks
+
+    times = list(hier.times())
+    if len(times) == 1:
+        time = times[0]
+    else:
+        raise ValueError('multiple time is not yet implemented')
+
+    # pks_ = np.array([])
+    # hgs_ = np.array([])
+
+    names_ = kwargs.pop("names", None)
+    if names_ is None:
+        names_ = list(hier.levels(time)[0].patches[0].patch_datas.keys())
+
+    # ph_ = kwargs.get('peak_heights', None)
+    # if ph_ is None:
+    #     raise ValueError("the kwarg 'peak_heights' is mandatory for now...")
+
+    for lvl in hier.levels(time).values():
+        for patch in lvl.patches:
+            for name in names_:
+                pdata = patch.patch_datas[name]
+                ds = np.asarray(pdata.dataset)
+                # pks = find_peaks(ds, **kwargs)
+                # for pk, hg in zip(pks[0], pks[1]['peak_heights']):
+                #     pks_ = np.append(pks_, np.add(np.multiply(pk, patch.dl), patch.origin))
+                #     hgs_ = np.append(hgs_, hg)
+                out = find_peaks(ds, **kwargs)
+
+    # return pks_, hgs_
+    return zip(*out)

pyphare/pyphare/pharesee/hierarchy/hierarchy.py

@@ -438,8 +438,9 @@
                 label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip)
                 marker = kwargs.get("marker", "")
                 ls = kwargs.get("ls", "--")
+                lw = kwargs.get("lw", 1)
                 color = kwargs.get("color", "k")
-                ax.plot(x, val, label=label, marker=marker, ls=ls, color=color)
+                ax.plot(x, val, label=label, marker=marker, ls=ls, lw=lw, color=color)
 
         ax.set_title(kwargs.get("title", ""))
         ax.set_xlabel(kwargs.get("xlabel", "x"))
@@ -614,6 +615,128 @@
 
         return final, dp(final, **kwargs)
 
+    def interpol(self, time, interp="nearest"):
+        from pyphare.pharesee.hierarchy.hierarchy_utils import flat_finest_field
+        from pyphare.pharesee.run.utils import build_interpolator
+
+        nbrGhosts = list(self.level(0, time).patches[0].patch_datas.values())[0].ghosts_nbr
+
+        interp_ = {}
+        for qty in self.quantities():
+            box = self.level(0, time).patches[0].patch_datas[qty].box
+            dl = self.level(0, time).patches[0].patch_datas[qty].dl
+            data, coords = flat_finest_field(self, qty, time=time)
+            interp_[qty] = build_interpolator(data, coords, interp, box, dl, qty, nbrGhosts)
+        return interp_
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        print(f"__array_function__ of PatchHierarchy called for {ufunc.__name__}")
+        if method != "__call__":
+            return NotImplemented
+
+        final = []
+
+        times = inputs[0].times()
+        for x in inputs:
+            assert(times == x.times())
+            if not isinstance(x, PatchHierarchy):
+                raise TypeError("this arg should be a PatchHierarchy")
+
+
+        # for  d in inputs[0].patch_levels:
+        #     print(type(d), d)
+
+
+        ils = [key for dict_ in inputs[0].patch_levels for key in dict_]
+        h_type = type(inputs[0])
+
+        all_ = []
+        for i in range(len(times)):
+            pls = []
+            for x in inputs:
+                pls_ = []
+                for plvl in x.patch_levels[i].values():
+                    pls_.append(plvl)
+                pls.append(pls_)
+
+            out = [getattr(ufunc, method)(*pl, **kwargs)  for pl in zip(*pls)]
+
+            # out est une liste de liste de patchlevel : indice  sur le levels (pour 1 temps donne)
+            # il faut les remettre dans un dict avec ilvl
+
+            final = {}
+            for il, pl in zip(ils, out):
+                final[il] = pl
+
+            all_.append(final)
+
+        h_ = PatchHierarchy(all_,
+                           domain_box=self.domain_box,
+                           refinement_box=self.refinement_ratio,
+                           times=self.times(),
+                           data_files=self.data_files,
+                           selection_box=self.selection_box)
+
+        from .scalarfield import ScalarField
+        from .vectorfield import VectorField
+
+        if h_type is ScalarField:
+            return ScalarField(h_)
+        elif h_type is VectorField:
+            return VectorField(h_)
+
+    def __array_function__(self, func, types, args, kwargs):
+        print(f"__array_function__ of Patch {func.__name__} called for {[getattr(a, 'name', a) for a in args]}")
+
+        final = []
+        others = []
+
+        times = args[0].times()
+        for x in args:
+            if isinstance(x, PatchHierarchy):
+                assert(times == x.times())
+            else:
+                others.append(x)
+
+        ils = [key for dict_ in args[0].patch_levels for key in dict_]
+        h_type = type(args[0])
+
+
+        all_ = []
+        for i in range(len(times)):
+            pls = []
+            for x in args:
+                if isinstance(x, PatchHierarchy):
+                    pls_ = []
+                    for plvl in x.patch_levels[i].values():
+                        pls_.append(plvl)
+                    pls.append(pls_)
+
+            out = []
+            for pl in zip(*pls):
+                out.append(func(*pl, *others, **kwargs))
+
+            final = {}
+            for il, pl in zip(ils, out):
+                final[il] = pl
+
+            all_.append(final)
+
+        h_ = PatchHierarchy(all_,
+                           domain_box=self.domain_box,
+                           refinement_box=self.refinement_ratio,
+                           times=self.times(),
+                           data_files=self.data_files,
+                           selection_box=self.selection_box)
+
+        from .scalarfield import ScalarField
+        from .vectorfield import VectorField
+
+        if h_type is ScalarField:
+            return ScalarField(h_)
+        elif h_type is VectorField:
+            return VectorField(h_)
+
 
 def finest_part_data(hierarchy, time=None):
     """