From 0f0ae08caa8fc07aaa398c2d5b7261d97dabd18f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Fri, 20 May 2016 16:46:03 +0200
Subject: [PATCH]  - Add WorkSpaceTools class for internal work array
 management - clean/rename ... for pep8 compat

---
 CMakeLists.txt                                |   2 +-
 hysop/domain/mesh.py                          |   4 +-
 hysop/domain/submesh.py                       |   4 +-
 hysop/domain/tests/test_submesh.py            |   6 +-
 hysop/mpi/bridge.py                           |   6 +-
 hysop/mpi/bridge_inter.py                     |   8 +-
 hysop/mpi/bridge_overlap.py                   |   6 +-
 hysop/mpi/topology.py                         |   6 +-
 hysop/tools/misc.py                           | 196 +++++++++++++++++-
 hysop/tools/numpywrappers.py                  |  41 +++-
 hysop/tools/problem2dot.py                    |  24 ++-
 hysop/tools/remeshing_formula_parsing.py      |  35 ++--
 .../tests/{test_parameters.py => test_io.py}  |   0
 13 files changed, 288 insertions(+), 50 deletions(-)
 rename hysop/tools/tests/{test_parameters.py => test_io.py} (100%)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 782725fa1..9e85aa59a 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -463,7 +463,7 @@ add_custom_target(pyclean COMMAND rm -f ${PYCFILES}
 # ============= Tests =============
 if(WITH_TESTS)
   if(NOT ENABLE_LONG_TESTS)
-    set(TEST_TIMEOUT 4 CACHE STRING "Default value for test runtime limit (in seconds)")
+    set(TEST_TIMEOUT 10 CACHE STRING "Default value for test runtime limit (in seconds)")
   else()
     set(TEST_TIMEOUT 500 CACHE STRING "Default value for test runtime limit (in seconds)")
   endif()
diff --git a/hysop/domain/mesh.py b/hysop/domain/mesh.py
index 2f181729a..4b0ec7576 100644
--- a/hysop/domain/mesh.py
+++ b/hysop/domain/mesh.py
@@ -11,7 +11,7 @@ from hysop.constants import debug
 import hysop.tools.numpywrappers as npw
 from hysop.tools.parameters import Discretization
 import numpy as np
-from hysop.tools.misc import utils
+from hysop.tools.misc import Utils
 
 
 class Mesh(object):
@@ -272,7 +272,7 @@ class Mesh(object):
 
         Returns the same kind of list but in local coordinates
         """
-        tmp = utils.intersl(sl, self.position)
+        tmp = Utils.intersl(sl, self.position)
         if tmp is not None:
             return [slice(tmp[i].start + self._shift[i],
                           tmp[i].stop + self._shift[i])
diff --git a/hysop/domain/submesh.py b/hysop/domain/submesh.py
index dc9f307c8..46b0735fe 100644
--- a/hysop/domain/submesh.py
+++ b/hysop/domain/submesh.py
@@ -5,7 +5,7 @@
 import hysop.tools.numpywrappers as npw
 from hysop.tools.parameters import Discretization
 import numpy as np
-from hysop.tools.misc import utils
+from hysop.tools.misc import Utils
 
 
 class SubMesh(object):
@@ -70,7 +70,7 @@ class SubMesh(object):
         # find its computational points. Warning:
         # the indices of computational points must be
         # relative to the parent mesh local grid!
-        sl = utils.intersl(global_position_in_parent, self.mesh.position)
+        sl = Utils.intersl(global_position_in_parent, self.mesh.position)
         # Bool to check if this process holds the end (in any direction)
         # of the domain. Useful for proper integration on this subset.
         is_last = [False, ] * self._dim
diff --git a/hysop/domain/tests/test_submesh.py b/hysop/domain/tests/test_submesh.py
index 5df5aeb85..9f427d7bc 100755
--- a/hysop/domain/tests/test_submesh.py
+++ b/hysop/domain/tests/test_submesh.py
@@ -4,7 +4,7 @@ from hysop.domain.box import Box
 from hysop.fields.continuous import Field
 from hysop.tools.parameters import Discretization
 from hysop.domain.submesh import SubMesh
-from hysop.tools.misc import utils
+from hysop.tools.misc import Utils
 import numpy as np
 from hysop.mpi import main_rank
 
@@ -43,10 +43,10 @@ def check_submesh(discr):
     # of the submesh and the position of the parent mesh.
     if subm.on_proc:
         assert subm.position_in_parent == [s for s in
-                                           utils.intersl(gpos, mesh.position)]
+                                           Utils.intersl(gpos, mesh.position)]
 
         r2 = mesh.convert2local([s for s
-                                 in utils.intersl(gpos, mesh.position)])
+                                 in Utils.intersl(gpos, mesh.position)])
         assert subm.iCompute == r2
         pos = subm.position_in_parent
         subpos = [slice(pos[d].start - gstart[d],
diff --git a/hysop/mpi/bridge.py b/hysop/mpi/bridge.py
index b9765d94d..f07f590e5 100644
--- a/hysop/mpi/bridge.py
+++ b/hysop/mpi/bridge.py
@@ -4,7 +4,7 @@ Tools to compute the intersection between
 two HySoP topologies.
 """
 from hysop.mpi.topology import Cartesian, TopoTools
-from hysop.tools.misc import utils
+from hysop.tools.misc import Utils
 
 
 class Bridge(object):
@@ -79,7 +79,7 @@ class Bridge(object):
         # defined by the slices above.
         current = indices_source[self._rank]
         for rk in indices_target:
-            inter = utils.intersl(current, indices_target[rk])
+            inter = Utils.intersl(current, indices_target[rk])
             if inter is not None:
                 self._send_indices[rk] = inter
         # Back to local indices
@@ -96,7 +96,7 @@ class Bridge(object):
         # the grid points defined by the slices above.
         current = indices_target[self._rank]
         for rk in indices_source:
-            inter = utils.intersl(current, indices_source[rk])
+            inter = Utils.intersl(current, indices_source[rk])
             if inter is not None:
                 self._recv_indices[rk] = inter
 
diff --git a/hysop/mpi/bridge_inter.py b/hysop/mpi/bridge_inter.py
index d761a05fb..ecc0af52f 100644
--- a/hysop/mpi/bridge_inter.py
+++ b/hysop/mpi/bridge_inter.py
@@ -4,7 +4,7 @@ Tools to compute the intersection between
 two HySoP topologies.
 """
 from hysop.mpi.topology import Cartesian, TopoTools
-from hysop.tools.misc import utils
+from hysop.tools.misc import Utils
 from hysop.mpi import MPI
 import hysop.tools.numpywrappers as npw
 
@@ -73,7 +73,7 @@ class BridgeInter(object):
         self._tranfer_indices = {}
         current = current_indices[self._topology.rank]
         for rk in remote_indices:
-            inter = utils.intersl(current, remote_indices[rk])
+            inter = Utils.intersl(current, remote_indices[rk])
             if inter is not None:
                 self._tranfer_indices[rk] = inter
 
@@ -116,8 +116,8 @@ class BridgeInter(object):
                 self.comm.bcast(current_indices, root=MPI.PROC_NULL)
 
         # Convert numpy arrays to dict of slices ...
-        current_indices = utils.arrayToDict(current_indices)
-        remote_indices = utils.arrayToDict(remote_indices)
+        current_indices = Utils.arrayToDict(current_indices)
+        remote_indices = Utils.arrayToDict(remote_indices)
 
         return current_indices, remote_indices
 
diff --git a/hysop/mpi/bridge_overlap.py b/hysop/mpi/bridge_overlap.py
index 4dfa0b904..e4f75ba06 100644
--- a/hysop/mpi/bridge_overlap.py
+++ b/hysop/mpi/bridge_overlap.py
@@ -5,7 +5,7 @@ two HySoP topologies defined on the same comm but for a
 different number of processes.
 """
 from hysop.mpi.topology import Cartesian, TopoTools
-from hysop.tools.misc import utils
+from hysop.tools.misc import Utils
 from hysop.mpi.bridge import Bridge
 
 
@@ -90,7 +90,7 @@ class BridgeOverlap(Bridge):
             current = [slice(None, None, None), ] * dimension
 
         for rk in indices_target:
-            inter = utils.intersl(current, indices_target[rk])
+            inter = Utils.intersl(current, indices_target[rk])
             if inter is not None:
                 self._send_indices[rk] = inter
 
@@ -105,7 +105,7 @@ class BridgeOverlap(Bridge):
         else:
             current = [slice(None, None, None), ] * dimension
         for rk in indices_source:
-            inter = utils.intersl(current, indices_source[rk])
+            inter = Utils.intersl(current, indices_source[rk])
             if inter is not None:
                 self._recv_indices[rk] = inter
 
diff --git a/hysop/mpi/topology.py b/hysop/mpi/topology.py
index fc91dd946..1ab614715 100644
--- a/hysop/mpi/topology.py
+++ b/hysop/mpi/topology.py
@@ -13,7 +13,7 @@ from hysop.mpi.main_var import MPI
 from hysop.tools.parameters import Discretization, MPIParams
 import numpy as np
 import hysop.tools.numpywrappers as npw
-from hysop.tools.misc import utils
+from hysop.tools.misc import Utils
 
 
 class Cartesian(object):
@@ -508,7 +508,7 @@ class TopoTools(object):
                         root=root)
 
         if toslice:
