major rewrite representing the name changes, mapping to CPython scipy…

….integrate, more background info

docs/scipy-integrate.ipynb

@@ -225,26 +225,25 @@
    "source": [
     "# scipy.integrate\n",
     "\n",
-    "This module defines four different numeric integration algorithms:\n",
+    "This module provides a simplified subset of CPython's `scipy.integrate` module. The algorithms were not ported from CPython's `scipy.integrate` for the sake of resource usage, but derived from a paper found in https://www.genivia.com/qthsh.html. There are four numerical integration algorithms:\n",
     "\n",
     "1. [scipy.integrate.quad](#quad)\n",
     "2. [scipy.integrate.romberg](#romberg)\n",
     "3. [scipy.integrate.simpson](#simpson)\n",
-    "4. [scipy.integrate.quadgk](#quadgk)\n",
-    "\n"
+    "4. [scipy.integrate.tanhsinh](#tanhsinh)\n"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Motivation\n",
+    "## Introduction\n",
     "\n",
-    "Why would we need multiple integration algorithms at all? The main reason is that all of these have their individual uses in numerical math. quad for example can handle integrands that have singularities which simpson cannot, while simpson is good enough for simple polynomials, if you choose the number of steps carefully. \n",
+    "Numerical integration works best with float64 math enabled. If you require float64 math, be sure to set `MICROPY_OBJ_REPR_A` and `MICROPY_FLOAT_IMPL_DOUBLE`. This being said, the modules work equally well using float32, albeit with reduced precision. The required error tolerance can be specified for each of the function calls using the \"eps=\" option, defaulting to the compiled in `etolerance` value (1e-14 for fp64, 1e-8 for fp32).\n",
     "\n",
-    "In any case, before you use numerical integration, be sure you understand the algorithms' individual strengths, weaknesses, and limits. A good starting point is the paper found in https://www.genivia.com/qthsh.html, from which all of these algorithms were ported. \n",
+    "The submodule can be enabled by setting `ULAB_SCIPY_HAS_INTEGRATE_MODULE` in `code/ulab.h`. As for the individual integration algorithms, you can select which to include by setting one or more of `ULAB_INTEGRATE_HAS_QUAD`, `ULAB_INTEGRATE_HAS_ROMBERG`, `ULAB_INTEGRATE_HAS_SIMPSON`, and `ULAB_INTEGRATE_HAS_TANHSINH`.\n",
     "\n",
-    "Also note that these algorithms do not support complex numbers. "
+    "Also note that these algorithms do not support complex numbers, although it is certainly possible to implement complex integration in MicroPython on top of this module, e.g. as in https://stackoverflow.com/questions/5965583/use-scipy-integrate-quad-to-integrate-complex-numbers.  "
    ]
   },
   {
@@ -255,16 +254,16 @@
     "\n",
     "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html \n",
     "\n",
-    "Implements an optimized Tanh-Sinh, Sinh-Sinh and Exp-Sinh quadrature algorithm. It is especially applied where singularities or infinite derivatives exist at one or both endpoints. The method uses hyperbolic functions in a change of variables to transform an integral on the interval x ∈ (−1, 1) to an integral on the entire real line t ∈ (−∞, ∞), the two integrals having the same value. After this transformation, the integrand decays with a double exponential rate, and thus, this method is also known as the double exponential (DE) formula (https://en.wikipedia.org/wiki/Tanh-sinh_quadrature). \n",
+    "In CPython `scipy.integrate`, `quad` is a wrapper implementing many algorithms based on the Fortran QUADPACK package. Gauss-Kronrod is just one of them, and it is useful for most general-purpose tasks. This particular function implements an Adaptive Gauss-Kronrod (G10,K21) quadrature algorithm. The Gauss–Kronrod quadrature formula is a variant of Gaussian quadrature, in which the evaluation points are chosen so that an accurate approximation can be computed by re-using the information produced by the computation of a less accurate approximation (https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula). \n",
     "\n",
-    "The function takes three to five positional arguments: \n",
+    "The function takes three to five arguments: \n",
     "\n",
-    "* f, a function defined e.g. with lambda,\n",
+    "* f, a callable,\n",
     "* a and b, the lower and upper integration limit, \n",
-    "* levels=, the number of loops taken to calculate (default 6)\n",
-    "* eps=, the error tolerance (default 1e-14) \n",
+    "* order=, the order of integration (default 5),\n",
+    "* eps=, the error tolerance (default etolerance) \n",
     "\n",
-    "The function returns the result and the final error tolerance as a tuple of mp_float_t. "
+    "The function returns the result and the error estimate as a tuple of floats.  "
    ]
   },
   {
@@ -288,10 +287,10 @@
    "source": [
     "%%micropython -unix 1\n",
     "\n",
-    "from ulab import scipy as sc\n",
+    "from ulab import scipy\n",
     "\n",
     "f = lambda x: x**2 + 2*x + 1\n",
-    "result = sc.integrate.quad(f, 0, 5, levels=6, eps=1e-14)\n",
+    "result = scipy.integrate.quad(f, 0, 5, order=5, eps=1e-10)\n",
     "print (f\"result = {result}\")\n"
    ]
   },
@@ -303,16 +302,16 @@
     "\n",
     "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romberg.html \n",
     "\n",
-    "Implements the Romberg quadrature algorithm.  Romberg's method is a Newton–Cotes formula – it evaluates the integrand at equally spaced points. The integrand must have continuous derivatives, though fairly good results may be obtained if only a few derivatives exist. If it is possible to evaluate the integrand at unequally spaced points, then other methods such as Gaussian quadrature and Clenshaw–Curtis quadrature are generally more accurate (https://en.wikipedia.org/wiki/Romberg%27s_method). \n",
+    "This function implements the Romberg quadrature algorithm.  Romberg's method is a Newton–Cotes formula – it evaluates the integrand at equally spaced points. The integrand must have continuous derivatives, though fairly good results may be obtained if only a few derivatives exist. If it is possible to evaluate the integrand at unequally spaced points, then other methods such as Gaussian quadrature and Clenshaw–Curtis quadrature are generally more accurate (https://en.wikipedia.org/wiki/Romberg%27s_method). \n",
     "\n",
-    "Please note: This function is deprecated as of SciPy 1.12.0 and will be removed in SciPy 1.15.0. Please use scipy.integrate.quad instead. \n",
+    "Please note: This function is deprecated as of SciPy 1.12.0 and will be removed in SciPy 1.15.0. Please use `scipy.integrate.quad` instead. \n",
     "\n",
-    "The function takes three to five positional arguments: \n",
+    "The function takes three to five arguments: \n",
     "\n",
-    "* f, a function defined e.g. with lambda,\n",
+    "* f, a callable,\n",
     "* a and b, the lower and upper integration limit, \n",
-    "* steps=, the number of steps taken to calculate (default 100)\n",
-    "* eps=, the error tolerance (default 1e-14) \n",
+    "* steps=, the number of steps taken to calculate (default 100),\n",
+    "* eps=, the error tolerance (default etolerance) \n",
     "\n",
     "The function returns the result as a float.\n"
    ]
@@ -333,10 +332,10 @@
    "source": [
     "%%micropython -unix 1\n",
     "\n",
-    "from ulab import scipy as sc\n",
+    "from ulab import scipy\n",
     "\n",
     "f = lambda x: x**2 + 2*x + 1\n",
-    "result = sc.integrate.romberg(f, 0, 5)\n",
+    "result = scipy.integrate.romberg(f, 0, 5)\n",
     "print (f\"result = {result}\")\n"
    ]
   },
@@ -348,78 +347,96 @@
     "\n",
     "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.simpson.html \n",
     "\n",
-    "Implements an Adaptive Simpson quadrature algorithm.  Adaptive Simpson's method, also called adaptive Simpson's rule, is a method of numerical integration proposed by G.F. Kuncir in 1962. It is probably the first recursive adaptive algorithm for numerical integration to appear in print, although more modern adaptive methods based on Gauss–Kronrod quadrature and Clenshaw–Curtis quadrature are now generally preferred (https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method). \n",
+    "This function is different from CPython's `simpson` method in that it does not take an array of function values but determines the optimal spacing of samples itself. Adaptive Simpson's method, also called adaptive Simpson's rule, is a method of numerical integration proposed by G.F. Kuncir in 1962. It is probably the first recursive adaptive algorithm for numerical integration to appear in print, although more modern adaptive methods based on Gauss–Kronrod quadrature and Clenshaw–Curtis quadrature are now generally preferred (https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method). \n",
     "\n",
-    "The function takes three to five positional arguments: \n",
+    "The function takes three to five arguments: \n",
     "\n",
-    "* f, a function defined e.g. with lambda,\n",
+    "* f, a callable,\n",
     "* a and b, the lower and upper integration limit, \n",
-    "* steps=, the number of steps taken to calculate (default 100)\n",
-    "* eps=, the error tolerance (default 1e-14) \n",
+    "* steps=, the number of steps taken to calculate (default 100),\n",
+    "* eps=, the error tolerance (default etolerance) \n",
     "\n",
