From 73ed8cc11f57c5c7605d97eaf5303687a2d21106 Mon Sep 17 00:00:00 2001 From: Harald Milz <149242962+h-milz@users.noreply.github.com> Date: Sun, 15 Dec 2024 18:49:08 +0100 Subject: [PATCH] Add scipy integration (#699) * adding scipy integrate, initial check-in * compile unix double-precision, select integrations algos * bumping ulab version number to 6.7.0 * adding documentation * documentation fix * documentation fix * documentation fix * rewritten in some places * complex number error handling * added test cases * resolved importing scipy.integrate * resolved importing scipy.integrate #2 * build integrate only when we have MICROPY_FLOAT_IMPL_DOUBLE * reverting commit a4c0c0b * re-pushing failed commit * Revert "re-pushing failed commit" This reverts commit a10e89fe14ce17306bb35e629509c01b1b6f8094. * improve tests using math.isclose() * enabled fp32 builds * removed conditional includes * adapted to new function names, corrected importing * function names similar to in CPython scipy.integrate, some minor corrections * major rewrite representing the name changes, mapping to CPython scipy.integrate, more background info --- code/micropython.mk | 1 + code/scipy/integrate/integrate.c | 701 +++++++++++++++++++++++++++++++ code/scipy/integrate/integrate.h | 34 ++ code/scipy/scipy.c | 6 +- code/ulab.c | 2 +- code/ulab.h | 22 + docs/scipy-integrate.ipynb | 510 ++++++++++++++++++++++ tests/2d/scipy/integrate.py | 28 ++ tests/2d/scipy/integrate.py.exp | 4 + 9 files changed, 1306 insertions(+), 2 deletions(-) create mode 100644 code/scipy/integrate/integrate.c create mode 100644 code/scipy/integrate/integrate.h create mode 100644 docs/scipy-integrate.ipynb create mode 100644 tests/2d/scipy/integrate.py create mode 100644 tests/2d/scipy/integrate.py.exp diff --git a/code/micropython.mk b/code/micropython.mk index 4aa6f615..e835d87b 100644 --- a/code/micropython.mk +++ b/code/micropython.mk @@ -2,6 +2,7 @@ USERMODULES_DIR := $(USERMOD_DIR) # Add all C files to SRC_USERMOD. +SRC_USERMOD += $(USERMODULES_DIR)/scipy/integrate/integrate.c SRC_USERMOD += $(USERMODULES_DIR)/scipy/linalg/linalg.c SRC_USERMOD += $(USERMODULES_DIR)/scipy/optimize/optimize.c SRC_USERMOD += $(USERMODULES_DIR)/scipy/signal/signal.c diff --git a/code/scipy/integrate/integrate.c b/code/scipy/integrate/integrate.c new file mode 100644 index 00000000..80def711 --- /dev/null +++ b/code/scipy/integrate/integrate.c @@ -0,0 +1,701 @@ +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2024 Harald Milz + * + * References: + * - Dr. Robert van Engelen, Improving the mp_float_t Exponential Quadrature Tanh-Sinh, Sinh-Sinh and Exp-Sinh Formulas, + * 2021, https://www.genivia.com/qthsh.html + * - Borwein, Bailey & Girgensohn, "Experimentation in Mathematics - Computational Paths to Discovery", A K Peters, + * 2003, pages 312-313 + * - Joren Vanherck, Bart Sorée, Wim Magnus, Tanh-sinh quadrature for single and multiple integration using + * floating-point arithmetic, 2020, https://arxiv.org/abs/2007.15057 + * - Tanh-Sinh quadrature, Wikipedia, https://en.wikipedia.org/wiki/Tanh-sinh_quadrature + * - Romberg's method, Wikipedia, https://en.wikipedia.org/wiki/Romberg%27s_method + * - Adaptive Simpson's method, Wikipedia, https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method + * - Gauss–Kronrod quadrature formula, Wikipedia, https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula + * + * This module provides four integration methods, and thus deviates from scipy.integrate a bit. + * As for the pros and cons of the different methods please consult the literature above. + * The code was ported to Micropython from Dr. Engelen's paper and used with his written kind permission + * - quad - Tanh-Sinh, Sinh-Sinh and Exp-Sinh quadrature + * - romberg - Romberg quadrature + * - simpson - Adaptive Simpson quadrature + * - quadgk - Adaptive Gauss-Kronrod (G10,K21) quadrature + */ + +#include +#include "py/obj.h" +#include "py/runtime.h" +#include "py/misc.h" +#include "py/objtuple.h" + +#include "../../ndarray.h" +#include "../../ulab.h" +#include "../../ulab_tools.h" +#include "integrate.h" + +#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE +ULAB_DEFINE_FLOAT_CONST(etolerance, MICROPY_FLOAT_CONST(1e-14), 0x283424dcUL, 0x3e901b2b29a4692bULL); +#define MACHEPS MICROPY_FLOAT_CONST(1e-17) +#else +ULAB_DEFINE_FLOAT_CONST(etolerance, MICROPY_FLOAT_CONST(1e-8), 0x358637cfUL, 0x3e7010c6f7d42d18ULL); +#define MACHEPS MICROPY_FLOAT_CONST(1e-8) +#endif + +#define ZERO MICROPY_FLOAT_CONST(0.0) +#define POINT_TWO_FIVE MICROPY_FLOAT_CONST(0.25) +#define ONE MICROPY_FLOAT_CONST(1.0) +#define TWO MICROPY_FLOAT_CONST(2.0) +#define FOUR MICROPY_FLOAT_CONST(4.0) +#define SIX MICROPY_FLOAT_CONST(6.0) +#define TEN MICROPY_FLOAT_CONST(10.0) +#define FIFTEEN MICROPY_FLOAT_CONST(15.0) +#define EPS_5 MICROPY_FLOAT_CONST(1e-5) + + +static mp_float_t integrate_python_call(const mp_obj_type_t *type, mp_obj_t fun, mp_float_t x, mp_obj_t *fargs, uint8_t nparams) { + // Helper function for calculating the value of f(x, a, b, c, ...), + // where f is defined in python. Takes a float, returns a float. + // The array of mp_obj_t type must be supplied, as must the number of parameters (a, b, c...) in nparams + fargs[0] = mp_obj_new_float(x); + return mp_obj_get_float(MP_OBJ_TYPE_GET_SLOT(type, call)(fun, nparams+1, 0, fargs)); +} + +// sign helper function +int sign(mp_float_t x) { + if (x >= ZERO) + return 1; + else + return -1; +} + + +#if ULAB_INTEGRATE_HAS_TANHSINH +// Tanh-Sinh, Sinh-Sinh and Exp-Sinh quadrature +// https://www.genivia.com/qthsh.html + +// return optimized Exp-Sinh integral split point d +mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t eps, mp_float_t d) { + const mp_obj_type_t *type = mp_obj_get_type(fun); + mp_obj_t fargs[1]; + mp_float_t h2 = integrate_python_call(type, fun, a + d/2, fargs, 0) - integrate_python_call(type, fun, (a + d*2)*4, fargs, 0); + int i = 1, j = 32; // j=32 is optimal to find r + if (isfinite(h2) && MICROPY_FLOAT_C_FUN(fabs)(h2) > EPS_5) { // if |h2| > 2^-16 + mp_float_t r, fl, fr, h, s = 0, lfl, lfr, lr = 2; + do { // find max j such that fl and fr are finite + j /= 2; + r = 1 << (i + j); + fl = integrate_python_call(type, fun, a + d/r, fargs, 0); + fr = integrate_python_call(type, fun, (a + d*r)*r*r, fargs, 0); + h = fl - fr; + } while (j > 1 && !isfinite(h)); + if (j > 1 && isfinite(h) && sign(h) != sign(h2)) { + lfl = fl; // last fl=f(a+d/r) + lfr = fr; // last fr=f(a+d*r)*r*r + do { // bisect in 4 iterations + j /= 2; + r = 1 << (i + j); + fl = integrate_python_call(type, fun, a + d/r, fargs, 0); + fr = integrate_python_call(type, fun, (a + d*r)*r*r, fargs, 0); + h = fl - fr; + if (isfinite(h)) { + s += MICROPY_FLOAT_C_FUN(fabs)(h); // sum |h| to remove noisy cases + if (sign(h) == sign(h2)) { + i += j; // search right half + } + else { // search left half + lfl = fl; // record last fl=f(a+d/r) + lfr = fr; // record last fl=f(a+d*r)*r*r + lr = r; // record last r + } + } + } while (j > 1); + if (s > eps) { // if sum of |h| > eps + h = lfl - lfr; // use last fl and fr before the sign change + r = lr; // use last r before the sign change + if (h != ZERO) // if last diff != 0, back up r by one step + r /= TWO; + if (MICROPY_FLOAT_C_FUN(fabs)(lfl) < MICROPY_FLOAT_C_FUN(fabs)(lfr)) + d /= r; // move d closer to the finite endpoint + else + d *= r; // move d closer to the infinite endpoint + } + } + } + return d; +} + + +// integrate function f, range a..b, max levels n, error tolerance eps +mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint16_t n, mp_float_t eps, mp_float_t *e) { + const mp_obj_type_t *type = mp_obj_get_type(fun); + mp_obj_t fargs[1]; + const mp_float_t tol = TEN*eps; + mp_float_t c = ZERO, d = ONE, s, sign = ONE, v, h = TWO; + int k = 0, mode = 0; // Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2 + if (b < a) { // swap bounds + v = b; + b = a; + a = v; + sign = -1; + } + if (isfinite(a) && isfinite(b)) { + c = (a+b)/TWO; + d = (b-a)/TWO; + v = c; + } + else if (isfinite(a)) { + mode = 1; // Exp-Sinh + d = exp_sinh_opt_d(fun, a, eps, d); + c = a; + v = a+d; + } + else if (isfinite(b)) { + mode = 1; // Exp-Sinh + // d = -d; + d = exp_sinh_opt_d(fun, b, eps, -d); + sign = -sign; + c = b; + v = b+d; + } + else { + mode = 2; // Sinh-Sinh + v = ZERO; + } + s = integrate_python_call(type, fun, v, fargs, 0); + do { + mp_float_t p = ZERO, q, fp = ZERO, fm = ZERO, t, eh; + h /= TWO; + t = eh = MICROPY_FLOAT_C_FUN(exp)(h); + if (k > ZERO) + eh *= eh; + if (mode == 0) { // Tanh-Sinh + do { + mp_float_t u = MICROPY_FLOAT_C_FUN(exp)(ONE / t - t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2 + mp_float_t r = TWO * u / (ONE + u); // = 1 - tanh(sinh(j*h)) + mp_float_t w = (t + ONE / t) * r / (ONE + u); // = cosh(j*h)/cosh(sinh(j*h))^2 + mp_float_t x = d*r; + if (a+x > a) { // if too close to a then reuse previous fp + mp_float_t y = integrate_python_call(type, fun, a+x, fargs, 0); + if (isfinite(y)) + fp = y; // if f(x) is finite, add to local sum + } + if (b-x < b) { // if too close to a then reuse previous fp + mp_float_t y = integrate_python_call(type, fun, b-x, fargs, 0); + if (isfinite(y)) + fm = y; // if f(x) is finite, add to local sum + } + q = w*(fp+fm); + p += q; + t *= eh; + } while (MICROPY_FLOAT_C_FUN(fabs)(q) > eps*MICROPY_FLOAT_C_FUN(fabs)(p)); + } + else { + t /= TWO; + do { + mp_float_t r = MICROPY_FLOAT_C_FUN(exp)(t - POINT_TWO_FIVE / t); // = exp(sinh(j*h)) + mp_float_t x, y, w = r; + q = ZERO; + if (mode == 1) { // Exp-Sinh + x = c + d/r; + if (x == c) // if x hit the finite endpoint then break + break; + y = integrate_python_call(type, fun, x, fargs, 0); + if (isfinite(y)) // if f(x) is finite, add to local sum + q += y/w; + } + else { // Sinh-Sinh + r = (r - ONE / r) / TWO; // = sinh(sinh(j*h)) + w = (w + ONE / w) / TWO; // = cosh(sinh(j*h)) + x = c - d*r; + y = integrate_python_call(type, fun, x, fargs, 0); + if (isfinite(y)) // if f(x) is finite, add to local sum + q += y*w; + } + x = c + d*r; + y = integrate_python_call(type, fun, x, fargs, 0); + if (isfinite(y)) // if f(x) is finite, add to local sum + q += y*w; + q *= t + POINT_TWO_FIVE / t; // q *= cosh(j*h) + p += q; + t *= eh; + } while (MICROPY_FLOAT_C_FUN(fabs)(q) > eps*MICROPY_FLOAT_C_FUN(fabs)(p)); + } + v = s-p; + s += p; + ++k; + } while (MICROPY_FLOAT_C_FUN(fabs)(v) > tol*MICROPY_FLOAT_C_FUN(fabs)(s) && k <= n); + // return the error estimate by reference + *e = MICROPY_FLOAT_C_FUN(fabs)(v)/(MICROPY_FLOAT_C_FUN(fabs)(s)+eps); + return sign*d*s*h; // result with estimated relative error e +} + +//| def tanhsinh( +//| fun: Callable[[float], float], +//| a: float, +//| b: float, +//| *, +//| levels: int = 6 +//| eps: float = etolerance +//| ) -> float: +//| """ +//| :param callable f: The function to integrate +//| :param float a: The lower integration limit +//| :param float b: The upper integration limit +//| :param float levels: The number of levels to perform (6..7 is optimal) +//| :param float eps: The error tolerance value +//| +//| Find a quadrature of the function ``f(x)`` on the interval +//| (``a``..``b``) using an optimized double exponential. The result is accurate to within +//| ``eps`` unless more than ``levels`` levels are required.""" +//| + + +static mp_obj_t integrate_tanhsinh(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_levels, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 6} }, + { MP_QSTR_eps, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = ULAB_REFERENCE_FLOAT_CONST(etolerance)} }, + }; + + mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; + mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); + + mp_obj_t fun = args[0].u_obj; + const mp_obj_type_t *type = mp_obj_get_type(fun); + if(!MP_OBJ_TYPE_HAS_SLOT(type, call)) { + mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a callable")); + } + + // iterate over args 1, 2, and 4 + // arg 3 will be handled by MP_ARG_INT above. + for (int i=1; i<=4; i*=2) { + type = mp_obj_get_type(args[i].u_obj); + if (type != &mp_type_float && type != &mp_type_int) { + mp_raise_msg_varg(&mp_type_TypeError, + MP_ERROR_TEXT("can't convert arg %d from %s to float"), i, mp_obj_get_type_str(args[i].u_obj)); + } + } + mp_float_t a = mp_obj_get_float(args[1].u_obj); + mp_float_t b = mp_obj_get_float(args[2].u_obj); + uint16_t n = (uint16_t)args[3].u_int; + if (n < 1) { + mp_raise_ValueError(MP_ERROR_TEXT("levels needs to be a positive integer")); + } + mp_float_t eps = mp_obj_get_float(args[4].u_obj); + + mp_obj_t res[2]; + mp_float_t e; + res[0] = mp_obj_new_float(tanhsinh(fun, a, b, n, eps, &e)); + res[1] = mp_obj_new_float(e); + return mp_obj_new_tuple(2, res); +} + +MP_DEFINE_CONST_FUN_OBJ_KW(integrate_tanhsinh_obj, 2, integrate_tanhsinh); +#endif /* ULAB_INTEGRATE_HAS_TANHSINH */ + +#if ULAB_INTEGRATE_HAS_ROMBERG +// Romberg quadrature +// 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. +// https://en.wikipedia.org/wiki/Romberg%27s_method, https://www.genivia.com/qthsh.html, +// https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romberg.html (which is different +// insofar as the latter expects an array of function values). + +mp_float_t qromb(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint16_t n, mp_float_t eps) { + const mp_obj_type_t *type = mp_obj_get_type(fun); + mp_obj_t fargs[1]; + mp_float_t R1[n], R2[n]; + mp_float_t *Ro = &R1[0], *Ru = &R2[0]; + mp_float_t h = b-a; + uint16_t i, j; + Ro[0] = (integrate_python_call(type, fun, a, fargs, 0) + integrate_python_call(type, fun, b, fargs, 0)) * h/2; + for (i = 1; i < n; ++i) { + unsigned long long k = 1UL << i; + unsigned long long s = 1; + mp_float_t sum = ZERO; + mp_float_t *Rt; + h /= TWO; + for (j = 1; j < k; j += 2) + sum += integrate_python_call(type, fun, a+j*h, fargs, 0); + Ru[0] = h*sum + Ro[0] / TWO; + for (j = 1; j <= i; ++j) { + s <<= 2; + Ru[j] = (s*Ru[j-1] - Ro[j-1])/(s-1); + } + if (i > 2 && MICROPY_FLOAT_C_FUN(fabs)(Ro[i-1]-Ru[i]) <= eps*MICROPY_FLOAT_C_FUN(fabs)(Ru[i])+eps) + return Ru[i]; + Rt = Ro; + Ro = Ru; + Ru = Rt; + } + return Ro[n-1]; +} + +//| def romberg( +//| fun: Callable[[float], float], +//| a: float, +//| b: float, +//| *, +//| steps: int = 100 +//| eps: float = etolerance +//| ) -> float: +//| """ +//| :param callable f: The function to integrate +//| :param float a: The lower integration limit +//| :param float b: The upper integration limit +//| :param float steps: The number of equidistant steps +//| :param float eps: The tolerance value +//| +//| Find a quadrature of the function ``f(x)`` on the interval +//| (``a``..``b``) using the Romberg method. The result is accurate to within +//| ``eps`` unless more than ``steps`` steps are required.""" +//| ... +//| + +static mp_obj_t integrate_romberg(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_steps, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 100} }, + { MP_QSTR_eps, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = ULAB_REFERENCE_FLOAT_CONST(etolerance)} }, + }; + + mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; + mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); + + mp_obj_t fun = args[0].u_obj; + const mp_obj_type_t *type = mp_obj_get_type(fun); + if(!MP_OBJ_TYPE_HAS_SLOT(type, call)) { + mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a callable")); + } + + // iterate over args 1, 2, and 4 + // arg 3 will be handled by MP_ARG_INT above. + for (int i=1; i<=4; i*=2) { + type = mp_obj_get_type(args[i].u_obj); + if (type != &mp_type_float && type != &mp_type_int) { + mp_raise_msg_varg(&mp_type_TypeError, + MP_ERROR_TEXT("can't convert arg %d from %s to float"), i, mp_obj_get_type_str(args[i].u_obj)); + } + } + mp_float_t a = mp_obj_get_float(args[1].u_obj); + mp_float_t b = mp_obj_get_float(args[2].u_obj); + uint16_t steps = (uint16_t)args[3].u_int; + if (steps < 1) { + mp_raise_ValueError(MP_ERROR_TEXT("steps needs to be a positive integer")); + } + mp_float_t eps = mp_obj_get_float(args[4].u_obj); + + return mp_obj_new_float(qromb(fun, a, b, steps, eps)); +} + +MP_DEFINE_CONST_FUN_OBJ_KW(integrate_romberg_obj, 2, integrate_romberg); +#endif /* ULAB_INTEGRATE_HAS_ROMBERG */ + +#if ULAB_INTEGRATE_HAS_SIMPSON +// Adaptive Simpson quadrature +// https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method, https://www.genivia.com/qthsh.html + +mp_float_t as(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, mp_float_t fa, mp_float_t fm, + mp_float_t fb, mp_float_t v, mp_float_t eps, int n, mp_float_t t) { + const mp_obj_type_t *type = mp_obj_get_type(fun); + mp_obj_t fargs[1]; + mp_float_t h = (b-a) / TWO; + mp_float_t f1 = integrate_python_call(type, fun, a + h / TWO, fargs, 0); + mp_float_t f2 = integrate_python_call(type, fun, b - h / TWO, fargs, 0); + mp_float_t sl = h*(fa + FOUR * f1 + fm) / SIX; + mp_float_t sr = h*(fm + FOUR * f2 + fb) / SIX; + mp_float_t s = sl+sr; + mp_float_t d = (s-v) / FIFTEEN; + mp_float_t m = a+h; + if (n <= 0 || MICROPY_FLOAT_C_FUN(fabs)(d) < eps) + return t + s + d; // note: fabs(d) can be used as error estimate + eps /= TWO; + --n; + t = as(fun, a, m, fa, f1, fm, sl, eps, n, t); + return as(fun, m, b, fm, f2, fb, sr, eps, n, t); +} + +mp_float_t qasi(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n, mp_float_t eps) { + const mp_obj_type_t *type = mp_obj_get_type(fun); + mp_obj_t fargs[1]; + mp_float_t fa = integrate_python_call(type, fun, a, fargs, 0); + mp_float_t fm = integrate_python_call(type, fun, (a+b)/2, fargs, 0); + mp_float_t fb = integrate_python_call(type, fun, b, fargs, 0); + mp_float_t v = (fa + FOUR * fm + fb) * (b-a) / SIX; + return as(fun, a, b, fa, fm, fb, v, eps, n, 0); +} + +//| def simpson( +//| fun: Callable[[float], float], +//| a: float, +//| b: float, +//| *, +//| steps: int = 100 +//| eps: float = etolerance +//| ) -> float: +//| """ +//| :param callable f: The function to integrate +//| :param float a: The lower integration limit +//| :param float b: The upper integration limit +//| :param float steps: The number of equidistant steps +//| :param float eps: The tolerance value +//| +//| Find a quadrature of the function ``f(x)`` on the interval +//| (``a``..``b``) using the Adaptive Simpson's method. The result is accurate to within +//| ``eps`` unless more than ``steps`` steps are required.""" +//| ... +//| + +static mp_obj_t integrate_simpson(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_steps, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 100} }, + { MP_QSTR_eps, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = ULAB_REFERENCE_FLOAT_CONST(etolerance)} }, + }; + + mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; + mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); + + mp_obj_t fun = args[0].u_obj; + const mp_obj_type_t *type = mp_obj_get_type(fun); + if(!MP_OBJ_TYPE_HAS_SLOT(type, call)) { + mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a function")); + } + + // iterate over args 1, 2, and 4 + // arg 3 will be handled by MP_ARG_INT above. + for (int i=1; i<=4; i*=2) { + type = mp_obj_get_type(args[i].u_obj); + if (type != &mp_type_float && type != &mp_type_int) { + mp_raise_msg_varg(&mp_type_TypeError, + MP_ERROR_TEXT("can't convert arg %d from %s to float"), i, mp_obj_get_type_str(args[i].u_obj)); + } + } + mp_float_t a = mp_obj_get_float(args[1].u_obj); + mp_float_t b = mp_obj_get_float(args[2].u_obj); + uint16_t steps = (uint16_t)args[3].u_int; + if (steps < 1) { + mp_raise_ValueError(MP_ERROR_TEXT("steps needs to be a positive integer")); + } + mp_float_t eps = mp_obj_get_float(args[4].u_obj); + + return mp_obj_new_float(qasi(fun, a, b, steps, eps)); +} + +MP_DEFINE_CONST_FUN_OBJ_KW(integrate_simpson_obj, 2, integrate_simpson); +#endif /* ULAB_INTEGRATE_HAS_SIMPSON */ + +#if ULAB_INTEGRATE_HAS_QUAD +// Adaptive Gauss-Kronrod (G10,K21) quadrature +// https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula, https://www.genivia.com/qthsh.html + +mp_float_t gk(mp_float_t (*fun)(mp_float_t), mp_float_t c, mp_float_t d, mp_float_t *err) { +// abscissas and weights pre-calculated with Legendre Stieltjes polynomials + static const mp_float_t abscissas[21] = { + MICROPY_FLOAT_CONST(0.00000000000000000e+00), + MICROPY_FLOAT_CONST(7.65265211334973338e-02), + MICROPY_FLOAT_CONST(1.52605465240922676e-01), + MICROPY_FLOAT_CONST(2.27785851141645078e-01), + MICROPY_FLOAT_CONST(3.01627868114913004e-01), + MICROPY_FLOAT_CONST(3.73706088715419561e-01), + MICROPY_FLOAT_CONST(4.43593175238725103e-01), + MICROPY_FLOAT_CONST(5.10867001950827098e-01), + MICROPY_FLOAT_CONST(5.75140446819710315e-01), + MICROPY_FLOAT_CONST(6.36053680726515025e-01), + MICROPY_FLOAT_CONST(6.93237656334751385e-01), + MICROPY_FLOAT_CONST(7.46331906460150793e-01), + MICROPY_FLOAT_CONST(7.95041428837551198e-01), + MICROPY_FLOAT_CONST(8.39116971822218823e-01), + MICROPY_FLOAT_CONST(8.78276811252281976e-01), + MICROPY_FLOAT_CONST(9.12234428251325906e-01), + MICROPY_FLOAT_CONST(9.40822633831754754e-01), + MICROPY_FLOAT_CONST(9.63971927277913791e-01), + MICROPY_FLOAT_CONST(9.81507877450250259e-01), + MICROPY_FLOAT_CONST(9.93128599185094925e-01), + MICROPY_FLOAT_CONST(9.98859031588277664e-01), + }; + static const mp_float_t weights[21] = { + MICROPY_FLOAT_CONST(7.66007119179996564e-02), + MICROPY_FLOAT_CONST(7.63778676720807367e-02), + MICROPY_FLOAT_CONST(7.57044976845566747e-02), + MICROPY_FLOAT_CONST(7.45828754004991890e-02), + MICROPY_FLOAT_CONST(7.30306903327866675e-02), + MICROPY_FLOAT_CONST(7.10544235534440683e-02), + MICROPY_FLOAT_CONST(6.86486729285216193e-02), + MICROPY_FLOAT_CONST(6.58345971336184221e-02), + MICROPY_FLOAT_CONST(6.26532375547811680e-02), + MICROPY_FLOAT_CONST(5.91114008806395724e-02), + MICROPY_FLOAT_CONST(5.51951053482859947e-02), + MICROPY_FLOAT_CONST(5.09445739237286919e-02), + MICROPY_FLOAT_CONST(4.64348218674976747e-02), + MICROPY_FLOAT_CONST(4.16688733279736863e-02), + MICROPY_FLOAT_CONST(3.66001697582007980e-02), + MICROPY_FLOAT_CONST(3.12873067770327990e-02), + MICROPY_FLOAT_CONST(2.58821336049511588e-02), + MICROPY_FLOAT_CONST(2.03883734612665236e-02), + MICROPY_FLOAT_CONST(1.46261692569712530e-02), + MICROPY_FLOAT_CONST(8.60026985564294220e-03), + MICROPY_FLOAT_CONST(3.07358371852053150e-03), + }; + static const mp_float_t gauss_weights[10] = { + MICROPY_FLOAT_CONST(1.52753387130725851e-01), + MICROPY_FLOAT_CONST(1.49172986472603747e-01), + MICROPY_FLOAT_CONST(1.42096109318382051e-01), + MICROPY_FLOAT_CONST(1.31688638449176627e-01), + MICROPY_FLOAT_CONST(1.18194531961518417e-01), + MICROPY_FLOAT_CONST(1.01930119817240435e-01), + MICROPY_FLOAT_CONST(8.32767415767047487e-02), + MICROPY_FLOAT_CONST(6.26720483341090636e-02), + MICROPY_FLOAT_CONST(4.06014298003869413e-02), + MICROPY_FLOAT_CONST(1.76140071391521183e-02), + }; + const mp_obj_type_t *type = mp_obj_get_type(fun); + mp_obj_t fargs[1]; + mp_float_t p = ZERO; // kronrod quadrature sum + mp_float_t q = ZERO; // gauss quadrature sum + mp_float_t fp, fm; + mp_float_t e; + int i; + fp = integrate_python_call(type, fun, c, fargs, 0); + p = fp * weights[0]; + for (i = 1; i < 21; i += 2) { + fp = integrate_python_call(type, fun, c + d * abscissas[i], fargs, 0); + fm = integrate_python_call(type, fun, c - d * abscissas[i], fargs, 0); + p += (fp + fm) * weights[i]; + q += (fp + fm) * gauss_weights[i/2]; + } + for (i = 2; i < 21; i += 2) { + fp = integrate_python_call(type, fun, c + d * abscissas[i], fargs, 0); + fm = integrate_python_call(type, fun, c - d * abscissas[i], fargs, 0); + p += (fp + fm) * weights[i]; + } + *err = MICROPY_FLOAT_C_FUN(fabs)(p - q); + e = MICROPY_FLOAT_C_FUN(fabs)(2 * p * MACHEPS); // optional, to take 1e-17 MachEps prec. into account + if (*err < e) + *err = e; + return p; +} + +mp_float_t qakro(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n, mp_float_t tol, mp_float_t eps, mp_float_t *err) { + mp_float_t c = (a+b) / TWO; + mp_float_t d = (b-a) / TWO; + mp_float_t e; + mp_float_t r = gk(fun, c, d, &e); + mp_float_t s = d*r; + mp_float_t t = MICROPY_FLOAT_C_FUN(fabs)(s*tol); + if (tol == ZERO) + tol = t; + if (n > 0 && t < e && tol < e) { + s = qakro(fun, a, c, n-1, t / TWO, eps, err); + s += qakro(fun, c, b, n-1, t / TWO, eps, &e); + *err += e; + return s; + } + *err = e; + return s; +} + + +//| def quad( +//| fun: Callable[[float], float], +//| a: float, +//| b: float, +//| *, +//| order: int = 5 +//| eps: float = etolerance +//| ) -> float: +//| """ +//| :param callable f: The function to integrate +//| :param float a: The lower integration limit +//| :param float b: The upper integration limit +//| :param float order: Order of quadrature integration. Default is 5. +//| :param float eps: The tolerance value +//| +//| Find a quadrature of the function ``f(x)`` on the interval +//| (``a``..``b``) using the Adaptive Gauss-Kronrod method. The result is accurate to within +//| ``eps`` unless a higher order than ``order`` is required.""" +//| ... +//| + +static mp_obj_t integrate_quad(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_order, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 5} }, + { MP_QSTR_eps, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = ULAB_REFERENCE_FLOAT_CONST(etolerance)} }, + }; + + mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; + mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); + + mp_obj_t fun = args[0].u_obj; + const mp_obj_type_t *type = mp_obj_get_type(fun); + if(!MP_OBJ_TYPE_HAS_SLOT(type, call)) { + mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a callable")); + } + + // iterate over args 1, 2, and 4 + // arg 3 will be handled by MP_ARG_INT above. + for (int i=1; i<=4; i*=2) { + type = mp_obj_get_type(args[i].u_obj); + if (type != &mp_type_float && type != &mp_type_int) { + mp_raise_msg_varg(&mp_type_TypeError, + MP_ERROR_TEXT("can't convert arg %d from %s to float"), i, mp_obj_get_type_str(args[i].u_obj)); + } + } + mp_float_t a = mp_obj_get_float(args[1].u_obj); + mp_float_t b = mp_obj_get_float(args[2].u_obj); + uint16_t order = (uint16_t)args[3].u_int; + if (order < 1) { + mp_raise_ValueError(MP_ERROR_TEXT("order needs to be a positive integer")); + } + mp_float_t eps = mp_obj_get_float(args[4].u_obj); + + mp_obj_t res[2]; + mp_float_t e; + res[0] = mp_obj_new_float(qakro(fun, a, b, order, 0, eps, &e)); + res[1] = mp_obj_new_float(e); + return mp_obj_new_tuple(2, res); +} + +MP_DEFINE_CONST_FUN_OBJ_KW(integrate_quad_obj, 2, integrate_quad); +#endif /* ULAB_INTEGRATE_HAS_QUAD */ + +static const mp_rom_map_elem_t ulab_scipy_integrate_globals_table[] = { + { MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_integrate) }, +#if ULAB_INTEGRATE_HAS_TANHSINH + { MP_ROM_QSTR(MP_QSTR_tanhsinh), MP_ROM_PTR(&integrate_tanhsinh_obj) }, +#endif +#if ULAB_INTEGRATE_HAS_ROMBERG + { MP_ROM_QSTR(MP_QSTR_romberg), MP_ROM_PTR(&integrate_romberg_obj) }, +#endif +#if ULAB_INTEGRATE_HAS_SIMPSON + { MP_ROM_QSTR(MP_QSTR_simpson), MP_ROM_PTR(&integrate_simpson_obj) }, +#endif +#if ULAB_INTEGRATE_HAS_QUAD + { MP_ROM_QSTR(MP_QSTR_quad), MP_ROM_PTR(&integrate_quad_obj) }, +#endif +}; + +static MP_DEFINE_CONST_DICT(mp_module_ulab_scipy_integrate_globals, ulab_scipy_integrate_globals_table); + +const mp_obj_module_t ulab_scipy_integrate_module = { + .base = { &mp_type_module }, + .globals = (mp_obj_dict_t*)&mp_module_ulab_scipy_integrate_globals, +}; +#if CIRCUITPY_ULAB +MP_REGISTER_MODULE(MP_QSTR_ulab_dot_scipy_dot_integrate, ulab_scipy_integrate_module); +#endif + diff --git a/code/scipy/integrate/integrate.h b/code/scipy/integrate/integrate.h new file mode 100644 index 00000000..ebfac2ea --- /dev/null +++ b/code/scipy/integrate/integrate.h @@ -0,0 +1,34 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2024 Harald Milz + * +*/ + +#ifndef _SCIPY_INTEGRATE_ +#define _SCIPY_INTEGRATE_ + +#include "../../ulab_tools.h" + +extern const mp_obj_module_t ulab_scipy_integrate_module; + +#if ULAB_INTEGRATE_HAS_TANHSINH +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_tanhsinh_obj); +#endif +#if ULAB_INTEGRATE_HAS_ROMBERG +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_romberg_obj); +#endif +#if ULAB_INTEGRATE_HAS_SIMPSON +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_simpson_obj); +#endif +#if ULAB_INTEGRATE_HAS_QUAD +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_quad_obj); +#endif + +#endif /* _SCIPY_INTEGRATE_ */ + diff --git a/code/scipy/scipy.c b/code/scipy/scipy.c index 6895127b..e4ad306c 100644 --- a/code/scipy/scipy.c +++ b/code/scipy/scipy.c @@ -1,4 +1,3 @@ - /* * This file is part of the micropython-ulab project, * @@ -20,6 +19,8 @@ #include "signal/signal.h" #include "special/special.h" #include "linalg/linalg.h" +#include "integrate/integrate.h" + #if ULAB_HAS_SCIPY @@ -28,6 +29,9 @@ static const mp_rom_map_elem_t ulab_scipy_globals_table[] = { { MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_scipy) }, + #if ULAB_SCIPY_HAS_INTEGRATE_MODULE + { MP_ROM_QSTR(MP_QSTR_integrate), MP_ROM_PTR(&ulab_scipy_integrate_module) }, + #endif #if ULAB_SCIPY_HAS_LINALG_MODULE { MP_ROM_QSTR(MP_QSTR_linalg), MP_ROM_PTR(&ulab_scipy_linalg_module) }, #endif diff --git a/code/ulab.c b/code/ulab.c index 12f37027..be148c45 100644 --- a/code/ulab.c +++ b/code/ulab.c @@ -33,7 +33,7 @@ #include "user/user.h" #include "utils/utils.h" -#define ULAB_VERSION 6.6.1 +#define ULAB_VERSION 6.7.0 #define xstr(s) str(s) #define str(s) #s diff --git a/code/ulab.h b/code/ulab.h index 78b8a1ba..3eb30131 100644 --- a/code/ulab.h +++ b/code/ulab.h @@ -398,6 +398,28 @@ #define ULAB_NUMPY_HAS_WHERE (1) #endif +// the integrate module; functions of the integrate module still have +// to be defined separately +#ifndef ULAB_SCIPY_HAS_INTEGRATE_MODULE +#define ULAB_SCIPY_HAS_INTEGRATE_MODULE (1) +#endif + +#ifndef ULAB_INTEGRATE_HAS_TANHSINH +#define ULAB_INTEGRATE_HAS_TANHSINH (1) +#endif + +#ifndef ULAB_INTEGRATE_HAS_ROMBERG +#define ULAB_INTEGRATE_HAS_ROMBERG (1) +#endif + +#ifndef ULAB_INTEGRATE_HAS_SIMPSON +#define ULAB_INTEGRATE_HAS_SIMPSON (1) +#endif + +#ifndef ULAB_INTEGRATE_HAS_QUAD +#define ULAB_INTEGRATE_HAS_QUAD (1) +#endif + // the linalg module; functions of the linalg module still have // to be defined separately #ifndef ULAB_NUMPY_HAS_LINALG_MODULE diff --git a/docs/scipy-integrate.ipynb b/docs/scipy-integrate.ipynb new file mode 100644 index 00000000..23220231 --- /dev/null +++ b/docs/scipy-integrate.ipynb @@ -0,0 +1,510 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "ExecuteTime": { + "end_time": "2021-01-12T16:11:12.111639Z", + "start_time": "2021-01-12T16:11:11.914041Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Notebook magic" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2022-01-29T20:50:20.813162Z", + "start_time": "2022-01-29T20:50:20.794562Z" + } + }, + "outputs": [], + "source": [ + "from IPython.core.magic import Magics, magics_class, line_cell_magic\n", + "from IPython.core.magic import cell_magic, register_cell_magic, register_line_magic\n", + "from IPython.core.magic_arguments import argument, magic_arguments, parse_argstring\n", + "import subprocess\n", + "import os" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "ExecuteTime": { + "end_time": "2022-01-29T20:50:21.613220Z", + "start_time": "2022-01-29T20:50:21.557819Z" + } + }, + "outputs": [], + "source": [ + "@magics_class\n", + "class PyboardMagic(Magics):\n", + " @cell_magic\n", + " @magic_arguments()\n", + " @argument('-skip')\n", + " @argument('-unix')\n", + " @argument('-pyboard')\n", + " @argument('-file')\n", + " @argument('-data')\n", + " @argument('-time')\n", + " @argument('-memory')\n", + " def micropython(self, line='', cell=None):\n", + " args = parse_argstring(self.micropython, line)\n", + " if args.skip: # doesn't care about the cell's content\n", + " print('skipped execution')\n", + " return None # do not parse the rest\n", + " if args.unix: # tests the code on the unix port. Note that this works on unix only\n", + " with open('/dev/shm/micropython.py', 'w') as fout:\n", + " fout.write(cell)\n", + " proc = subprocess.Popen([\"../micropython/ports/unix/micropython-2\", \"/dev/shm/micropython.py\"], \n", + " stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n", + " print(proc.stdout.read().decode(\"utf-8\"))\n", + " print(proc.stderr.read().decode(\"utf-8\"))\n", + " return None\n", + " if args.file: # can be used to copy the cell content onto the pyboard's flash\n", + " spaces = \" \"\n", + " try:\n", + " with open(args.file, 'w') as fout:\n", + " fout.write(cell.replace('\\t', spaces))\n", + " printf('written cell to {}'.format(args.file))\n", + " except:\n", + " print('Failed to write to disc!')\n", + " return None # do not parse the rest\n", + " if args.data: # can be used to load data from the pyboard directly into kernel space\n", + " message = pyb.exec(cell)\n", + " if len(message) == 0:\n", + " print('pyboard >>>')\n", + " else:\n", + " print(message.decode('utf-8'))\n", + " # register new variable in user namespace\n", + " self.shell.user_ns[args.data] = string_to_matrix(message.decode(\"utf-8\"))\n", + " \n", + " if args.time: # measures the time of executions\n", + " pyb.exec('import utime')\n", + " message = pyb.exec('t = utime.ticks_us()\\n' + cell + '\\ndelta = utime.ticks_diff(utime.ticks_us(), t)' + \n", + " \"\\nprint('execution time: {:d} us'.format(delta))\")\n", + " print(message.decode('utf-8'))\n", + " \n", + " if args.memory: # prints out memory information \n", + " message = pyb.exec('from micropython import mem_info\\nprint(mem_info())\\n')\n", + " print(\"memory before execution:\\n========================\\n\", message.decode('utf-8'))\n", + " message = pyb.exec(cell)\n", + " print(\">>> \", message.decode('utf-8'))\n", + " message = pyb.exec('print(mem_info())')\n", + " print(\"memory after execution:\\n========================\\n\", message.decode('utf-8'))\n", + "\n", + " if args.pyboard:\n", + " message = pyb.exec(cell)\n", + " print(message.decode('utf-8'))\n", + "\n", + "ip = get_ipython()\n", + "ip.register_magics(PyboardMagic)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## pyboard" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": { + "ExecuteTime": { + "end_time": "2020-05-07T07:35:35.126401Z", + "start_time": "2020-05-07T07:35:35.105824Z" + } + }, + "outputs": [], + "source": [ + "import pyboard\n", + "pyb = pyboard.Pyboard('/dev/ttyACM0')\n", + "pyb.enter_raw_repl()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "ExecuteTime": { + "end_time": "2020-05-19T19:11:18.145548Z", + "start_time": "2020-05-19T19:11:18.137468Z" + } + }, + "outputs": [], + "source": [ + "pyb.exit_raw_repl()\n", + "pyb.close()" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "metadata": { + "ExecuteTime": { + "end_time": "2020-05-07T07:35:38.725924Z", + "start_time": "2020-05-07T07:35:38.645488Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], + "source": [ + "%%micropython -pyboard 1\n", + "\n", + "import utime\n", + "import ulab as np\n", + "\n", + "def timeit(n=1000):\n", + " def wrapper(f, *args, **kwargs):\n", + " func_name = str(f).split(' ')[1]\n", + " def new_func(*args, **kwargs):\n", + " run_times = np.zeros(n, dtype=np.uint16)\n", + " for i in range(n):\n", + " t = utime.ticks_us()\n", + " result = f(*args, **kwargs)\n", + " run_times[i] = utime.ticks_diff(utime.ticks_us(), t)\n", + " print('{}() execution times based on {} cycles'.format(func_name, n, (delta2-delta1)/n))\n", + " print('\\tbest: %d us'%np.min(run_times))\n", + " print('\\tworst: %d us'%np.max(run_times))\n", + " print('\\taverage: %d us'%np.mean(run_times))\n", + " print('\\tdeviation: +/-%.3f us'%np.std(run_times)) \n", + " return result\n", + " return new_func\n", + " return wrapper\n", + "\n", + "def timeit(f, *args, **kwargs):\n", + " func_name = str(f).split(' ')[1]\n", + " def new_func(*args, **kwargs):\n", + " t = utime.ticks_us()\n", + " result = f(*args, **kwargs)\n", + " print('execution time: ', utime.ticks_diff(utime.ticks_us(), t), ' us')\n", + " return result\n", + " return new_func" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__END_OF_DEFS__" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# scipy.integrate\n", + "\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.tanhsinh](#tanhsinh)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Introduction\n", + "\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", + "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, 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. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## quad\n", + "\n", + "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html \n", + "\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 arguments: \n", + "\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 etolerance) \n", + "\n", + "The function returns the result and the error estimate as a tuple of floats. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2020-06-19T20:24:10.529668Z", + "start_time": "2020-06-19T20:24:10.520389Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "UsageError: Cell magic `%%micropython` not found.\n" + ] + } + ], + "source": [ + "%%micropython -unix 1\n", + "\n", + "from ulab import scipy\n", + "\n", + "f = lambda x: x**2 + 2*x + 1\n", + "result = scipy.integrate.quad(f, 0, 5, order=5, eps=1e-10)\n", + "print (f\"result = {result}\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## romberg\n", + "\n", + "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romberg.html \n", + "\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", + "\n", + "The function takes three to five arguments: \n", + "\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 etolerance) \n", + "\n", + "The function returns the result as a float.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "UsageError: Cell magic `%%micropython` not found.\n" + ] + } + ], + "source": [ + "%%micropython -unix 1\n", + "\n", + "from ulab import scipy\n", + "\n", + "f = lambda x: x**2 + 2*x + 1\n", + "result = scipy.integrate.romberg(f, 0, 5)\n", + "print (f\"result = {result}\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## simpson\n", + "\n", + "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.simpson.html \n", + "\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 arguments: \n", + "\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 etolerance) \n", + "\n", + "The function returns the result as a float." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "UsageError: Cell magic `%%micropython` not found.\n" + ] + } + ], + "source": [ + "%%micropython -unix 1\n", + "\n", + "from ulab import scipy\n", + "\n", + "f = lambda x: x**2 + 2*x + 1\n", + "result = scipy.integrate.simpson(f, 0, 5)\n", + "print (f\"result = {result}\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 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", + "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", + "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 arguments: \n", + "\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: etolerance)\n", + "\n", + "The function returns the result and the error estimate as a tuple of floats.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "UsageError: Cell magic `%%micropython` not found.\n" + ] + } + ], + "source": [ + "%%micropython -unix 1\n", + "\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": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": { + "height": "calc(100% - 180px)", + "left": "10px", + "top": "150px", + "width": "382.797px" + }, + "toc_section_display": true, + "toc_window_display": true + }, + "varInspector": { + "cols": { + "lenName": 16, + "lenType": 16, + "lenVar": 40 + }, + "kernels_config": { + "python": { + "delete_cmd_postfix": "", + "delete_cmd_prefix": "del ", + "library": "var_list.py", + "varRefreshCmd": "print(var_dic_list())" + }, + "r": { + "delete_cmd_postfix": ") ", + "delete_cmd_prefix": "rm(", + "library": "var_list.r", + "varRefreshCmd": "cat(var_dic_list()) " + } + }, + "types_to_exclude": [ + "module", + "function", + "builtin_function_or_method", + "instance", + "_Feature" + ], + "window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/tests/2d/scipy/integrate.py b/tests/2d/scipy/integrate.py new file mode 100644 index 00000000..1d0edb7b --- /dev/null +++ b/tests/2d/scipy/integrate.py @@ -0,0 +1,28 @@ +import sys +from math import * + +try: + from ulab import scipy +except ImportError: + import scipy + +f = lambda x: x * sin(x) * exp(x) +a=1 +b=2 + +(res, err) = scipy.integrate.tanhsinh(f, a, b) +if isclose (res, 7.11263821415851) and isclose (err, 5.438231077315757e-14): + print (res, err) + +res = scipy.integrate.romberg(f, a, b) +if isclose (res, 7.112638214158507): + print (res) + +res = scipy.integrate.simpson(f, a, b) +if isclose (res, 7.112638214158494): + print (res) + +(res, err) = scipy.integrate.quad(f, a, b) +if isclose (res, 7.112638214158507) and isclose (err, 7.686723611780195e-14): + print (res, err) + diff --git a/tests/2d/scipy/integrate.py.exp b/tests/2d/scipy/integrate.py.exp new file mode 100644 index 00000000..10426e6c --- /dev/null +++ b/tests/2d/scipy/integrate.py.exp @@ -0,0 +1,4 @@ +7.11263821415851 5.438231077315757e-14 +7.112638214158507 +7.112638214158494 +7.112638214158507 7.686723611780195e-14