-            return utils.arrayToDict(iglob_res)
+            return Utils.arrayToDict(iglob_res)
         else:
             return iglob_res
 
@@ -572,7 +572,7 @@ class TopoTools(object):
                 comm.Gather([iglob[:, rank], MPI.INT], [iglob_res, MPI.INT],
                             root=root)
             if toslice:
-                return utils.arrayToDict(iglob_res)
+                return Utils.arrayToDict(iglob_res)
             else:
                 return iglob_res
 
diff --git a/hysop/tools/misc.py b/hysop/tools/misc.py
index a5e102317..a065daba8 100644
--- a/hysop/tools/misc.py
+++ b/hysop/tools/misc.py
@@ -1,10 +1,23 @@
-"""Some utilities to deal with slices and other python objects
+"""Some utilities to deal with workspaces, slices and other python objects
+
+* :class:`~hysop.tools.misc.Utils` : some utilities to handle slices in hysop.
+
+* :class:`~hysop.tools.misc.WorkSpaceTools` to deal with operators' local
+  work spaces, see :ref:`work_spaces`.
+
+
 """
+import numpy as np
+import hysop.tools.numpywrappers as npw
+from hysop.constants import HYSOP_REAL, HYSOP_INTEGER
+
 
-class utils(object):
+class Utils(object):
+    """tools to handle array and slices.
+    """
 
     @staticmethod
-    def arrayToDict(inarray):
+    def array_to_dict(inarray):
         """
         convert an array into a dictionnary,
         keys being the column numbers in array
@@ -45,3 +58,180 @@ class utils(object):
                 return None
             res[d] = slice(start, stop)
         return tuple(res)
+
+    @staticmethod
+    def is_on_proc(sl):
+        """True if sl represent a non empty set of indices.
+        """
+        return (np.asarray([ind.stop != ind.start
+                            for ind in sl])).all()
+
+
+class WorkSpaceTools(object):
+    """Tools to deal with internal work arrays for operators
+    """
+
+    @staticmethod
+    def check_work_array(lwork, subshape, work=None,
+                         data_type=HYSOP_REAL):
+        """Check properties of existing working array or
+        allocate some new buffers complient with some properties.
+
+        Parameters
+        ----------
+        lwork : int
+            required number of working arrays
+        subshape : list of tuples of int
+            required shape for work array
+            subshape[i] == expected shape for work[i]
+        work : list of numpy arrays
+        data_type : either HYSOP_REAL or HYSOP_INTEGER
+
+        Notes
+        -----
+        work arrays are 1D arrays of size prod(subshape) that are 'reshaped'
+        according to the specific needs of each operator. That means
+        that one memory location (the 1D array) may be shared
+        between several operators thanks to the numpy reshape function.
+
+        """
+        result = []
+        #print subshape, len(subshape), np.prod(subshape[0])
+        if isinstance(subshape, list):
+            subsize = [np.prod(subshape[i]) for i in xrange(len(subshape))]
+        else:
+            subsize = [np.prod(subshape), ] * lwork
+            subshape = [subshape, ] * lwork
+        if work is None:
+            for i in xrange(lwork):
+                result.append(npw.zeros(subsize[i],
+                                        dtype=data_type).reshape(subshape[i]))
+        else:
+            assert isinstance(work, list), 'rwork must be a list.'
+            msg = 'Work arrays list must be at least of size '
+            msg += str(lwork) + '.'
+            assert len(work) >= lwork, msg
+            msg1 = 'Work array size is too small.'
+            msg2 = 'Work array must be a flat array (1D).'
+            try:
+                for i in xrange(lwork):
+                    wk = work[i]
+                    assert wk.size >= subsize[i], msg1
+                    assert len(np.where(
+                        np.asarray(wk.shape) > 1)[0]) == 1, msg2
+                    result.append(wk.ravel()[:subsize[i]].reshape(subshape[i]))
+
+            except AttributeError:
+                # Work array has been replaced by an OpenCL Buffer
+                # Testing the buffer size instead of shape
+                for i in xrange(lwork):
+                    wk = work[i]
+                    s = wk.size / subsize[i]
+                    WorkSpaceTools._check_ocl_buffer(s, data_type)
+
+                result = work
+        return result
+
+    @staticmethod
+    def _check_ocl_buffer(s, dtype):
+        """check if opencl buffer size is complient with input type.
+        """
+        if dtype is HYSOP_REAL:
+            assert (HYSOP_REAL is np.float32 and s == 4) or \
+                (HYSOP_REAL is np.float64 and s == 8)
+        elif dtype is HYSOP_INTEGER:
+            assert (HYSOP_INTEGER is np.int16 and s == 2) or \
+                (HYSOP_INTEGER is np.int32 and s == 4) or \
+                (HYSOP_INTEGER is np.int64 and s == 8)
+
+    @staticmethod
+    def find_common_workspace(operators, array_type='rwork'):
+        """Allocate a list of common workspaces able to work
+        with some given operators
+
+        Parameters
+        ----------
+        operators : a list of dictionnary
+            Each component of the list or value of the dictionnary
+            must be an object that may need some internal workspaces
+            and with a 'get_work_properties' function, hysop operators indeed.
+        array_type : string, optional
+            between 'rwork' (real arrays) or 'iwork' (integer arrays),
+            depending on the required work.
+            Default is 'rwork'.
+        Returns
+        -------
+        work : list of numpy arrays
+            workspaces common to all operators.
+
+        Example
+        -------
+
+        # op1, op2 some predifined operators
+        op1.discretize()
+        op2.discretize()
+        rwork = Utils.find_common_workspace([op1, op2])
+        iwork = Utils.find_common_workspace([op1, op2], 'iwork')
+        op1.setup(rwork=work, iwork=iwork)
+        op2.setup(rwork=work)
+        # work is a common internal list of workspaces that can be used by both
+        # operators.
+        """
+        properties = []
+        if isinstance(operators, dict):
+            oplist = operators.values()
+        elif isinstance(operators, list):
+            oplist = operators
+        else:
+            raise AttributeError('operators must be a list or a dict')
+
+        for op in oplist:
+            properties.append(op.get_work_properties()[array_type])
+
+        return WorkSpaceTools.allocate_common_workspaces(properties)
+
+    @staticmethod
+    def allocate_common_workspaces(work_properties):
+        """Allocate a list of common workspaces according to
+        some given properties.
+
+        Parameters
+        ----------
+        work_properties : list of lists of tuples
+            properties of workspaces. Each component of this list
+            must be the return value of a get_work_properties function
+            (for instance from an operator)
+
+        Returns
+        -------
+        work : list of numpy arrays
+            workspaces fitting with properties requirements.
+
+        Example
+        -------
+
+        # op1, op2 some predifined operators
+        op1.discretize()
+        op2.discretize()
+        wk_prop = []
+        wk_prop.append(op1.get_work_properties()['rwork'])
+        wk_prop.append(op2.get_work_properties()['rwork'])
+        work = Utils.find_common_workspace(wk_prop)
+
+        op1.setup(rwork=work)
+        op2.setup(rwork=work)
+        # work is a common internal list of workspaces that can be used by both
+        # operators.
+        """
+        assert isinstance(work_properties, list)
+        # Max number of required working arrays:
+        properties = [p for p in work_properties if p is not None]
+        lwork = max([len(prop) for prop in properties])
+        # Then find the max required workspace for work array
+        shapes = [(0,), ] * lwork
+        for prop in properties:
+            lp = len(prop)
+            shapes[:lp] = np.maximum(shapes[:lp], prop)
+        work = [npw.zeros(shape) for shape in shapes]
+
+        return work
diff --git a/hysop/tools/numpywrappers.py b/hysop/tools/numpywrappers.py
index b0868f27e..b31e7100c 100644
--- a/hysop/tools/numpywrappers.py
+++ b/hysop/tools/numpywrappers.py
@@ -1,8 +1,8 @@
 # -*- coding: utf-8 -*-
-"""
-@file numpywrappers.py
+"""Interface to numpy arrays, with hysop predifined types for int, real ...
 
-Tools to build numpy arrays based on hysop setup (float type ...)
+Those functions are useful to enforce hysop predefined types in numpy arrays
+and to keep hysop predefined data layout.
 """
 from hysop.constants import HYSOP_REAL, ORDER, HYSOP_INTEGER,\
     HYSOP_DIM
@@ -31,6 +31,7 @@ def zeros_like(tab):
     """
     return np.zeros_like(tab, dtype=tab.dtype, order=ORDER)
 
+
 def ones_like(tab):
     """
     Wrapper to numpy.ones_like, force order to hysop.constants.ORDER
@@ -42,7 +43,7 @@ def reshape(tab, shape):
     """
     Wrapper to numpy.reshape, force order to hysop.constants.ORDER
     """
-    return np.reshape(tab, shape, dtype=tab.dtype, order=ORDER)
+    return np.reshape(tab, shape, order=ORDER)
 
 
 def realempty(tab):
@@ -160,6 +161,8 @@ def dim_zeros(shape):
 
 
 def equal(a, b):
+    """compare two arrays and ensure data types are the same.
+    """
     msg = 'You try to compare two values of different '
     msg += 'types : ' + str(np.asarray(a).dtype) + 'and '
     msg += str(np.asarray(b).dtype) + '.'
@@ -195,7 +198,7 @@ def add(a, b, c, dtype=HYSOP_REAL):
     return np.add(a, b, c, dtype=dtype)
 
 
-def writeToFile(fname, data, mode='a'):
+def write_to_file(fname, data, mode='a'):
     """
     write data (numpy array) to file fname
     """
@@ -216,3 +219,31 @@ def unlock(tab):
     set tab as a writeable array
     """
     tab.flags.writeable = True
+
+
+def random(shape):
+    """Return an array filled with random values, of type HYSOP_REAL
+
+    Parameters
+    ----------
+    shape : tuple
+        shape of the array
+    """
+    return asrealarray(np.random.random(shape))
+
+
+# Some useful functions, from http://ipython-books.github.io/featured-01/
+def get_data_base(arr):
+    """For a given Numpy array, finds the
+    base array that "owns" the actual data."""
+    base = arr
+    while isinstance(base.base, np.ndarray):
+        base = base.base
+    return base
+
+
+def arrays_share_data(x, y):
+    """Returns true if x and y arrays share the same
+    memory location
+    """
+    return get_data_base(x) is get_data_base(y)
diff --git a/hysop/tools/problem2dot.py b/hysop/tools/problem2dot.py
index abf8d58a7..f0949d8a4 100644
--- a/hysop/tools/problem2dot.py
+++ b/hysop/tools/problem2dot.py
@@ -1,6 +1,4 @@
-"""
-@file problem2dot.py
-Converts a problem instance to a graph throw dot syntax.
+"""Converts a problem instance into a 'dot' graph.
 
 """
 from hysop.operator.advection import Advection
@@ -21,19 +19,29 @@ colors = [
 
 
 def get_shape(op):
+    """Return graph shape depending on the type of mpi communicator"""
     if isinstance(op, Redistribute) or isinstance(op, RedistributeInter):
         return 'octagon'
     else:
         return 'box'
 
 
-def toDot(pb, filename='graph.pdf'):
+def to_dot(pb, filename='graph.pdf'):
+    """Convert problem into graph using dot and
+    save result into a pdf
+
+    Parameters
+    ----------
+    pb : :class:`~hysop.problem.problem.Problem`
+    filename : string, optional
+        output file, default = graph.pdf
+    """
     if main_rank == 0:
         all_ops = []
         all_vars = []
         tasks = []
         for op in pb.operators:
-            if isinstance(op, Advection) and not op.advecDir is None:
+            if isinstance(op, Advection) and op.advecDir is not None:
                 for ad_op in op.advecDir:
                     all_ops.append(ad_op)
                     tasks.append(ad_op.task_id)
@@ -52,8 +60,8 @@ def toDot(pb, filename='graph.pdf'):
                 for req in op.wait_list():
                     if v in req.output:
                         out = req
-                i = op_id-1
-                while(out is None):
+                i = op_id - 1
+                while out is None:
                     if v in all_ops[i].output:
                         if isinstance(all_ops[i], RedistributeInter):
                             if op == all_ops[i].opTo:
@@ -61,7 +69,7 @@ def toDot(pb, filename='graph.pdf'):
                         else:
                             if not isinstance(all_ops[i], Redistribute):
                                 out = all_ops[i]
-                    i = i-1
+                    i = i - 1
                 if (out, op) in all_edges.keys():
                     all_edges[(out, op)].append(v.name)
                 else:
diff --git a/hysop/tools/remeshing_formula_parsing.py b/hysop/tools/remeshing_formula_parsing.py
index cb158d933..279d5767c 100644
--- a/hysop/tools/remeshing_formula_parsing.py
+++ b/hysop/tools/remeshing_formula_parsing.py
@@ -1,10 +1,9 @@
-"""
-@file remeshing_formula_parsing.py
-
-Functions to parse some remeshing formula code (given as strings from Maple
+"""Functions to parse some remeshing formula code (given as strings from Maple
 or Sympy for instance). Result is a formula usable in the
 hysop.numerics.remeshing module or in OpenCL code in
 hysop/gpu/cl_src/remeshing/weights*
+
+To work properly, this module needs sympy python package.
 """
 import re
 try:
@@ -24,14 +23,24 @@ try:
 
     def parse(f, fac=1, vec=False, toOpenCL=True,
               CLBuiltins=False, keep=False):
-        """
-        Parsing function.
-        @param f : functions to parse as string
-        @param fac : numeric factor for all formulas
-        @param vec : is OpenCL output is generated to use vector builtin types
-        @param toOpenCL : is OpenCL output
-        @param CLBuiltins : is OpenCL output uses fma builtin function
-        @param keep : low parsing
+        """Parsing function.
+
+
+        Parameters
+        ----------
+        f : python functions
+            to be parsed as string
+        fac : real
+            numeric factor for all formulas
+        vec : boolean
+            true if we use vector builtin types
+        toOpenCL : boolean
+            true if OpenCL output
+        CLBuiltins : boolean
+            true if OpenCL output uses fma builtin function
+        keep : boolean
+            true low parsing
+
         """
         msg = 'Vector works only in OpenCL parsing'
         assert not (vec and not toOpenCL), msg
@@ -121,7 +130,7 @@ except:
 def createFMA(s):
     """
     Function to handle fma replacements in formula.
-    @param s : formula to parse
+    :param s : formula to parse
 
     \code
     >>> createFMA("(y)")
diff --git a/hysop/tools/tests/test_parameters.py b/hysop/tools/tests/test_io.py
similarity index 100%
rename from hysop/tools/tests/test_parameters.py
rename to hysop/tools/tests/test_io.py
-- 
GitLab