Create plot_shi_xie_xuan_nocedal_general_benchmark

doc/examples/examples.rst

@@ -45,3 +45,4 @@ Shi, Xie, Xuan & Nocedal
    ../auto_example/plot_shi_xie_xuan_nocedal_forward
    ../auto_example/plot_shi_xie_xuan_nocedal_forward_benchmark
    ../auto_example/plot_shi_xie_xuan_nocedal_general
+   ../auto_example/plot_shi_xie_xuan_nocedal_general_benchmark

doc/examples/plot_shi_xie_xuan_nocedal_forward_benchmark.py

@@ -26,7 +26,7 @@ class on a collection of test problems.
 class ShiXieXuanNocedalForwardMethod:
     def __init__(self, relative_precision, initial_step):
         """
-        Create a ShiXieXuanNocedal method to compute the approximate first derivative
+        Create a ShiXieXuanNocedalForward method to compute the approximate first derivative
 
         Parameters
         ----------
@@ -229,9 +229,9 @@ def compute_first_derivative(self, function, x):
 
 # %%
 # Notice that the method cannot perform correctly for the sin function
-# at the point
+# at the point :math:`x = \pm \pi`.
 # Indeed, this function is such that :math:`f''(x) = 0` if :math:`x = \pm \pi`.
-# In this case, the test ratio and the method cannot work.
+# In this case, the test ratio is zero and the method cannot work.
 # Therefore, we make so that the points :math:`\pm \pi` are excluded from the benchmark.
 # The same problem appears at the point :math:`x = 0`.
 # This point is not included in the test set if the number of points is even

doc/examples/plot_shi_xie_xuan_nocedal_general_benchmark.py

@@ -0,0 +1,265 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# Copyright 2024 - Michaël Baudin.
+"""
+Benchmark Shi, Xie, Xuan & Nocedal's general method
+===================================================
+
+The goal of this example is to benchmark the :class:`~numericalderivative.ShiXieXuanNocedalGeneral`
+class on a collection of test problems.
+These problems are created by the :meth:`~numericalderivative.build_benchmark()` 
+static method, which returns a list of problems.
+
+"""
+
+# %%
+import numpy as np
+import pylab as pl
+import tabulate
+import numericalderivative as nd
+
+# %%
+# Compute the first derivative
+# ----------------------------
+
+
+class ShiXieXuanNocedalGeneralMethod:
+    def __init__(
+        self,
+        relative_precision,
+        initial_step,
+        differentiation_order=1,
+        formula_accuracy=2,
+        direction="central",
+    ):
+        """
+        Create a ShiXieXuanNocedalGeneral method to compute the approximate first derivative
+
+        Parameters
+        ----------
+        relative_precision : float, > 0, optional
+            The relative precision of evaluation of f.
+        initial_step : float, > 0
+            The initial step in the algorithm.
+        """
+        self.relative_precision = relative_precision
+        self.initial_step = initial_step
+        self.differentiation_order = differentiation_order
+        self.formula_accuracy = formula_accuracy
+        self.direction = direction
+
+    def compute_derivative(self, function, x):
+        """
+        Compute the first derivative using ShiXieXuanNocedal
+
+        Parameters
+        ----------
+        function : function
+            The function
+        x : float
+            The test point
+
+        Returns
+        -------
+        f_derivative_approx : float
+            The approximate value of the first derivative of the function at point x
+        number_of_function_evaluations : int
+            The number of function evaluations.
+        """
+        formula = nd.GeneralFiniteDifference(
+            function,
+            x,
+            self.differentiation_order,
+            self.formula_accuracy,
+            self.direction,
+        )
+        absolute_precision = abs(function(x)) * self.relative_precision
+        algorithm = nd.ShiXieXuanNocedalGeneral(
+            formula,
+            absolute_precision,
+        )
+        step, _ = algorithm.find_step(self.initial_step)
+        f_derivative_approx = algorithm.compute_derivative(step)
+        number_of_function_evaluations = algorithm.get_number_of_function_evaluations()
+        return f_derivative_approx, number_of_function_evaluations
+
+
+# %%
+# The next example computes the approximate derivative on the
+# :class:`~numericalderivative.ExponentialProblem`.
+
+# %%
+initial_step = 1.0e0
+problem = nd.ExponentialProblem()
+print(problem)
+function = problem.get_function()
+x = problem.get_x()
+differentiation_order = 1  # First derivative
+formula_accuracy = 2  # Order 2
+formula = nd.GeneralFiniteDifference(
+    function,
+    x,
+    differentiation_order,
+    formula_accuracy,
+    direction="central",  # Central formula
+)
+algorithm = nd.ShiXieXuanNocedalGeneral(
+    formula,
+    verbose=True,
+)
+third_derivative = problem.get_third_derivative()
+third_derivative_value = third_derivative(x)
+optimal_step, absolute_error = nd.FirstDerivativeCentral.compute_step(
+    third_derivative_value
+)
+print("Exact h* = %.3e" % (optimal_step))
+function = problem.get_function()
+first_derivative = problem.get_first_derivative()
+x = 1.0
+relative_precision = 1.0e-15
+absolute_precision = abs(function(x)) * relative_precision
+method = ShiXieXuanNocedalGeneralMethod(absolute_precision, initial_step)
+(
+    f_prime_approx,
+    number_of_function_evaluations,
+) = method.compute_derivative(function, x)
+first_derivative_value = first_derivative(x)
+absolute_error = abs(f_prime_approx - first_derivative_value)
+print(
+    "x = %.3f, error = %.3e, Func. eval. = %d"
+    % (x, absolute_error, number_of_function_evaluations)
+)
+
+# %%
+# Perform the benchmark
+# ---------------------
+
+
+# %%
+print("+ Benchmark on several points")
+number_of_test_points = 21  # This number of test points is odd
+problem = nd.SinProblem()
+print(problem)
+interval = problem.get_interval()
+function = problem.get_function()
+first_derivative = problem.get_first_derivative()
+initial_step = 1.0e2
+test_points = np.linspace(interval[0], interval[1], number_of_test_points)
+relative_precision = 1.0e-15
+absolute_precision = abs(function(x)) * relative_precision
+method = ShiXieXuanNocedalGeneralMethod(absolute_precision, initial_step)
+average_relative_error, average_feval, data = nd.benchmark_method(
+    function,
+    first_derivative,
+    test_points,
+    method.compute_derivative,
+    verbose=True,
+)
+print("Average relative error =", average_relative_error)
+print("Average number of function evaluations =", average_feval)
+tabulate.tabulate(data, headers=["x", "Rel. err.", "F. Eval."], tablefmt="html")
+
+# %%
+# Notice that the method does not perform correctly for the point
+# :math:`x = 0` for the polynomial problem.
+# This test point appears only if the number of test points is odd,
+# because the test interval is symmetric with respect to :math:`x = 0`.
+# For this problem, the method does not perform correctly
+# because the value of the function is zero at :math:`x = 0`.
+# The method can perform correctly in this case, if it is
+# provided a consistent value of the absolute error of the function
+# value.
+# Here, we compute the absolute error depending on the relative
+# error and the absolute value of the value of the function.
+# If the value of the function is zero, then the computed
+# absolute error is zero, which produces a failure of the method.
+
+# %%
+# Map from the problem name to the initial step.
+
+# %%
+initial_step_map = {
+    "polynomial": 1.0,
+    "inverse": 1.0e0,
+    "exp": 1.0e-1,
+    "log": 1.0e-3,  # x > 0
+    "sqrt": 1.0e-3,  # x > 0
+    "atan": 1.0e0,
+    "sin": 1.0e0,
+    "scaled exp": 1.0e5,
+    "GMSW": 1.0e0,
+    "SXXN1": 1.0e0,
+    "SXXN2": 1.0e0,
+    "SXXN3": 1.0e0,
+    "SXXN4": 1.0e0,
+    "Oliver1": 1.0e0,
+    "Oliver2": 1.0e0,
+    "Oliver3": 1.0e-3,
+}
+
+# %%
+# The next script evaluates a collection of benchmark problems
+# using the :class:`~numericalderivative.ShiXieXuanNocedalGeneral` class.
+
+# %%
+number_of_test_points = 100  # This value can significantly change the results
+data = []
+function_list = nd.build_benchmark()
+number_of_functions = len(function_list)
+average_absolute_error_list = []
+average_feval_list = []
+relative_precision = 1.0e-15
+delta_x = 1.0e-9
+for i in range(number_of_functions):
+    problem = function_list[i]
+    name = problem.get_name()
+    initial_step = initial_step_map[name]
+    function = problem.get_function()
+    first_derivative = problem.get_first_derivative()
+    interval = problem.get_interval()
+    lower_x_bound, upper_x_bound = problem.get_interval()
+    print(f"Function #{i}, {name}")
+    if name == "sin":
+        # Change the lower and upper bound so that the points +/-pi
+        # are excluded (see below for details).
+        lower_x_bound += delta_x
+        upper_x_bound -= delta_x
+    test_points = np.linspace(lower_x_bound, upper_x_bound, number_of_test_points)
+    method = ShiXieXuanNocedalGeneralMethod(relative_precision, initial_step)
+    average_relative_error, average_feval, _ = nd.benchmark_method(
+        function, first_derivative, test_points, method.compute_derivative
+    )
+    average_absolute_error_list.append(average_relative_error)
+    average_feval_list.append(average_feval)
+    data.append(
+        (
+            name,
+            average_relative_error,
+            average_feval,
+        )
+    )
+data.append(
+    [
+        "Average",
+        np.nanmean(average_absolute_error_list),
+        np.nanmean(average_feval_list),
+    ]
+)
+tabulate.tabulate(
+    data,
+    headers=["Name", "Average rel. error", "Average func. eval"],
+    tablefmt="html",
+)
+
+# %%
+# Notice that the method cannot perform correctly for the sin function
+# at the point :math:`x = \pm \pi`.
+# Indeed, this function is such that :math:`f''(x) = 0` if :math:`x = \pm \pi`.
+# In this case, the test ratio is zero and the method cannot work.
+# Therefore, we make so that the points :math:`\pm \pi` are excluded from the benchmark.
+# The same problem appears at the point :math:`x = 0`.
+# This point is not included in the test set if the number of points is even
+# (e.g. with `number_of_test_points = 100`), but it might appear if the
+# number of test points is odd.
+
+# %%

doc/examples/run_examples.sh

@@ -15,3 +15,4 @@ python3 plot_stepleman_winarsky.py
 python3 plot_shi_xie_xuan_nocedal_forward.py
 python3 plot_shi_xie_xuan_nocedal_forward_benchmark.py
 python3 plot_shi_xie_xuan_nocedal_general.py
+python3 plot_shi_xie_xuan_nocedal_general_benchmark.py