Added CallableNumericalModel (#182)

* Added CallableNumericalModel, which is similar to CallableModel but whose dict components are python callables rather than sympy Expr.

* Added a test for constraints in conjunction with CallableNumericalModel

* Added Powell minimizer, and a small bugfix for this method.

* __neg__ dind't negate yet, fixed.

* Cleaned up inheritance structure to make the changes minimal, added docstrings.

* Added tests for multiple components.

* Fixed a problem with multidimensional data.

Was caused by a small bug in finiti differences, which assumed 1D data.

* Made the fix < py3.5 compatible

* Updated tests to include a mixed model and wrong argument name

* Enabled cashing on models whose components are mixed between pythonic and sympy.

* Also check that 2D is the same as 1D after flattening.

* Added example to the docs for CallableNumericalModel

docs/examples/ex_CallableNumericalModel.rst

@@ -0,0 +1,20 @@
+Example: CallableNumericalModel
+===============================
+
+Below is an example of how to use the
+:class:`symfit.core.fit.CallableNumericalModel`. This class allows you to
+provide custom callables as your model, while still allowing clean interfacing
+with the :mod:`symfit` API.
+
+These models also accept a mixture of symbolic and callable components, as will
+be demonstrated below. This allows the power-user great flexibility, since it is
+still easy to interface with :mod:`symfit`'s constraints, minimizers, etc.
+
+.. literalinclude:: ../../examples/callable_numerical_model.py
+    :language: python
+
+This is the resulting fit:
+
+.. figure:: ../_static/callable_numerical_model.png
+   :width: 500px
+   :alt: Custom Callable Model

docs/examples/index.rst

@@ -17,6 +17,7 @@ paper.
     ex_piecewise
     ex_poly_surface_fit
     ex_ODEModel
+    ex_CallableNumericalModel
 
 Interactive Guess Module
 ------------------------

examples/callable_numerical_model.py

@@ -0,0 +1,45 @@
+from symfit import variables, parameters, Fit, D, ODEModel, CallableNumericalModel
+import numpy as np
+import matplotlib.pyplot as plt
+
+def nonanalytical_func(x, a, b):
+    """
+    This can be any pythonic function which should be fitted, typically one
+    which is not easily written or supported as an analytical expression.
+    """
+    # Do your non-trivial magic here. In this case a Piecewise, although this
+    # could also be done symbolically.
+    y = np.zeros_like(x)
+    y[x > b] = (a * (x - b) + b)[x > b]
+    y[x <= b] = b
+    return y
+
+x, y1, y2 = variables('x, y1, y2')
+a, b = parameters('a, b')
+
+mixed_model = CallableNumericalModel(
+    {y1: nonanalytical_func, y2: x ** a},
+    independent_vars=[x],
+    params=[a, b]
+)
+
+# Generate data
+xdata = np.linspace(0, 10)
+y1data, y2data = mixed_model(x=xdata, a=1.3, b=4)
+y1data = np.random.normal(y1data, 0.1 * y1data)
+y2data = np.random.normal(y2data, 0.1 * y2data)
+
+# Perform the fit
+b.value = 3.5
+fit = Fit(mixed_model, x=xdata, y1=y1data, y2=y2data)
+fit_result = fit.execute()
+print(fit_result)
+
+# Plotting, irrelevant to the symfit part.
+y1_fit, y2_fit, = mixed_model(x=xdata, **fit_result.params)
+plt.scatter(xdata, y1data)
+plt.plot(xdata, y1_fit, label=r'$y_1$')
+plt.scatter(xdata, y2data)
+plt.plot(xdata, y2_fit, label=r'$y_2$')
+plt.legend(loc=0)
+plt.show()

symfit/api.py

@@ -5,7 +5,7 @@
 from symfit.core.fit import (
     Fit, Model, Constraint,
     LinearLeastSquares, NonLinearLeastSquares,
-    TaylorModel, ODEModel, ModelError
+    TaylorModel, ODEModel, ModelError, CallableModel, CallableNumericalModel
 )
 from symfit.core.fit_results import FitResults
 from symfit.core.argument import Variable, Parameter

symfit/core/fit.py

@@ -35,10 +35,12 @@ class ModelError(Exception):
     """
     pass
 
+
 class BaseModel(Mapping):
     """
     ABC for ``Model``'s. Makes sure models are iterable.
-    Models can be initiated from Mappings or Iterables of Expressions, or from an expression directly.
+    Models can be initiated from Mappings or Iterables of Expressions,
+    or from an expression directly.
     Expressions are not enforced for ducktyping purposes.
     """
     def __init__(self, model):
@@ -191,38 +193,109 @@ def shared_parameters(self):
             else:
                 return False
 
-    @abstractmethod
     def __str__(self):
         """
-        A representation upon printing has to provided.
+        Printable representation of a Mapping model.
+
+        :return: str
         """
+        template = "{}({}; {}) = {}"
+        parts = []
+        for var, expr in self.items():
+            parts.append(template.format(
+                    var,
+                    ", ".join(arg.name for arg in self.independent_vars),
+                    ", ".join(arg.name for arg in self.params),
+                    expr
+                )
+            )
+        return "\n".join(parts)
 
 
-class CallableModel(BaseModel):
+class BaseNumericalModel(BaseModel):
     """
-    Defines a callable model. The usual rules apply to the ordering of the arguments:
+    ABC for Numerical Models. These are models whose components are generic
+    python callables.
+    """
+    def __init__(self, model, independent_vars, params):
+        """
 
-    * first independent variables, then dependent variables, then parameters.
-    * within each of these groups they are ordered alphabetically.
+        :param model: dict of ``callable``, where dependent variables are the
+            keys. If instead of a dict a (sequence of) ``callable`` is provided,
+            it will be turned into a dict automatically.
+        :param independent_vars: The independent variables of the model.
+        :param params: The parameters of the model.
+        """
+        self.independent_vars = sorted(independent_vars, key=str)
+        self.params = sorted(params, key=str)
+        super(BaseNumericalModel, self).__init__(model)
+
+    def _init_from_dict(self, model_dict):
+        """
+        Initiate self from a model_dict which has python callables as its values.
+        This init makes sure self.dependent_vars is available, the other are
+        set in __init__.
+
+        :param model_dict: dict of (dependent_var, callable) pairs.
+        """
+        sort_func = lambda symbol: str(symbol)
+        self.model_dict = OrderedDict(sorted(model_dict.items(), key=lambda i: sort_func(i[0])))
+        self.dependent_vars = sorted(model_dict.keys(), key=sort_func)
+
+        # Make Variable object corresponding to each var.
+        self.sigmas = {var: Variable(name='sigma_{}'.format(var.name)) for var in self.dependent_vars}
+
+    def __eq__(self, other):
+        raise NotImplementedError(
+            'Equality checking for {} is ambiguous.'.format(self.__class__.__name__)
+        )
+
+    def __neg__(self):
+        """
+        :return: new model with opposite sign. Does not change the model in-place,
+            but returns a new copy.
+        """
+        new_model_dict = {}
+        for key, callable_expr in self.model_dict.values():
+            new_model_dict[key] = lambda *args, **kwargs: - callable_expr(*args, **kwargs)
+        return self.__class__(new_model_dict)
+
+    @property
+    def shared_parameters(self):
+        raise NotImplementedError(
+            'Shared parameters can not be inferred for {}'.format(self.__class__.__name__)
+        )
+
+
+class BaseCallableModel(BaseModel):
+    """
+    Baseclass for callable models.
     """
-    @abstractmethod
     def eval_components(self, *args, **kwargs):
         """
-        Evaluate the components of the model with the given data.
-        Used for numerical evaluation.
+        :return: lambda functions of each of the components in model_dict, to be used in numerical calculation.
         """
+        return [expr(*args, **kwargs) for expr in self.numerical_components]
 
+    @abstractmethod
+    def numerical_components(self):
+        """
+        :return: A list of callables corresponding to each of the components
+            of the model.
+        """
 
     def _make_signature(self):
         # Handle args and kwargs according to the allowed names.
-        parameters = [  # Note that these are inspect_sig.Parameter's, not symfit parameters!
-            inspect_sig.Parameter(arg.name, inspect_sig.Parameter.POSITIONAL_OR_KEYWORD)
-                for arg in self.independent_vars + self.params
+        parameters = [
+            # Note that these are inspect_sig.Parameter's, not symfit parameters!
+            inspect_sig.Parameter(arg.name,
+                                  inspect_sig.Parameter.POSITIONAL_OR_KEYWORD)
+            for arg in self.independent_vars + self.params
         ]
         return inspect_sig.Signature(parameters=parameters)
 
     def _init_from_dict(self, model_dict):
-        super(CallableModel, self)._init_from_dict(model_dict)
+        super(BaseCallableModel, self)._init_from_dict(model_dict)
         self.__signature__ = self._make_signature()
 
     def __call__(self, *args, **kwargs):
@@ -240,16 +313,8 @@ def __call__(self, *args, **kwargs):
         """
         bound_arguments = self.__signature__.bind(*args, **kwargs)
         Ans = namedtuple('Ans', [var.name for var in self])
-        # return Ans(*[component(**bound_arguments.arguments) for component in self.numerical_components])
         return Ans(*self.eval_components(**bound_arguments.arguments))
 
-    @cached_property
-    def numerical_components(self):
-        """
-        :return: lambda functions of each of the components in model_dict, to be used in numerical calculation.
-        """
-        return [sympy_to_py(expr, self.independent_vars, self.params) for expr in self.values()]
-
     @keywordonly(dx=1e-8)
     def finite_difference(self, *args, **kwargs):
         """
@@ -302,15 +367,17 @@ def finite_difference(self, *args, **kwargs):
                     # a list of arrays.
                     for comp_idx in range(len(self)):
                         try:
-                            data_len = len(up[comp_idx])
+                            len(up[comp_idx])
                         except TypeError:  # output[comp_idx] is a number
-                            data_len = 1
+                            data_shape = (1,)
+                        else:
+                            data_shape = up[comp_idx].shape
                         # Initialize at 0 so we can += all the contributions
-                        param_grad = np.zeros((len(self.params), data_len), dtype=float)
+                        param_grad = np.zeros([len(self.params)] + list(data_shape), dtype=float)
                         out.append(param_grad)
                 for comp_idx in range(len(self)):
                     diff = up[comp_idx] - down[comp_idx]
-                    out[comp_idx][param_idx, :] += factor * diff/(2*h[param_idx])
+                    out[comp_idx][param_idx, :] += factor * diff / (2 * h[param_idx])
         return out
 
     def eval_jacobian(self, *args, **kwargs):
@@ -320,6 +387,58 @@ def eval_jacobian(self, *args, **kwargs):
         return self.finite_difference(*args, **kwargs)
 
 
+class CallableNumericalModel(BaseCallableModel, BaseNumericalModel):
+    """
+    Callable model, whose components are callables provided by the user.
+    This allows the user to provide the components directly.
+
+    Example::
+
+        x, y = variables('x, y')
+        a, b = parameters('a, b')
+        numerical_model = CallableNumericalModel(
+            {y: lambda x, a, b: a * x + b},
+            independent_vars=[x],
+            params=[a, b]
+        )
+
+    This is identical in functionality to the more traditional::
+
+        x, y = variables('x, y')
+        a, b = parameters('a, b')
+        model = CallableModel({y: a * x + b})
+
+    but allows power-users a lot more freedom while still interacting
+    seamlessly with the :mod:`symfit` API.
+
+    .. note:: All of the callables must accept all of the ``independent_vars``
+        and  ``params`` of the model as arguments, even if not all of them are
+        used by every callable.
+    """
+    @cached_property
+    def numerical_components(self):
+        return [expr if not isinstance(expr, sympy.Expr) else
+                sympy_to_py(expr, self.independent_vars, self.params)
+                for expr in self.values()]
+
+
+class CallableModel(BaseCallableModel):
+    """
+    Defines a callable model. The usual rules apply to the ordering of the
+    arguments:
+
+    * first independent variables, then dependent variables, then parameters.
+    * within each of these groups they are ordered alphabetically.
+    """
+    @cached_property
+    def numerical_components(self):
+        """
+        :return: lambda functions of each of the analytical components in
+        model_dict, to be used in numerical calculation.
+        """
+        return [sympy_to_py(expr, self.independent_vars, self.params) for expr in self.values()]
+
+
 class Model(CallableModel):
     """
     Model represents a symbolic function and all it's derived properties such as sum of squares, jacobian etc.
@@ -336,24 +455,6 @@ class Model(CallableModel):
     Models are also iterable, behaving as their internal model_dict. In the example above,
     a[y] returns x**2, len(a) == 1, y in a == True, etc.
     """
-    def __str__(self):
-        """
-        Printable representation of this model.
-
-        :return: str
-        """
-        template = "{}({}; {}) = {}"
-        parts = []
-        for var, expr in self.items():
-            parts.append(template.format(
-                    var,
-                    ", ".join(arg.name for arg in self.independent_vars),
-                    ", ".join(arg.name for arg in self.params),
-                    expr
-                )
-            )
-        return "\n".join(parts)
-
     @property
     def jacobian(self):
         """
@@ -490,14 +591,6 @@ def eval_jacobian(self, *args, **kwargs):
             jac[idx] = np.array(jac[idx], dtype=float)
         return jac
 
-    def eval_components(self, *args, **kwargs):
-        """
-        :return: lambda functions of each of the components in model_dict, to be used in numerical calculation.
-        """
-        return [expr(*args, **kwargs) for expr in self.numerical_components]
-        # return [sympy_to_py(expr, self.independent_vars, self.params)(*args, **kwargs) for expr in self.values()]
-
-
 
 class TaylorModel(Model):
     """

symfit/core/minimizers.py

@@ -293,7 +293,7 @@ def _pack_output(self, ans):
         }
 
         best_vals = []
-        found = iter(ans.x)
+        found = iter(np.atleast_1d(ans.x))
         for param in self.parameters:
             if param.fixed:
                 best_vals.append(param.value)
@@ -461,6 +461,12 @@ def method_name(cls):
         return 'Nelder-Mead'
 
 
+class Powell(ScipyMinimize, BaseMinimizer):
+    """
+    Wrapper around :func:`scipy.optimize.minimize`'s Powell algorithm.
+    """
+
+
 class BasinHopping(ScipyMinimize, BaseMinimizer):
     """
     Wrapper around :func:`scipy.optimize.basinhopping`'s basin-hopping algorithm.