-    "The function returns the result as a mp_float_t."
+    "The function returns the result as a float."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "UsageError: Cell magic `%%micropython` not found.\n"
+     ]
+    }
+   ],
    "source": [
     "%%micropython -unix 1\n",
     "\n",
-    "from ulab import scipy as sc\n",
+    "from ulab import scipy\n",
     "\n",
     "f = lambda x: x**2 + 2*x + 1\n",
-    "result = sc.integrate.simpson(f, 0, 5)\n",
+    "result = scipy.integrate.simpson(f, 0, 5)\n",
     "print (f\"result = {result}\")\n"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## quadgk\n",
+    "## tanhsinh\n",
+    "\n",
+    "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html \n",
+    "\n",
+    "In CPython `scipy.integrate`, `tanhsinh` is written in Python (https://github.com/scipy/scipy/blob/main/scipy/integrate/_tanhsinh.py). It is used in cases where Newton-Cotes, Gauss-Kronrod, and other formulae do not work due to properties of the integrand or the integration limits. (In SciPy v1.14.1, it is not a public function but it has been marked as public in SciPy v1.15.0rc1). \n",
     "\n",
-    "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quadgk.html \n",
+    "This particular function implements an optimized Tanh-Sinh, Sinh-Sinh and Exp-Sinh quadrature algorithm. It is especially applied where singularities or infinite derivatives exist at one or both endpoints. The method uses hyperbolic functions in a change of variables to transform an integral on the interval x ∈ (−1, 1) to an integral on the entire real line t ∈ (−∞, ∞), the two integrals having the same value. After this transformation, the integrand decays with a double exponential rate, and thus, this method is also known as the double exponential (DE) formula (https://en.wikipedia.org/wiki/Tanh-sinh_quadrature). \n",
     "\n",
-    "Implements an Adaptive Gauss-Kronrod (G10,K21) quadrature algorithm. The Gauss–Kronrod quadrature formula is a variant of Gaussian quadrature, in which the evaluation points are chosen so that an accurate approximation can be computed by re-using the information produced by the computation of a less accurate approximation (https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula). \n",
+    "As opposed to the three algorithms mentioned before, it also supports integrals with infinite limits like the Gaussian integral (https://en.wikipedia.org/wiki/Gaussian_integral), as shown below. \n",
     "\n",
-    "The function takes three to five positional arguments: \n",
+    "The function takes three to five  arguments: \n",
     "\n",
-    "* f, a function defined e.g. with lambda,\n",
+    "* f, a callable,\n",
     "* a and b, the lower and upper integration limit, \n",
-    "* order=, the order of integration (default 5)\n",
-    "* eps=, the error tolerance (default 1e-14) \n",
+    "* levels=, the number of loops taken to calculate (default 6),\n",
+    "* eps=, the error tolerance (default: etolerance)\n",
     "\n",
-    "The function returns the result and the final error tolerance as a tuple of mp_float_t. "
+    "The function returns the result and the error estimate as a tuple of floats.\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 5,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "UsageError: Cell magic `%%micropython` not found.\n"
+     ]
+    }
+   ],
    "source": [
     "%%micropython -unix 1\n",
     "\n",
-    "from ulab import scipy as sc\n",
-    "\n",
-    "f = lambda x: x**2 + 2*x + 1\n",
-    "result = sc.integrate.quadgk(f, 0, 5)\n",
-    "print (f\"result = {result}\")\n"
+    "from ulab import scipy, numpy as np\n",
+    "from math import *\n",
+    "f = lambda x: exp(- x**2)\n",
+    "result = scipy.integrate.tanhsinh(f, -np.inf, np.inf)\n",
+    "print (f\"result = {result}\")\n",
+    "exact = sqrt(pi)   # which is the exact value\n",
+    "print (f\"exact value = {exact}\")\n"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
-   "source": [
-    "# Note to Package Builders\n",
-    "\n",
-    "Numerical integration makes little sense with float32 math. For compilation, please be sure to set `MICROPY_OBJ_REPR_A` and `MICROPY_FLOAT_IMPL_DOUBLE`, otherwise compilation may fail. \n",
-    "\n",
-    "The submodule can be enabled by setting `ULAB_SCIPY_HAS_INTEGRATE_MODULE`. As for the individual integration algorithms, you can select which to include by setting one or more of `ULAB_INTEGRATE_HAS_QUAD`, `ULAB_INTEGRATE_HAS_ROMBERG`, `ULAB_INTEGRATE_HAS_SIMPSON`, and `ULAB_INTEGRATE_HAS_QUADGK`. The latter will only be built with `MICROPY_FLOAT_IMPL_DOUBLE` enabled. \n"
-   ]
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {