Repository: fitbenchmarking/fitbenchmarking

Issue:
Create a Jacobian and residual evaluation function in the fitting problem
When adding RALFit (#215) you have to supply Jacobians and residuals to the software. Currently the Jacobians is done by using `SciPy` in the RALFit controller. We should create two new functions `eval_j` and `eval_r` which evaluate the Jacobian and residual (`y - eval_f(x)`) in the fitting problem, similar to `eval_f`. 

Blocked by #283  

Discussion:

Patch:
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 96dba92b..c39956f4 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -16,10 +16,9 @@
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
 #
-# import os
-# import sys
-# sys.path.insert(0, os.path.abspath("."))
-
+import os
+import sys
+sys.path.insert(0, os.path.abspath('../../'))
 
 # -- General configuration ------------------------------------------------
 
@@ -31,7 +30,8 @@
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
 extensions = ['sphinx.ext.autodoc',
-    'sphinx.ext.todo']
+              'sphinx.ext.todo',
+              'sphinx.ext.imgmath']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
diff --git a/docs/source/contributors/extending_fitbenchmarking.rst b/docs/source/contributors/extending_fitbenchmarking.rst
index dd964efa..74dc2c31 100644
--- a/docs/source/contributors/extending_fitbenchmarking.rst
+++ b/docs/source/contributors/extending_fitbenchmarking.rst
@@ -112,6 +112,9 @@ written so that the fitting is separated from the preparation wherever possible
 in order to give accurate timings for the fitting. Examples of these
 controllers can be found in ``fitbenchmarking/fitting/controllers``.
 
+Fitting software requirements
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
 In order to add a new controller, you will need to:
 
 1. Give the software a name `<software_name>` this will be used by users when
@@ -142,5 +145,18 @@ In order to add a new controller, you will need to:
    Unless the new controller is more complicated than the currently available
    controllers, this can be done by following the example of the others.
 
-Note: For ease of maintainance, please add new controllers in alphabetical
-      order, whenever it is added to a list of softwares.
\ No newline at end of file
+.. note::
+   For ease of maintainance, please add new controllers to a list of software in alphabetical order.
+   
+   
+The :class:`~fitbenchmarking.parsing.fitting_problem.FittingProblem` class
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+When adding new minimizers, you will find it helpful to make use of the
+following members of the :class:`~fitbenchmarking.parsing.fitting_problem.FittingProblem` class:
+
+.. currentmodule:: fitbenchmarking.parsing.fitting_problem
+.. autoclass:: fitbenchmarking.parsing.fitting_problem.FittingProblem
+   :members: eval_f, eval_r, eval_r_norm, eval_j,
+	     data_x, data_y, data_e, starting_values
+	  
diff --git a/fitbenchmarking/fitting/controllers/dfogn_controller.py b/fitbenchmarking/fitting/controllers/dfogn_controller.py
index 1ae8f73f..e56c6c1a 100644
--- a/fitbenchmarking/fitting/controllers/dfogn_controller.py
+++ b/fitbenchmarking/fitting/controllers/dfogn_controller.py
@@ -30,20 +30,13 @@ class DFOGNController(Controller):
         """
         self._pinit = np.asarray(self.initial_params)
 
-    def _prediction_error(self, p):
-        f = self.data_y - self.problem.eval_f(params=p)
-        if self.use_errors:
-            f = f/self.data_e
-
-        return f
-
     def fit(self):
         """
         Run problem with DFO-GN.
         """
         self.success = False
 
-        self._soln = dfogn.solve(self._prediction_error,
+        self._soln = dfogn.solve(self.problem.eval_r,
                                  self._pinit)
 
         if (self._soln.flag == 0):
diff --git a/fitbenchmarking/fitting/controllers/minuit_controller.py b/fitbenchmarking/fitting/controllers/minuit_controller.py
index 50a82c59..94772446 100644
--- a/fitbenchmarking/fitting/controllers/minuit_controller.py
+++ b/fitbenchmarking/fitting/controllers/minuit_controller.py
@@ -37,19 +37,11 @@ class MinuitController(Controller):
         self._initial_step = 0.1 * np.array(self.initial_params)
         # set small steps to something sensible(?)
         self._initial_step[self._initial_step < 1e-12] = 1e-12
-        self._minuit_problem = Minuit.from_array_func(self._prediction_error,
+        self._minuit_problem = Minuit.from_array_func(self.problem.eval_r_norm,
                                                       self.initial_params,
                                                       error=self._initial_step,
                                                       errordef=1)
 
-    def _prediction_error(self, p):
-        f = self.problem.eval_f(params=p)
-        f = f - self.data_y
-        if self.use_errors:
-            f = f / self.data_e
-
-        return np.dot(f, f)
-
     def fit(self):
         """
         Run problem with Minuit
diff --git a/fitbenchmarking/fitting/controllers/ralfit_controller.py b/fitbenchmarking/fitting/controllers/ralfit_controller.py
index b76971b3..abd8fa4f 100644
--- a/fitbenchmarking/fitting/controllers/ralfit_controller.py
+++ b/fitbenchmarking/fitting/controllers/ralfit_controller.py
@@ -4,7 +4,6 @@ https://github.com/ralna/RALFit
 """
 
 import ral_nlls
-from scipy.optimize._numdiff import approx_derivative
 
 from fitbenchmarking.fitting.controllers.base_controller import Controller
 
@@ -43,27 +42,14 @@ class RALFitController(Controller):
         else:
             raise RuntimeError("An undefined RALFit minmizer was selected")
 
-    def _prediction_error(self, p):
-        f = self.problem.eval_f(params=p)
-        f = f - self.data_y
-        if self.use_errors:
-            f = f / self.data_e
-
-        return f
-
-    def _jac(self, p):
-        j = approx_derivative(self._prediction_error,
-                              p)
-        return j
-
     def fit(self):
         """
         Run problem with RALFit.
         """
         self.success = False
         self._popt = ral_nlls.solve(self.initial_params,
-                                    self._prediction_error,
-                                    self._jac,
+                                    self.problem.eval_r,
+                                    self.problem.eval_j,
                                     options=self._options)[0]
 
         self.success = (self._popt is not None)
diff --git a/fitbenchmarking/parsing/fitting_problem.py b/fitbenchmarking/parsing/fitting_problem.py
index 06c2dfeb..2cf4b83f 100644
--- a/fitbenchmarking/parsing/fitting_problem.py
+++ b/fitbenchmarking/parsing/fitting_problem.py
@@ -11,51 +11,66 @@ except ImportError:
     # python3
     from itertools import zip_longest as izip_longest
 import numpy as np
+from scipy.optimize._numdiff import approx_derivative
 
 
 class FittingProblem:
-    """
-    Definition of a fitting problem, normally populated by a parser from a
+    r"""
+    Definition of a fitting problem, which will be populated by a parser from a
     problem definition file.
 
-    Types of data:
-        - strings: name, equation
-        - floats: start_x, end_x
-        - numpy arrays: data_x, data_y, data_e
-        - arrays: starting_values, value_ranges, functions
+    This defines a fitting problem where, given a set of :math:`n` data points
+    :math:`(x_i,y_i)`, associated errors :math:`e_i`, and a model
+    function :math:`f(x,p)`, we find the optimal parameters in the
+    least-squares sense by solving:
+
+    .. math:: \min_p \sum_{i=1}^n \left( \frac{y_i - f(x_i, p)}{e_i} \right)^2
+
+    where :math:`p` is a vector of length :math:`m`, and we start from a given
+    intial guess for the optimal parameters.
     """
 
     def __init__(self):
 
-        # Name (title) of the fitting problem
-        # str
+        #: *string* Name (title) of the fitting problem
         self.name = None
 
-        # Equation (function or model) to fit against data
-        # string
+        #: *string* Equation (function or model) to fit against data
         self.equation = None
 
-        # Define range to fit model data over if different from entire range
-        # of data
-        # floats
+        #: *float* The start of the range to fit model data over
+        #: (if different from entire range)
         self.start_x = None
+
+        #: *float* The end of the range to fit model data over
+        #: (if different from entire range) (/float/)
         self.end_x = None
 
-        # The data
-        # numpy array of floats
+        #: *numpy array* The x-data
         self.data_x = None
+
+        #: *numpy array* The y-data
         self.data_y = None
+
+        #: *numpy array* The errors
         self.data_e = None
 
-        # Starting values of the fitting parameters
-        # list of dict -> [{p1_name: p1_val1, p2_name: p2_val1, ...},
-        #                  {p1_name: p1_val2, ...},
-        #                  ...]
+        #: *list of dict*
+        #: Starting values of the fitting parameters
+        #:
+        #: e.g.
+        #: :code:`[{p1_name: p1_val1, p2_name: p2_val1, ...},
+        #: {p1_name: p1_val2, ...}, ...]`
         self.starting_values = None
-        # dict -> {p1_name: [p1_min, p1_max], ...}
+
+        #: *dict*
+        #: Smallest and largest values of interest in the data
+        #:
+        #: e.g.
+        #: :code:`{p1_name: [p1_min, p1_max], ...}`
         self.value_ranges = None
 
-        # Callable function
+        #: Callable function
         self.function = None
 
         self._param_names = None
@@ -95,6 +110,73 @@ class FittingProblem:
             x = self.data_x
         return self.function(x, *params)
 
+    def eval_r(self, params, x=None, y=None, e=None):
+        """
+        Calculate residuals and weight them if using errors
+
+        :param params: The parameters to calculate residuals for
+        :type params: list
+        :param x: x data points, defaults to self.data_x
+        :type x: numpy array, optional
+        :param y: y data points, defaults to self.data_y
+        :type y: numpy array, optional
+        :param e: error at each data point, defaults to self.data_e
+        :type e: numpy array, optional
+
+        :return: The residuals for the datapoints at the given parameters
+        :rtype: numpy array
+        """
+
+        if x is None and y is None and e is None:
+            x = self.data_x
+            y = self.data_y
+            e = self.data_e
+        elif x is None or y is None:
+            raise ValueError('Residuals could not be computed with only one'
+                             'of x and y')
+
+        result = y - self.eval_f(params=params, x=x)
+        if e is not None:
+            result = result / e
+        return result
+
+    def eval_r_norm(self, params, x=None, y=None, e=None):
+        """
+        Evaluate the square of the L2 norm of the residuals
+
+        :param params: The parameters to calculate residuals for
+        :type params: list
+        :param x: x data points, defaults to self.data_x
+        :type x: numpy array, optional
+        :param y: y data points, defaults to self.data_y
+        :type y: numpy array, optional
+        :param e: error at each data point, defaults to self.data_e
+        :type e: numpy array, optional
+
+        :return: The sum of squares of residuals for the datapoints at the
+                 given parameters
+        :rtype: numpy array
+        """
+        r = self.eval_r(params=params, x=x, y=y, e=e)
+        return np.dot(r, r)
+
+    def eval_j(self, params, func=None, **kwargs):
+        """
+        Approximate the Jacobian using scipy for a given function at a given
+        point.
+
+        :param params: The parameter values to find the Jacobian at
+        :type params: list
+        :param func: Function to find the Jacobian for, defaults to self.eval_r
+        :type func: Callable, optional
+        :return: Approximation of the Jacobian
+        :rtype: numpy array
+        """
+        if func is None:
+            func = self.eval_r
+
+        return approx_derivative(func, params, kwargs=kwargs)
+
     def eval_starting_params(self, param_set):
         """
         Evaluate the function using the starting parameters.


Test patch:
diff --git a/fitbenchmarking/parsing/tests/test_fitting_problem.py b/fitbenchmarking/parsing/tests/test_fitting_problem.py
index ae1aeb8e..df6b5403 100644
--- a/fitbenchmarking/parsing/tests/test_fitting_problem.py
+++ b/fitbenchmarking/parsing/tests/test_fitting_problem.py
@@ -69,6 +69,83 @@ class TestFittingProblem(TestCase):
         eval_result = fitting_problem.eval_f(params=[5])
         self.assertTrue(all(eval_result == np.array([25, 26, 27])))
 
+    def test_eval_r(self):
+        """
+        Test that eval_r is correct
+        """
+        fitting_problem = FittingProblem()
+        self.assertRaises(ValueError,
+                          fitting_problem.eval_r,
+                          params=[1, 2, 3],
+                          x=2)
+
+        fitting_problem.function = lambda x, p1: x + p1
+        x_val = np.array([1, 8, 11])
+        y_val = np.array([6, 10, 20])
+
+        eval_result = fitting_problem.eval_r(x=x_val,
+                                             y=y_val,
+                                             params=[5])
+        self.assertTrue(all(eval_result == np.array([0, -3, 4])))
+
+        e_val = np.array([2, 4, 1])
+        eval_result = fitting_problem.eval_r(x=x_val,
+                                             y=y_val,
+                                             e=e_val,
+                                             params=[5])
+        self.assertTrue(all(eval_result == np.array([0, -0.75, 4])))
+
+        fitting_problem.data_x = np.array([20, 21, 22])
+        fitting_problem.data_y = np.array([20, 30, 35])
+        eval_result = fitting_problem.eval_r(params=[5])
+        self.assertTrue(all(eval_result == np.array([-5, 4, 8])))
+
+        fitting_problem.data_e = np.array([2, 5, 10])
+        eval_result = fitting_problem.eval_r(params=[5])
+        self.assertTrue(all(eval_result == np.array([-2.5, 0.8, 0.8])))
+
+    def test_eval_r_norm(self):
+        """
+        Test that eval_r_norm is correct
+        """
+        fitting_problem = FittingProblem()
+        fitting_problem.function = lambda x, p1: x + p1
+        x_val = np.array([1, 8, 11])
+        y_val = np.array([6, 10, 20])
+        e_val = np.array([0.5, 10, 0.1])
+
+        eval_result = fitting_problem.eval_r_norm(params=[5],
+                                                  x=x_val,
+                                                  y=y_val,
+                                                  e=e_val)
+        self.assertEqual(eval_result, 1600.09)
+
+        fitting_problem.data_x = x_val
+        fitting_problem.data_y = y_val
+        eval_result = fitting_problem.eval_r_norm(params=[5])
+        self.assertEqual(eval_result, 25)
+
+    def test_eval_j(self):
+        """
+        Test that eval_j is correct
+        """
+        def f(x, p1, p2):
+            return p1*np.exp(p2*x)
+
+        def J(x, p):
+            return np.column_stack((-np.exp(p[1]*x),
+                                    -x*p[0]*np.exp(p[1]*x)))
+
+        fitting_problem = FittingProblem()
+        fitting_problem.function = f
+        fitting_problem.data_x = np.array([1, 2, 3, 4, 5])
+        fitting_problem.data_y = np.array([1, 2, 4, 8, 16])
+        params = [6, 0.1]
+
+        eval_result = fitting_problem.eval_j(params=params)
+        actual = J(x=fitting_problem.data_x, p=params)
+        self.assertTrue(np.isclose(actual, eval_result).all())
+
     def test_eval_starting_params(self):
         """
         Test that eval_starting_params returns the correct result