Added parser for inputting sparse jacobian

docs/source/users/options/jacobian_option.rst

@@ -9,13 +9,16 @@ The Jacobian section allows you to control which methods for computing Jacobians
 Analytic (:code:`analytic`)
 ---------------------------
 
-Analytic Jacobians can only be used for specific :ref:`problem_def`.
-Currently the supported formats are CUTEst, Mantid, and NIST.
-The only option is:
+Functions for analytic Jacobians are available for specific :ref:`problem_def`. Specifically, dense jacobian functions are
+available for cutest, NIST and Mantid, while sparse jacobians are available for cutest only. However, the user can specify custom
+jacobian functions for any :ref:`problem_def` except from NIST and cutest, through a problem definition file. They can specify a dense
+and a sparse jacobian, and select which one to use by the following options:
 
-* ``default`` - use the analytic derivative provided by a supported format.
+* ``default`` - use the analytic (dense) derivative provided by a supported format or the dense jacobian function provided by the user.
+* ``sparse`` - use a sparse jacobian function, either available or provided by the user, to compute the derivative.
 
-Default is ``default``
+Default is ``default``. When dealing with problems of supported formats, if a jacobian function is specified, this will
+have priority over the analytic derivative provided by the format.
 
 .. code-block:: rst
 
@@ -38,6 +41,7 @@ options are:
 * ``2-point`` - use the first order accuracy forward or backward difference.
 * ``3-point`` - use central difference in interior points and the second order accuracy forward or backward difference near the boundary.
 * ``cs`` - use a complex-step finite difference scheme. This assumes that the user function is real-valued and can be analytically continued to the complex plane. Otherwise, produces bogus results.
+* ``2-point_sparse`` - use 2-point, with a sparsity pattern.
 
 Default is ``2-point``
 

docs/source/users/problem_definition_files/native.rst

@@ -25,7 +25,7 @@ keys are described below:
 
 software
   Either 'IVP', 'Mantid', 'SasView', or 'Horace' (case insensitive).
-  
+
   This defines whether to use an IVP format, Mantid, or SasView to generate the model.
   The 'Mantid' software also supports Mantid's MultiFit functionality, which
   requires the parameters listed here to be defined slightly differently.
@@ -72,7 +72,7 @@ input_file
   ``examples/benchmark_problems/Data_Assimilation/data_files/lorentz.txt``
 
 plot_scale
-  The scale of the x and y axis for the plots. The options are 'loglog', 'logy', 'logx' and 'linear'. If this 
+  The scale of the x and y axis for the plots. The options are 'loglog', 'logy', 'logx' and 'linear'. If this
   is not set it will default to 'linear'.
 
 function
@@ -121,10 +121,32 @@ function
 
   SASView functions can be any of
   `these <http://www.sasview.org/docs/user/qtgui/Perspectives/Fitting/models/index.html>`__.
-  
+
   **Horace**
 
-  The Horace functions are defined here :ref:`horace_format`
+  The Horace functions are defined here :ref:`horace_format` .
+
+jacobian
+  This is to define a dense jacobian function, or a sparse jacobian function, or both.
+  The parser for this function allows the user to define ``g`` in the following:
+
+  .. math:: \nabla_p f(x, *args) = g(x, *args)
+
+  To do this we use a python module, where the user can define a dense jacobian function
+  and a sparse jacobian function (or just one of the two). As in the above formula,
+  both functions can take the following arguments:
+
+  - *x* (np.array): A value for x to evaluate at
+  - *\*args* (floats): The parameters to fit
+
+  To link to this functions we use a string with the following parameters:
+
+  - *module*: The path to the module
+  - *dense_func*: The name of the dense jacobian function within the module
+  - *sparse_func*: The name of the sparse jacobian function within the module
+
+The sparse jacobian function provided must return a matrix in sparse format
+(e.g. coo, csr, crs), otherwise an error will be thrown.
 
 fit_ranges
   This specifies the region to be fit.
@@ -133,7 +155,7 @@ fit_ranges
   is the minimum in the range and the second is the maximum.
 
 parameter_ranges
-  An optional setting which specifies upper and lower bounds for 
+  An optional setting which specifies upper and lower bounds for
   parameters in the problem.
 
   Similarly to ``fit_ranges``, it takes the form where the first number

fitbenchmarking/cost_func/hellinger_nlls_cost_func.py

@@ -1,7 +1,7 @@
 """
 Implements the root non-linear least squares cost function
 """
-from numpy import array, matmul, ravel, sqrt
+from numpy import array, ravel, sqrt
 
 from fitbenchmarking.cost_func.nlls_base_cost_func import BaseNLLSCostFunc
 from fitbenchmarking.utils.exceptions import (CostFuncError,
@@ -85,7 +85,7 @@ def hes_res(self, params, **kwargs):
 
         for i in range(len(x)):
             jac_i = array([jac[i]])
-            hes[:, :, i] = matmul(jac_i.T, jac_i) / (4 * f[i] ** (3/2)) \
+            hes[:, :, i] = jac_i.T.dot(jac_i) / (4 * f[i] ** (3/2)) \
                 - hes[:, :, i] / (2 * f[i] ** (1/2))
         return hes, self.jac_res(params, **kwargs)
 

fitbenchmarking/cost_func/nlls_base_cost_func.py

@@ -83,7 +83,7 @@ def jac_cost(self, params, **kwargs):
         r = self.eval_r(params, **kwargs)
         J = self.jac_res(params, **kwargs)
 
-        return 2.0 * matmul(J.T, r)
+        return 2.0 * J.T.dot(r)
 
     def hes_cost(self, params, **kwargs):
         """

fitbenchmarking/cost_func/poisson_cost_func.py

@@ -118,7 +118,7 @@ def hes_res(self, params, **kwargs):
         for i in range(len(x)):
             jac_i = np.array([jac[i]])
             hes[:, :, i] = hes[:, :, i] - y[i] / f[i] * \
-                (hes[:, :, i] - np.matmul(jac_i.T, jac_i) / f[i])
+                (hes[:, :, i] - jac_i.T.dot(jac_i) / f[i])
         return hes, self.jac_res(params, **kwargs)
 
     def hes_cost(self, params, **kwargs):

fitbenchmarking/jacobian/analytic_jacobian.py

@@ -1,8 +1,12 @@
 """
 Module which acts as a analytic Jacobian calculator
 """
+from scipy.sparse import issparse
+
 from fitbenchmarking.jacobian.base_jacobian import Jacobian
-from fitbenchmarking.utils.exceptions import NoJacobianError
+from fitbenchmarking.utils.exceptions import (NoJacobianError,
+                                              NoSparseJacobianError,
+                                              SparseJacobianIsDenseError)
 
 
 class Analytic(Jacobian):
@@ -27,6 +31,22 @@ def eval(self, params, **kwargs):
         :rtype: numpy array
         """
         x = kwargs.get("x", self.problem.data_x)
+
+        if self.method == "sparse":
+
+            if self.problem.sparse_jacobian is None:
+                raise NoSparseJacobianError(
+                    f'The selected method is {self.method} but the '
+                    'sparse_jacobian function is None. Please provide a '
+                    'sparse jacobian function.')
+
+            sparse_jac = self.problem.sparse_jacobian(x, params)
+
+            if not issparse(sparse_jac):
+                raise SparseJacobianIsDenseError()
+
+            return sparse_jac
+
         return self.problem.jacobian(x, params)
 
     def name(self) -> str:

fitbenchmarking/jacobian/scipy_jacobian.py

@@ -3,8 +3,12 @@
 """
 import numpy as np
 from scipy.optimize._numdiff import approx_derivative
+from scipy.sparse import issparse
 
 from fitbenchmarking.jacobian.base_jacobian import Jacobian
+from fitbenchmarking.utils.exceptions import (NoSparseJacobianError,
+                                              SparseJacobianIsDenseError)
+from fitbenchmarking.utils.log import get_logger
 
 
 class Scipy(Jacobian):
@@ -15,6 +19,11 @@ class Scipy(Jacobian):
     # Problem formats that are incompatible with certain Scipy Jacobians
     INCOMPATIBLE_PROBLEMS = {"cs": ["mantid"]}
 
+    def __init__(self, problem):
+        super().__init__(problem)
+        self.jac_pattern = None
+        self.equiv_np_method = None
+
     def eval(self, params, **kwargs):
         """
         Evaluates Jacobian of problem.eval_model
@@ -25,9 +34,40 @@ def eval(self, params, **kwargs):
         :return: Approximation of the Jacobian
         :rtype: numpy array
         """
+
+        self.equiv_np_method = self.method
+
+        LOGGER = get_logger()
+
+        if self.method.endswith("_sparse"):
+
+            if self.problem.sparse_jacobian is None:
+                raise NoSparseJacobianError(
+                    f'The selected method is {self.method} but the '
+                    'sparse_jacobian function is None. Please provide a '
+                    'sparse jacobian function.')
+
+            self.jac_pattern = self.problem.\
+                sparse_jacobian(self.problem.data_x, params)
+
+            if not issparse(self.jac_pattern):
+                raise SparseJacobianIsDenseError()
+
+            # Remove the "_sparse" at the end
+            self.equiv_np_method = self.method[:-7]
+
+        else:
+            if self.problem.sparse_jacobian is not None:
+                LOGGER.info('Sparse_jacobian function found, but it '
+                            'will not be used as the selected method is %s.',
+                            self.method)
+
         func = self.problem.eval_model
-        jac = approx_derivative(func, params, method=self.method,
+        jac = approx_derivative(func, params,
+                                method=self.equiv_np_method,
                                 rel_step=None,
                                 bounds=(-np.inf, np.inf),
-                                kwargs=kwargs)
+                                kwargs=kwargs,
+                                sparsity=self.jac_pattern)
+
         return jac

fitbenchmarking/jacobian/tests/test_jacobian.py

@@ -4,6 +4,8 @@
 from unittest import TestCase
 
 import numpy as np
+from scipy import sparse
+from scipy.sparse import issparse
 
 from fitbenchmarking.cost_func.hellinger_nlls_cost_func import \
     HellingerNLLSCostFunc
@@ -18,6 +20,8 @@
 from fitbenchmarking.jacobian.scipy_jacobian import Scipy
 from fitbenchmarking.parsing.fitting_problem import FittingProblem
 from fitbenchmarking.utils import exceptions
+from fitbenchmarking.utils.exceptions import (NoSparseJacobianError,
+                                              SparseJacobianIsDenseError)
 from fitbenchmarking.utils.options import Options
 
 
@@ -54,6 +58,21 @@ def j(x, p):
                             x * p[0] * np.exp(p[1] * x)))
 
 
+def j_sparse(x, p):
+    """
+    Sparse Jacobian evaluation
+
+    :param x: x data points, defaults to self.data_x
+    :type x: numpy array, optional
+    :param p: list of parameters to fit
+    :type p: list
+
+    :return: Sparse Jacobian evaluation
+    :rtype: scipy.sparse.csr_matrix
+    """
+    return sparse.csr_matrix(j(x, p))
+
+
 def J(x, p):
     """
     Analytic Jacobian residuals evaluation
@@ -145,6 +164,7 @@ def setUp(self):
         self.fitting_problem = FittingProblem(options)
         self.fitting_problem.function = f
         self.fitting_problem.jacobian = j
+        self.fitting_problem.sparse_jacobian = j_sparse
         self.fitting_problem.data_x = np.array([1, 2, 3, 4, 5])
         self.fitting_problem.data_y = np.array([1, 2, 4, 8, 16])
         self.cost_func = NLLSCostFunc(self.fitting_problem)
@@ -173,6 +193,40 @@ def test_scipy_eval(self):
             eval_result = jac.eval(params=self.params)
             self.assertTrue(np.isclose(self.actual, eval_result).all())
 
+    def test_scipy_eval_returns_correct_sparse_jac(self):
+        """
+        Test whether Scipy evaluation is correct when
+        using sparsity
+        """
+        jac = Scipy(self.cost_func.problem)
+        jac.method = '2-point_sparse'
+        eval_result = jac.eval(params=self.params)
+
+        self.assertTrue(np.isclose(self.actual, eval_result.todense()).all())
+        self.assertTrue(issparse(eval_result))
+
+    def test_scipy_eval_raises_error_sparsej_dense(self):
+        """
+        Test that Scipy evaluation raises error when
+        result of sparse jacobian is not in sparse format
+        """
+        self.fitting_problem.sparse_jacobian = j
+        jac = Scipy(self.cost_func.problem)
+        jac.method = '2-point_sparse'
+        with self.assertRaises(SparseJacobianIsDenseError):
+            jac.eval(params=self.params)
+
+    def test_scipy_eval_raises_error_no_sparsej(self):
+        """
+        Test that Scipy evaluation raises error when
+        result of sparse jacobian is None
+        """
+        self.fitting_problem.sparse_jacobian = None
+        jac = Scipy(self.cost_func.problem)
+        jac.method = '2-point_sparse'
+        with self.assertRaises(NoSparseJacobianError):
+            jac.eval(params=self.params)
+
     def test_numdifftools_eval(self):
         """
         Test whether numdifftools evaluation is correct
@@ -198,6 +252,43 @@ def test_analytic_cutest_no_errors(self):
         actual = J(self.fitting_problem.data_x, self.params)
         self.assertTrue(np.isclose(actual, eval_result).all())
 
+    def test_analytic_eval_returns_correct_sparse_jac(self):
+        """
+        Test whether analytic jac evaluation is correct
+        when using sparsity
+        """
+        jac = Analytic(self.cost_func.problem)
+        jac.method = 'sparse'
+        self.cost_func.jacobian = jac
+        eval_result = self.cost_func.jac_res(params=self.params)
+        actual = J(self.fitting_problem.data_x, self.params)
+        self.assertTrue(np.isclose(actual, eval_result.todense()).all())
+        self.assertTrue(issparse(eval_result))
+
+    def test_analytic_eval_raises_error_no_sparsej(self):
+        """
+        Test that analytic jac evaluation raises error when
+        result of sparse jacobian is None
+        """
+        self.fitting_problem.sparse_jacobian = None
+        jac = Analytic(self.cost_func.problem)
+        jac.method = 'sparse'
+        self.cost_func.jacobian = jac
+        with self.assertRaises(NoSparseJacobianError):
+            self.cost_func.jac_res(params=self.params)
+
+    def test_analytic_eval_raises_error_sparsej_dense(self):
+        """
+        Test that analytic jac evaluation raises error when
+        result of sparse jacobian is not in sparse format
+        """
+        self.fitting_problem.sparse_jacobian = j
+        jac = Analytic(self.cost_func.problem)
+        jac.method = 'sparse'
+        self.cost_func.jacobian = jac
+        with self.assertRaises(SparseJacobianIsDenseError):
+            self.cost_func.jac_res(params=self.params)
+
     def test_analytic_cutest_weighted(self):
         """
         Test analytic Jacobian