Repository: tBuLi/symfit

Issue:
Global LogLikelihood fails
In principle LogLikelihood fitting should also work when fitting to multiple data sets, however this [SO ](https://stackoverflow.com/questions/56006357/global-distribution-fitting-sharing-some-parameters-without-any-specification-of)question shows that it doesn't.

The issue seems to be that numpy as not able to take the log if the model components have different sizes. We should therefore take this log per component before summing.

Discussion:

Patch:
diff --git a/docs/examples/ex_multidataset_likelihood.rst b/docs/examples/ex_multidataset_likelihood.rst
new file mode 100644
index 0000000..9cd8b2d
--- /dev/null
+++ b/docs/examples/ex_multidataset_likelihood.rst
@@ -0,0 +1,13 @@
+Example: Global Likelihood fitting
+==================================
+
+In this example, we shall perform likelihood fitting to two data sets at the
+same time to show how likelihood fitting can be used to fit multiple data sets
+at the same time.
+
+First we will fit using a different parameter for each data set, and then with
+the same for both data sets. This shows the ease with which we can experiment
+with different models.
+
+.. literalinclude:: ../../examples/multidataset_likelihood.py
+    :language: python
diff --git a/docs/examples/index.rst b/docs/examples/index.rst
index b5dd3d8..abd5852 100644
--- a/docs/examples/index.rst
+++ b/docs/examples/index.rst
@@ -29,6 +29,7 @@ standard least-squares objective, such as likelihood fitting.
     :maxdepth: 1
 
     ex_bivariate_likelihood
+    ex_multidataset_likelihood
 
 Interactive Guess Module
 ------------------------
diff --git a/docs/fitting_types.rst b/docs/fitting_types.rst
index c76f9f5..3adfeeb 100644
--- a/docs/fitting_types.rst
+++ b/docs/fitting_types.rst
@@ -115,6 +115,10 @@ Example::
 ``fit_result`` is a normal :class:`~symfit.core.fit_results.FitResults` object.
 As always, bounds on parameters and even constraints are supported.
 
+Multiple data sets can be likelihood fitted simultaneously by merging this
+example with that of global fitting, see `Example: Global Likelihood fitting`
+in the example section.
+
 .. _minimize_maximize:
 
 Minimize/Maximize
diff --git a/docs/module_docs.rst b/docs/module_docs.rst
index 8228e4b..f08b782 100644
--- a/docs/module_docs.rst
+++ b/docs/module_docs.rst
@@ -68,6 +68,15 @@ Support
    :exclude-members: __weakref__
    :show-inheritance:
 
+Printing
+--------
+
+.. automodule:: symfit.core.printing
+   :members:
+   :special-members:
+   :exclude-members: __weakref__
+   :show-inheritance:
+
 Distributions
 -------------
 
diff --git a/examples/multidataset_likelihood.py b/examples/multidataset_likelihood.py
new file mode 100644
index 0000000..d56e91c
--- /dev/null
+++ b/examples/multidataset_likelihood.py
@@ -0,0 +1,31 @@
+import numpy as np
+from symfit import variables, parameters, Fit, exp, Model
+from symfit.core.objectives import LogLikelihood
+
+# Draw samples from a bivariate distribution
+np.random.seed(42)
+data1 = np.random.exponential(5.5, 1000)
+data2 = np.random.exponential(6, 2000)
+
+# Define the model for an exponential distribution (numpy style)
+a, b = parameters('a, b')
+x1, y1, x2, y2 = variables('x1, y1, x2, y2')
+model = Model({
+    y1: (1 / a) * exp(-x1 / a),
+    y2: (1 / b) * exp(-x2 / b)
+})
+print(model)
+
+fit = Fit(model, x1=data1, x2=data2, objective=LogLikelihood)
+fit_result = fit.execute()
+print(fit_result)
+
+# Instead, we could also fit with only one parameter to see which works best
+model = Model({
+    y1: (1 / a) * exp(-x1 / a),
+    y2: (1 / a) * exp(-x2 / a)
+})
+
+fit = Fit(model, x1=data1, x2=data2, objective=LogLikelihood)
+fit_result = fit.execute()
+print(fit_result)
\ No newline at end of file
diff --git a/symfit/core/fit.py b/symfit/core/fit.py
index d79621c..c713e54 100644
--- a/symfit/core/fit.py
+++ b/symfit/core/fit.py
@@ -1095,6 +1095,34 @@ class TakesData(object):
         return np.array([param.value for param in self.model.params])
 
 
+class BaseFit(TakesData):
+    """
+    Abstract base class for all fitting objects.
+    """
+    def execute(self, *args, **kwargs):
+        """
+        Every fit object has to define an execute method.
+        Any * and ** arguments will be passed to the fitting module that is being wrapped, e.g. leastsq.
+
+        :args kwargs:
+        :return: Instance of FitResults
+        """
+        raise NotImplementedError('Every subclass of BaseFit must have an execute method.')
+
+    def error_func(self, *args, **kwargs):
+        """
+        Every fit object has to define an error_func method, giving the function to be minimized.
+        """
+        raise NotImplementedError('Every subclass of BaseFit must have an error_func method.')
+
+    def eval_jacobian(self, *args, **kwargs):
+        """
+        Every fit object has to define an eval_jacobian method, giving the jacobian of the
+        function to be minimized.
+        """
+        raise NotImplementedError('Every subclass of BaseFit must have an eval_jacobian method.')
+
+
 class HasCovarianceMatrix(TakesData):
     """
     Mixin class for calculating the covariance matrix for any model that has a
diff --git a/symfit/core/objectives.py b/symfit/core/objectives.py
index 584defd..76601e3 100644
--- a/symfit/core/objectives.py
+++ b/symfit/core/objectives.py
@@ -381,7 +381,9 @@ class LogLikelihood(HessianObjective):
             ordered_parameters, **parameters
         )
 
-        ans = - np.nansum(np.log(evaluated_func))
+        ans = - np.nansum(
+            [np.nansum(np.log(component)) for component in evaluated_func]
+        )
         return ans
 
     @keywordonly(apply_func=np.nansum)
@@ -403,15 +405,17 @@ class LogLikelihood(HessianObjective):
         )
 
         result = []
-        for jac_comp in evaluated_jac:
+        for component, jac_comp in zip(evaluated_func, evaluated_jac):
+            component_sums = []
             for df in jac_comp:
-                result.append(
+                component_sums.append(
                     - apply_func(
-                        df.flatten() / evaluated_func
+                        df / component
                     )
                 )
-        else:
-            return np.atleast_1d(np.squeeze(np.array(result)))
+            result.append(component_sums)
+        result = np.sum(result, axis=0)
+        return np.atleast_1d(np.squeeze(np.array(result)))
 
     def eval_hessian(self, ordered_parameters=[], **parameters):
         """
@@ -440,8 +444,8 @@ class LogLikelihood(HessianObjective):
             # We sum away everything except the matrices in the axes 0 & 1.
             axes = tuple(range(2, len(dd_logf.shape)))
             result += np.sum(dd_logf, axis=axes, keepdims=False)
-        else:
-            return np.atleast_2d(np.squeeze(np.array(result)))
+
+        return np.atleast_2d(np.squeeze(np.array(result)))
 
 
 class MinimizeModel(HessianObjective):
diff --git a/symfit/core/printing.py b/symfit/core/printing.py
index 3f95fd4..45dfea2 100644
--- a/symfit/core/printing.py
+++ b/symfit/core/printing.py
@@ -1,9 +1,18 @@
-from sympy.printing.pycode import NumPyPrinter
+"""
+``symfit`` occasionally updates the printing of ``sympy`` objects, such that they
+print into their ``numpy``/``scipy`` equivalent. This is done because sometimes
+such printing has not been implemented in ``sympy`` yet, or because we want
+slightly different behavior from the standard one.
+
+Users using both ``symfit`` and ``sympy`` should be aware of this.
+"""
+
+from sympy import HadamardProduct, MatMul, MatPow, Idx, Inverse
 from sympy.printing.codeprinter import CodePrinter
 
-#########################################################
-# Monkeypatch the printer until this is merged upstream #
-#########################################################
+##########################################################
+# Monkeypatch the printers until this is merged upstream #
+##########################################################
 
 class DontDeleteMe(object):
     def __init__(self, default_value):
@@ -22,61 +31,50 @@ class DontDeleteMe(object):
 CodePrinter._not_supported = DontDeleteMe(set())
 CodePrinter._number_symbols = DontDeleteMe(set())
 
-#########################################################
-# End of monkeypatch                                    #
-#########################################################
-
-class SymfitNumPyPrinter(NumPyPrinter):
+def _print_MatMul(self, printer):
     """
-    Our own NumpyPrinter subclass, in case we need to print certain numpy
-    features which are not yet supported in SymPy.
+    Matrix multiplication printer. The sympy one turns everything into a
+    dot product without type-checking.
     """
-    def _print_MatMul(self, expr):
-        """
-        Matrix multiplication printer. The sympy one turns everything into a
-        dot product without type-checking.
-        """
-        from sympy import MatrixExpr
-        links = []
-        for i, j in zip(expr.args[1:], expr.args[:-1]):
-            if isinstance(i, MatrixExpr) and isinstance(j, MatrixExpr):
-                links.append(').dot(')
-            else:
-                links.append('*')
-        printouts = [self._print(i) for i in expr.args]
-        result = [printouts[0]]
-        for link, printout in zip(links, printouts[1:]):
-            result.extend([link, printout])
-        return '({0})'.format(''.join(result))
+    from sympy import MatrixExpr
+    links = []
+    for i, j in zip(self.args[1:], self.args[:-1]):
+        if isinstance(i, MatrixExpr) and isinstance(j, MatrixExpr):
+            links.append(').dot(')
+        else:
+            links.append('*')
+    printouts = [printer.doprint(i) for i in self.args]
+    result = [printouts[0]]
+    for link, printout in zip(links, printouts[1:]):
+        result.extend([link, printout])
+    return '({0})'.format(''.join(result))
+MatMul._numpycode = _print_MatMul
 
-    def _print_Inverse(self, expr):
-        return "%s(%s)" % (self._module_format('numpy.linalg.inv'),
-                           self._print(expr.args[0]))
+def _print_Inverse(self, printer):
+    return "%s(%s)" % (printer._module_format('numpy.linalg.inv'),
+                       printer.doprint(self.args[0]))
+Inverse._numpycode = _print_Inverse
 
-    def _print_HadamardProduct(self, expr):
-        return "%s*%s" % (self._print(expr.args[0]),
-                          self._print(expr.args[1]))
+def _print_HadamardProduct(self, printer):
+    return "%s*%s" % (printer.doprint(self.args[0]),
+                      printer.doprint(self.args[1]))
+HadamardProduct._numpycode = _print_HadamardProduct
 
-    def _print_DiracDelta(self, expr):
-        """
-        Replaces a DiracDelta(x) by np.inf if x == 0, and 0 otherwise. This is
-        wrong, but the only thing we can do by the time we are printing. To
-        prevent mistakes, integrate before printing.
-        """
-        return "{0}({1}, [{1} == 0 , {1} != 0], [{2}, 0])".format(
-                                        self._module_format('numpy.piecewise'),
-                                        self._print(expr.args[0]),
-                                        self._module_format('numpy.inf'))
+def _print_Idx(self, printer):
+    """
+    Print ``Idx`` objects.
+    """
+    return "{0}".format(printer.doprint(self.args[0]))
+Idx._numpycode = _print_Idx
 
-    def _print_Idx(self, expr):
-        """
-        Print ``Idx`` objects.
-        """
-        return "{0}".format(self._print(expr.args[0]))
+def _print_MatPow(self, printer):
+    if self.shape == (1, 1):
+        # Scalar, so we can take a normal power.
+        return printer._print_Pow(self)
+    else:
+        return printer._print_MatPow(self)
+MatPow._numpycode = _print_MatPow
 
-    def _print_MatPow(self, expr):
-        if expr.shape == (1, 1):
-            # Scalar, so we can take a normal power.
-            return self._print_Pow(expr)
-        else:
-            return super(SymfitNumPyPrinter, self)._print_MatPow(expr)
+#########################################################
+# End of monkeypatch                                    #
+#########################################################
\ No newline at end of file
diff --git a/symfit/core/support.py b/symfit/core/support.py
index 5659e87..f4d1d04 100644
--- a/symfit/core/support.py
+++ b/symfit/core/support.py
@@ -19,7 +19,7 @@ from sympy import symbols, MatrixExpr
 from sympy.core.expr import Expr
 
 from symfit.core.argument import Parameter, Variable
-from symfit.core.printing import SymfitNumPyPrinter
+import symfit.core.printing  # Overwrites some numpy printing
 
 if sys.version_info >= (3,0):
     import inspect as inspect_sig
@@ -120,8 +120,7 @@ def sympy_to_py(func, args):
     func = func.xreplace(derivatives)
     args = [derivatives[var] if isinstance(var, sympy.Derivative) else var
             for var in args]
-    lambdafunc = lambdify(args, func, printer=SymfitNumPyPrinter,
-                          dummify=False)
+    lambdafunc = lambdify(args, func, dummify=False)
     # Check if the names of the lambda function are what we expect
     signature = inspect_sig.signature(lambdafunc)
     sig_parameters = OrderedDict(signature.parameters)
diff --git a/symfit/distributions.py b/symfit/distributions.py
index 2336405..4adb8c2 100644
--- a/symfit/distributions.py
+++ b/symfit/distributions.py
@@ -10,7 +10,12 @@ import sympy
 
 def Gaussian(x, mu, sig):
     """
+    .. math::
+
+        f(x) = \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} e^{- \\frac{(x - \\mu)^2}{2 \\sigma^2}}
+
     Gaussian pdf.
+
     :param x: free variable.
     :param mu: mean of the distribution.
     :param sig: standard deviation of the distribution.
@@ -20,7 +25,8 @@ def Gaussian(x, mu, sig):
 
 def BivariateGaussian(x, y, mu_x, mu_y, sig_x, sig_y, rho):
     """
-    Bivariate Gaussian pdf.
+    `Bivariate Gaussian pdf
+    <http://mathworld.wolfram.com/BivariateNormalDistribution.html>`_.
 
     :param x: :class:`symfit.core.argument.Variable`
     :param y: :class:`symfit.core.argument.Variable`
@@ -41,7 +47,12 @@ def BivariateGaussian(x, y, mu_x, mu_y, sig_x, sig_y, rho):
 
 def Exp(x, l):
     """
+    .. math::
+
+        f(x) = l e^{- l x}
+
     Exponential Distribution pdf.
+
     :param x: free variable.
     :param l: rate parameter.
     :return: sympy.Expr for an Exponential Distribution pdf.


Test patch:
diff --git a/tests/test_auto_fit.py b/tests/test_auto_fit.py
index 53a2056..4f88207 100644
--- a/tests/test_auto_fit.py
+++ b/tests/test_auto_fit.py
@@ -5,7 +5,7 @@ import numpy as np
 
 from symfit import (
     variables, parameters, Fit, Parameter, Variable,
-    Equality, Model
+    Equality, Model, GradientModel
 )
 from symfit.core.minimizers import BFGS, MINPACK, SLSQP, LBFGSB
 from symfit.distributions import Gaussian
@@ -250,7 +250,7 @@ class TestAutoFit(unittest.TestCase):
         y = Variable('y')
         g = Variable('g')
 
-        model = Model({g: A * Gaussian(x, x0, sig_x) * Gaussian(y, y0, sig_y)})
+        model = GradientModel({g: A * Gaussian(x, x0, sig_x) * Gaussian(y, y0, sig_y)})
         fit = Fit(model, x=xx, y=yy, g=ydata)
         fit_result = fit.execute()
 
@@ -291,7 +291,7 @@ class TestAutoFit(unittest.TestCase):
         y = Variable('y')
         g = Variable('g')
 
-        model = Model({g: A * Gaussian(x, x0, sig_x) * Gaussian(y, y0, sig_y) + b})
+        model = GradientModel({g: A * Gaussian(x, x0, sig_x) * Gaussian(y, y0, sig_y) + b})
 
         # ydata, = model(x=xx, y=yy, x0=mean[0], y0=mean[1], sig_x=np.sqrt(cov[0][0]), sig_y=np.sqrt(cov[1][1]), A=1, b=3.0)
         fit = Fit(model, x=xx, y=yy, g=ydata)
diff --git a/tests/test_constrained.py b/tests/test_constrained.py
index 8174a8d..8429343 100644
--- a/tests/test_constrained.py
+++ b/tests/test_constrained.py
@@ -10,7 +10,7 @@ from scipy.optimize import NonlinearConstraint, minimize
 from symfit import (
     variables, Variable, parameters, Parameter, ODEModel,
     Fit, Equality, D, Model, log, FitResults, GreaterThan, Eq, Ge, Le,
-    CallableNumericalModel, HadamardProduct
+    CallableNumericalModel, HadamardProduct, GradientModel
 )
 from symfit.distributions import Gaussian
 from symfit.core.minimizers import (
@@ -524,7 +524,7 @@ class TestConstrained(unittest.TestCase):
         x = Variable('x')
         y = Variable('y')
         g = Variable('g')
-        model = Model({g: A * Gaussian(x, x0, sig_x) * Gaussian(y, y0, sig_y)})
+        model = GradientModel({g: A * Gaussian(x, x0, sig_x) * Gaussian(y, y0, sig_y)})
         fit = Fit(model, x=xx, y=yy, g=ydata)
         fit_result = fit.execute()
 
@@ -732,7 +732,7 @@ class TestConstrained(unittest.TestCase):
         i = Idx('i', (0, 1000))
         sig.min = 0.0
 
-        model = Model({y: A * Gaussian(x, mu=mu, sig=sig)})
+        model = GradientModel({y: A * Gaussian(x, mu=mu, sig=sig)})
 
         # Generate data, 100 samples from a N(1.2, 2) distribution
         np.random.seed(2)
diff --git a/tests/test_general.py b/tests/test_general.py
index 91d75cc..620157b 100644
--- a/tests/test_general.py
+++ b/tests/test_general.py
@@ -428,7 +428,7 @@ class Tests(unittest.TestCase):
         sig = Parameter('sig')
         A = Parameter('A')
         x = Variable('x')
-        g = A * Gaussian(x, x0, sig)
+        g = GradientModel(A * Gaussian(x, x0, sig))
 
         fit = Fit(g, xdata, ydata)
         fit_result = fit.execute()
@@ -486,13 +486,13 @@ class Tests(unittest.TestCase):
         A_2 = Parameter()
         g_2 = A_2 * Gaussian(x, x0_2, sig_x_2) * Gaussian(y, y0_2, sig_y_2)
 
-        model = g_1 + g_2
+        model = GradientModel(g_1 + g_2)
         fit = Fit(model, xx, yy, ydata)
         fit_result = fit.execute()
 
         self.assertIsInstance(fit.minimizer, LBFGSB)
 
-        img = model(x=xx, y=yy, **fit_result.params)
+        img = model(x=xx, y=yy, **fit_result.params)[0]
         img_g_1 = g_1(x=xx, y=yy, **fit_result.params)
         img_g_2 = g_2(x=xx, y=yy, **fit_result.params)
         np.testing.assert_array_equal(img, img_g_1 + img_g_2)
@@ -600,7 +600,7 @@ class Tests(unittest.TestCase):
         sig.value = 3.0
         mu.value = 50.
         x = Variable()
-        pdf = Gaussian(x, mu, sig)
+        pdf = GradientModel(Gaussian(x, mu, sig))
 
         np.random.seed(10)
         xdata = np.random.normal(51., 3.5, 10000)
diff --git a/tests/test_global_opt.py b/tests/test_global_opt.py
index e8ffc18..d6ce039 100644
--- a/tests/test_global_opt.py
+++ b/tests/test_global_opt.py
@@ -7,8 +7,7 @@ import sys
 import numpy as np
 
 from symfit import (
-    variables, parameters, Fit, Parameter, Variable,
-    Equality, Model
+    Fit, Parameter, Variable, Model, GradientModel
 )
 from symfit.core.minimizers import BFGS, DifferentialEvolution
 from symfit.distributions import Gaussian
@@ -48,7 +47,7 @@ class TestGlobalOptGaussian(unittest.TestCase):
         g_1 = self.A_1 * Gaussian(x, self.x0_1, self.sig_x_1) *\
                 Gaussian(y, self.y0_1, self.sig_y_1)
 
-        self.model = Model(g_1)
+        self.model = GradientModel(g_1)
     
     def test_diff_evo(self):
         """
diff --git a/tests/test_objectives.py b/tests/test_objectives.py
index e7635ed..6fcc1fa 100644
--- a/tests/test_objectives.py
+++ b/tests/test_objectives.py
@@ -8,17 +8,16 @@ import numpy as np
 from symfit import (
     Variable, Parameter, Eq, Ge, Le, Lt, Gt, Ne, parameters, ModelError, Fit,
     Model, FitResults, variables, CallableNumericalModel, Idx,
-    IndexedBase, symbols, Sum, log
+    IndexedBase, symbols, Sum, log, exp, cos, pi, besseli
 )
 from symfit.core.objectives import (
     VectorLeastSquares, LeastSquares, LogLikelihood, MinimizeModel
 )
 from symfit.core.fit_results import FitResults
-from symfit.core.printing import SymfitNumPyPrinter
 from symfit.distributions import Exp
 
 # Overwrite the way Sum is printed by numpy just while testing. Is not
-# general enough to be moved to SymfitNumPyPrinter, but has to be used
+# general enough to be moved to symfit.core.printing, but has to be used
 # in this test. This way of summing completely ignores the summation indices and
 # the dimensions, and instead just flattens everything to a scalar. Only used
 # in this test to build the analytical equivalents of our LeastSquares
@@ -28,11 +27,9 @@ class FlattenSum(Sum):
     Just a sum which is printed differently: by flattening the whole array and
     summing it. Used in tests only.
     """
-
-def _print_FlattenSum(self, expr):
-    return "%s(%s)" % (self._module_format('numpy.sum'),
-                       self._print(expr.function))
-SymfitNumPyPrinter._print_FlattenSum = _print_FlattenSum
+    def _numpycode(self, printer):
+        return "%s(%s)" % (printer._module_format('numpy.sum'),
+                           printer.doprint(self.function))
 
 
 class TestObjectives(unittest.TestCase):
@@ -189,7 +186,54 @@ class TestObjectives(unittest.TestCase):
         self.assertAlmostEqual(fit_exact_result.stdev(a), fit_num_result.stdev(a))
         self.assertAlmostEqual(fit_exact_result.stdev(b), fit_num_result.stdev(b))
 
+    def test_LogLikelihood_global(self):
+        """
+        This is a test for global likelihood fitting to multiple data sets.
+        Based on SO question 56006357.
+        """
+        # creating the data
+        mu1, mu2 = .05, -.05
+        sigma1, sigma2 = 3.5, 2.5
+        n1, n2 = 80, 90
+        np.random.seed(42)
+        x1 = np.random.vonmises(mu1, sigma1, n1)
+        x2 = np.random.vonmises(mu2, sigma2, n2)
+
+        n = 2  # number of components
+        xs = variables(
+            'x,' + ','.join('x_{}'.format(i) for i in range(1, n + 1)))
+        x, xs = xs[0], xs[1:]
+        ys = variables(','.join('y_{}'.format(i) for i in range(1, n + 1)))
+        mu, kappa = parameters('mu, kappa')
+        kappas = parameters(','.join('k_{}'.format(i) for i in range(1, n + 1)),
+                            min=0, max=10)
+        mu.min, mu.max = - np.pi, np.pi
+
+        template = exp(kappa * cos(x - mu)) / (2 * pi * besseli(0, kappa))
+
+        model = Model(
+            {y_i: template.subs({kappa: k_i, x: x_i}) for y_i, x_i, k_i in
+             zip(ys, xs, kappas)}
+        )
+
+        all_data = {xs[0]: x1, xs[1]: x2, ys[0]: None, ys[1]: None}
+        all_params = {'mu': 1}
+        all_params.update({k_i.name: 1 for k_i in kappas})
+
+        # Evaluate the loglikelihood and its jacobian and hessian
+        logL = LogLikelihood(model, data=all_data)
+        eval_numerical = logL(**all_params)
+        jac_numerical = logL.eval_jacobian(**all_params)
+        hess_numerical = logL.eval_hessian(**all_params)
+
+        # Test the types and shapes of the components.
+        self.assertIsInstance(eval_numerical, float)
+        self.assertIsInstance(jac_numerical, np.ndarray)
+        self.assertIsInstance(hess_numerical, np.ndarray)
 
+        self.assertEqual(eval_numerical.shape, tuple())  # Empty tuple -> scalar
+        self.assertEqual(jac_numerical.shape, (3,))
+        self.assertEqual(hess_numerical.shape, (3, 3,))
 
 if __name__ == '__main__':
     try: