Enh/20 new potential pipeline model

[sannant]

Simplified Potential Implementations

This is a draft for the moment.

The support for arrays of work of adhesion needs passing of masks and complicates the use of "evaluate" and naive pot.

One source of complication is that we support the distances that are scalars or arrays of any shape, and this clashes with the use of a mask.

I suggest that we now only support two dimensional arrays in the functions that are called during a simulation.

For other purposes, like determining the parameters of a derived potential, we can provide a more flexible convenience interface for the energy, the gradient and the curvature (each separately) since we don't care about performance in these cases. This flexible convenience interface can be implemented once for all in the parent class.

[pastewka]

Hi @antoinesimtek, do you want me to look at this?

[sannant]

I wanted to make a draft pull request, it's far from ready. Maybe the best is that we shortly discuss this live today.

[pastewka]

Can you prefix draft pull requests with WIP?

[pastewka]

I just did this, remove the WIP prefix once this is ready

[sannant]

Here a proposal how the simplified class structure could look like.

evaluate is a wrapper to _evaluate implementing area_scale. The only reason evaluate is split in evaluate and _evaluate is that it avoids to need to implement handling the area_scale thing in every potential.

We could probably get rid of it if the user rescales the work of adhesion according to the pixel size or if we implement a rescale pipeline potential . See #12

class Potential(SoftWall, metaclass=abc.ABCMeta):
    """ Describes the minimum interface to interaction potentials for
        Adhesion. These potentials are purely 1D, which allows for a few
        simplifications. For instance, the potential energy and forces can be
        computed at any point in the problem from just the one-dimensional gap
        (h(x,y)-z(x,y)) at that point
    """
    _functions={}
    name = "generic_potential"

    class PotentialError(Exception):
        "umbrella exception for potential-related issues"
        pass

    class SliceableNone(object):
        """small helper class to remedy numpy's lack of views on
        index-sliced array views. This construction avoid the computation
        of all interactions as with np.where, and copies"""
        # pylint: disable=too-few-public-methods
        __slots__ = ()

        def __setitem__(self, index, val):
            pass

        def __getitem__(self, index):
            pass

    @classmethod
    def register_function(cls, name, function):
        cls._functions.update({name: function})

    def __getattr__(self, name):
        if name in self._functions:
            def func(*args, **kwargs):
                return self._functions[name](self, *args, **kwargs)

            func.__doc__ = self._functions[name].__doc__
            return func
        else:
            raise AttributeError(
                "Unkown attribute '{}' and no analysis or pipeline function "
                "of this name registered (class {}). Available functions: {}"
                .format(name, self.__class__.__name__,
                        ', '.join(self._functions.keys())))

    def __dir__(self):
        return sorted(super().__dir__() + [*self._functions])

    @abc.abstractmethod
    def __init__(self, communicator=MPI.COMM_WORLD):
        super().__init__(communicator)
        self.curvature = None

    @abc.abstractmethod
    def __repr__(self):
        return ("Potential '{0.name}'").format(self)

    @property
    def has_cutoff(self):
        return False

    def compute(self, gap, potential=True, gradient=False, curvature=False,
                area_scale=1.):
        """

        Updates the values of self.energy, self.gradient, self.curvature
        according to the provided gap, by calling evaluate()

        Parameters
        ----------
        gap:
            array of distances between the two surfaces
        potential: bool (default True)
            if true, returns potential energy
        gradient: bool, (default False)
            if true, returns gradient
        curvature: bool, (default False)
            if true, returns second derivative
        area_scale: float (default 1.)
            scale by this.
            (Interaction quantities are supposed to be expressed per unit
            area, so systems need to be able to scale their response for their
            nb_grid_pts)
        Returns
        -------
            Nothing

        """
        energy, self.gradient, self.curvature = self.evaluate(
            gap, potential, gradient, curvature, area_scale=area_scale)
        self.energy = self.pnp.sum(energy) if potential else None

    def evaluate(self, gap, potential=True, gradient=False, curvature=False,
                 area_scale=1.):
        """Evaluates the potential and its derivatives

        and scales the values by area_scale.

        Parameters:
        -----------
        gap:
            array of distances between the two surfaces
        potential: bool (default True)
            if true, returns potential energy
        gradient: bool, (default False)
            if true, returns gradient
        curvature: bool, (default False)
            if true, returns second derivative
        area_scale: float (default 1.)
            scale by this.
            (Interaction quantities are supposed to be expressed per unit
            area, so systems need to be able to scale their response for their
            nb_grid_pts)
        """

        # This is a wrapper to _evaluate implementing area_scale.
        # The only reason evaluate is split in evaluate and _evaluate is that
        # it avoids to need to implement handling the area_scale thing in
        # every potential.

        if np.isscalar(gap): # TODO: eventually cancel support for scalars
            gap = np.asarray(gap)
        if gap.shape == ():
            gap.shape = (1,)

        V, dV, ddV = self._evaluate(gap, potential, gradient, curvature)

        return (area_scale * V if potential else None,
                area_scale * dV if gradient else None,
                area_scale * ddV if curvature else None)

    @abc.abstractmethod
    def _evaluate(self, r, potential=True, gradient=False, curvature=False,
                  mask=None):
        """Evaluates the potential and its derivatives
        Keyword Arguments:
        r          -- array of distances
        pot        -- (default True) if true, returns potential energy
        gradient   -- (default False) if true, returns gradient
        curvature  -- (default False) if true, returns second derivative
        """
        if mask is None:
            mask = (slice(None),) * len(r.shape)
        return self.naive_pot(r, potential, gradient, curvature, mask=mask)

    @abc.abstractproperty
    def r_min(self):
        """
        convenience function returning the location of the energy minimum
        """
        raise NotImplementedError()

    @abc.abstractproperty
    def r_infl(self):
        """
        convenience function returning the location of the potential's
        inflection point (if applicable)
        """
        raise NotImplementedError()

    @property
    def max_tensile(self):
        """
        convenience function returning the value of the maximum stress
        (at r_infl)
        """
        max_tensile = self.evaluate(self.r_infl, gradient=True)[1]
        return max_tensile.item() if np.prod(max_tensile.shape) == 1 \
            else max_tensile

    @property
    def v_min(self):
        """ convenience function returning the value of the energy minimum
        """
        return float(self.evaluate(self.r_min)[0])

    def ancestor_potential(self):
        return self

class ChildPotential(Potential):
    def __init__(self, parent_potential):
        self.parent_potential = parent_potential
        self.pnp = parent_potential.pnp
        self.communicator = parent_potential.communicator

    def __getattr__(self, item):
        if item[:2] == "__" and item[-2:] == "__":
            raise AttributeError
        else:
            return getattr(self.parent_potential, item)

    def __getstate__(self):
        state = super().__getstate__(), self.parent_potential
        return state

    def __setstate__(self, state):
        superstate, self.parent_potential = state
        super().__setstate__(superstate)

    def ancestor_potential(self):
        return self.parent_potential.ancestor_potential()

[sannant]

Update of the proposal

We discussed in #12 that area_scale will be handled by the system, so we will only need the function evaluate instead of _evaluate and evaluate.

Also, we should probably move the storage of the last evaluated energy and force to the system, making the function compute obsolete. If we do so, I think there will not be any parallelisation in these interaction potentials anymore and we can get rid of the communicator.

@pastewka, what is your opinion ?

class SoftWall(Interaction):
    """base class for smooth contact mechanics"""
    def __init__(self):
        pass

    def __deepcopy__(self, memo):
        """
        makes a deepcopy of all the attributes except self.pnp, where it stores the same reference

        Parameters
        ----------
        memo

        Returns
        -------
        result SoftWall instance
        """

        result = self.__class__.__new__(self.__class__)
        memo[id(self)] = result

        keys = set(self.__dict__.keys())

        for k in keys:
            setattr(result, k, copy.deepcopy(getattr(self, k), memo))
        return result

    def __getstate__(self):
        return self.energy, self.gradient

    def __setstate__(self, state):
        self.energy, self.gradient = state

    def evaluate(self, gap, potential=True, gradient=False, area_scale=1.):
        """
        computes and returns the interaction energy and/or forces based on the
        as fuction of the gap
        Parameters:
        gap        -- array containing the point-wise gap values
        potential        -- (default True) whether the energy should be evaluated
        forces     -- (default False) whether the forces should be evaluated
        area_scale -- (default 1.) scale by this. (Interaction quantities are
                      supposed to be expressed per unit area, so systems need
                      to be able to scale their response for their nb_grid_pts))
        """
        raise NotImplementedError()

    @property
    def force(self): # This is more a legacy
        if self.gradient is not None:
            return - self.gradient
        return None

class Potential(SoftWall, metaclass=abc.ABCMeta):
    """ Describes the minimum interface to interaction potentials for
        Adhesion. These potentials are purely 1D, which allows for a few
        simplifications. For instance, the potential energy and forces can be
        computed at any point in the problem from just the one-dimensional gap
        (h(x,y)-z(x,y)) at that point
    """
    _functions={}
    name = "generic_potential"

    class PotentialError(Exception):
        "umbrella exception for potential-related issues"
        pass

    class SliceableNone(object):
        """small helper class to remedy numpy's lack of views on
        index-sliced array views. This construction avoid the computation
        of all interactions as with np.where, and copies"""
        # pylint: disable=too-few-public-methods
        __slots__ = ()

        def __setitem__(self, index, val):
            pass

        def __getitem__(self, index):
            pass

    @classmethod
    def register_function(cls, name, function):
        cls._functions.update({name: function})

    def __getattr__(self, name):
        if name in self._functions:
            def func(*args, **kwargs):
                return self._functions[name](self, *args, **kwargs)

            func.__doc__ = self._functions[name].__doc__
            return func
        else:
            raise AttributeError(
                "Unkown attribute '{}' and no analysis or pipeline function "
                "of this name registered (class {}). Available functions: {}"
                .format(name, self.__class__.__name__,
                        ', '.join(self._functions.keys())))

    def __dir__(self):
        return sorted(super().__dir__() + [*self._functions])

    @abc.abstractmethod
    def __init__(self, ):
         pass

    @abc.abstractmethod
    def __repr__(self):
        return ("Potential '{0.name}'").format(self)

    @property
    def has_cutoff(self):
        return False

    def evaluate(self, gap, potential=True, gradient=False, curvature=False,
                 area_scale=1., mask=None):
        """Evaluates the potential and its derivatives

        and scales the values by area_scale.

        Parameters:
        -----------
        gap:
            array of distances between the two surfaces
        potential: bool (default True)
            if true, returns potential energy
        gradient: bool, (default False)
            if true, returns gradient
        curvature: bool, (default False)
            if true, returns second derivative
        area_scale: float (default 1.)
            scale by this.
            (Interaction quantities are supposed to be expressed per unit
            area, so systems need to be able to scale their response for their
            nb_grid_pts)
        """
        raise NotImplementedError()

    @abc.abstractproperty
    def r_min(self):
        """
        convenience function returning the location of the energy minimum
        """
        raise NotImplementedError()

    @abc.abstractproperty
    def r_infl(self):
        """
        convenience function returning the location of the potential's
        inflection point (if applicable)
        """
        raise NotImplementedError()

    @property
    def max_tensile(self):
        """
        convenience function returning the value of the maximum stress
        (at r_infl)
        """
        max_tensile = self.evaluate(self.r_infl, gradient=True)[1]
        return max_tensile.item() if np.prod(max_tensile.shape) == 1 \
            else max_tensile

    @property
    def v_min(self):
        """ convenience function returning the value of the energy minimum
        """
        return float(self.evaluate(self.r_min)[0])

    def ancestor_potential(self):
        return self

class ChildPotential(Potential):
    def __init__(self, parent_potential):
        self.parent_potential = parent_potential

    def __getattr__(self, item):
        if item[:2] == "__" and item[-2:] == "__":
            raise AttributeError
        else:
            return getattr(self.parent_potential, item)

    def __getstate__(self):
        state = super().__getstate__(), self.parent_potential
        return state

    def __setstate__(self, state):
        superstate, self.parent_potential = state
        super().__setstate__(superstate)

    def ancestor_potential(self):
        return self.parent_potential.ancestor_potential()

[pastewka]

You proposal sounds good. We should not get rid of the communicator entirely, because we could envision solving the Poisson/Laplace equation for the electrostatic problem which would require parallelization. There may be other more complication interaction potentials that also require this. I agree that all of the present interaction models don't require communication.

[sannant]

I suggest to remove the docstring of inheritingclass(Potential).evaluate methods, so that the docstring of the parent class "Potential" is inherited. That way we don't have to write it multiple times in an identical way.

What do you think about this @pastewka

[pastewka]

Are you sure the docstring is automatically inherited. I don't think so.

[sannant]

I tried in the ipython console:

Signature: pot.evaluate(r, potential=True, gradient=False, curvature=False, mask=None)
Docstring:
Evaluates the potential and its derivatives

Parameters:
-----------
gap: array_like
    array of distances between the two surfaces
potential: bool (default True)
    if true, returns potential energy
gradient: bool, (default False)
    if true, returns gradient
curvature: bool, (default False)
    if true, returns second derivative
mask: boolean array, optional
    potential is only evaluated on gap[mask]
    this property is used by the child potential
Source:   
    def evaluate(self, r, potential=True, gradient=False, curvature=False, mask=None):
        r = np.asarray(r)

        V = dV = ddV = None
        sig_r3 = (self.sig/r)**3
        sig_r9 = sig_r3**3

        if potential:
            V = self.eps*(2./15*sig_r9 - sig_r3)
        if gradient or curvature:
            eps_r = self.eps/r
        if gradient:
            dV = - eps_r*(6./5*sig_r9 - 3*sig_r3)
        if curvature:
            ddV = 12*eps_r/r*(sig_r9 - sig_r3)

[sannant]

But it could be that other things like sphinx don't understand this, since pot.evaluate.__doc__ is empty

Note that what I suggest means that we have to move the formulas from evaluate to the class docstring (I think this is fine). example for formulas in evaluate:

Signature:
pot.evaluate(
    gap,
    potential=True,
    gradient=False,
    curvature=False,
    mask=(slice(None, None, None), slice(None, None, None)),
)
Source:   
    def evaluate(self, gap, potential=True, gradient=False, curvature=False,
                 mask=(slice(None), slice(None))):
        """ Evaluates the potential and its derivatives.
            These have been collected in a single method to reuse the
            computed vdW terms for efficiency

                           A      C_sr
            vdW(r)  = - ─────── + ────
                            2       8
                        12⋅r ⋅π    r

                        A      8⋅C_sr
            vdW'(r) = ────── - ──────
                         3        9
                      6⋅r ⋅π     r

                          A      72⋅C_sr
            vdW"(r) = - ────── + ───────
                           4        10
                        2⋅r ⋅π     r
            Parameters:
            -------------
            gap: array of distances
            potential    -- (default True) if true, returns potential energy
            gradient -- (default False) if true, returns gradient
            curvature   -- (default False) if true, returns second derivative
:

[pastewka]

Do these formulas show in Sphinx? You can also use some sort of LaTeX notation for formulas in Sphinx.

[sannant]

Yes this is on the todolist to convert them to LateX

[sannant]

The inherited docstring also shows up in sphinx

[sannant]

The ASCII-Art Formulas are ugly.

[sannant]

I removed the evaluate docstrings of the child Potentials.

[sannant]

@pastewka , I now replaced the use of ancestor potential by pot.pipeline(), in analogy to what is done in SurfaceTopography.

With this, I think this branch is merge ready.