From 7c13c24b1c4fd9f5c2e3ce89b4dd76d1e10d1d18 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Thu, 5 Dec 2024 12:31:49 +0100 Subject: [PATCH 01/22] adding scipy integrate, initial check-in --- code/scipy/integrate/integrate.c | 649 +++++++++++++++++++++++++++++++ code/scipy/integrate/integrate.h | 52 +++ code/scipy/scipy.c | 6 +- code/ulab.h | 23 ++ 4 files changed, 729 insertions(+), 1 deletion(-) create mode 100644 code/scipy/integrate/integrate.c create mode 100644 code/scipy/integrate/integrate.h diff --git a/code/scipy/integrate/integrate.c b/code/scipy/integrate/integrate.c new file mode 100644 index 00000000..f9de7e90 --- /dev/null +++ b/code/scipy/integrate/integrate.c @@ -0,0 +1,649 @@ +/* + * This file is not 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" + +ULAB_DEFINE_FLOAT_CONST(etolerance, MICROPY_FLOAT_CONST(1e-14), 0x283424dcUL, 0x3e901b2b29a4692bULL); + +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 >= 0) + return 1; + else + return -1; +} + + +#if ULAB_INTEGRATE_HAS_QUAD +// 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 (MICROPY_FLOAT_C_FUN(isfinite)(h2) && MICROPY_FLOAT_C_FUN(fabs)(h2) > 1e-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 && !MICROPY_FLOAT_C_FUN(isfinite)(h)); + if (j > 1 && MICROPY_FLOAT_C_FUN(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 (MICROPY_FLOAT_C_FUN(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 != 0) // if last diff != 0, back up r by one step + r /= 2; + 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 quad(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 = 10*eps; + mp_float_t c = 0, d = 1, s, sign = 1, v, h = 2; + 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 (MICROPY_FLOAT_C_FUN(isfinite)(a) && MICROPY_FLOAT_C_FUN(isfinite)(b)) { + c = (a+b)/2; + d = (b-a)/2; + v = c; + } + else if (MICROPY_FLOAT_C_FUN(isfinite)(a)) { + mode = 1; // Exp-Sinh + d = exp_sinh_opt_d(fun, a, eps, d); + c = a; + v = a+d; + } + else if (MICROPY_FLOAT_C_FUN(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 = 0; + } + s = integrate_python_call(type, fun, v, fargs, 0); + do { + mp_float_t p = 0, q, fp = 0, fm = 0, t, eh; + h /= 2; + t = eh = MICROPY_FLOAT_C_FUN(exp)(h); + if (k > 0) + eh *= eh; + if (mode == 0) { // Tanh-Sinh + do { + mp_float_t u = MICROPY_FLOAT_C_FUN(exp)(1/t-t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2 + mp_float_t r = 2*u/(1+u); // = 1 - tanh(sinh(j*h)) + mp_float_t w = (t+1/t)*r/(1+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 (MICROPY_FLOAT_C_FUN(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 (MICROPY_FLOAT_C_FUN(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 /= 2; + do { + mp_float_t r = MICROPY_FLOAT_C_FUN(exp)(t-.25/t); // = exp(sinh(j*h)) + mp_float_t x, y, w = r; + q = 0; + 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 (MICROPY_FLOAT_C_FUN(isfinite)(y)) // if f(x) is finite, add to local sum + q += y/w; + } + else { // Sinh-Sinh + r = (r-1/r)/2; // = sinh(sinh(j*h)) + w = (w+1/w)/2; // = cosh(sinh(j*h)) + x = c - d*r; + y = integrate_python_call(type, fun, x, fargs, 0); + if (MICROPY_FLOAT_C_FUN(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 (MICROPY_FLOAT_C_FUN(isfinite)(y)) // if f(x) is finite, add to local sum + q += y*w; + q *= t+.25/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 quad( +//| fun: Callable[[float], float], +//| a: float, +//| b: float, +//| *, +//| levels: int = 6 +//| eps: float = 1e-14, +//| ) -> float: +//| """ +//| :param callable f: The function to integrate +//| :param float a: The left side of the interval +//| :param float b: The right side of the interval +//| :param float levels: The number of levels to perform (6..7 is optimal) +//| :param float eps: The error tolerance value +//| +//| Find a solution (zero) of the function ``f(x)`` on the interval +//| (``a``..``b``) using the bisection method. The result is accurate to within +//| ``xtol`` unless more than ``maxiter`` steps are 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_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 function")); + } + + 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 < 0) { + mp_raise_ValueError(MP_ERROR_TEXT("levels should be > 0")); + } + 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(quad(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_quad_obj, 2, integrate_quad); +#endif /* ULAB_INTEGRATE_HAS_QUAD */ + +#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 it 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 = 0; + mp_float_t *Rt; + h /= 2; + for (j = 1; j < k; j += 2) + sum += integrate_python_call(type, fun, a+j*h, fargs, 0); + Ru[0] = h*sum + Ro[0]/2; + 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 = 1e-14, +//| ) -> float: +//| """ +//| :param callable f: The function to integrate +//| :param float a: The left side of the interval +//| :param float b: The right side of the interval +//| :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 function")); + } + + 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 < 0) { + mp_raise_ValueError(MP_ERROR_TEXT("steps should be > 0")); + } + 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)/2; + mp_float_t f1 = integrate_python_call(type, fun, a + h/2, fargs, 0); + mp_float_t f2 = integrate_python_call(type, fun, b - h/2, fargs, 0); + mp_float_t sl = h*(fa + 4*f1 + fm)/6; + mp_float_t sr = h*(fm + 4*f2 + fb)/6; + mp_float_t s = sl+sr; + mp_float_t d = (s-v)/15; + 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 /= 2; + --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+4*fm+fb)*(b-a)/6; + 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 = 1e-14, +//| ) -> float: +//| """ +//| :param callable f: The function to integrate +//| :param float a: The left side of the interval +//| :param float b: The right side of the interval +//| :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")); + } + + 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 < 0) { + mp_raise_ValueError(MP_ERROR_TEXT("steps should be > 0")); + } + 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_QUADGK +// 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] = { + 0.00000000000000000e+00, + 7.65265211334973338e-02, + 1.52605465240922676e-01, + 2.27785851141645078e-01, + 3.01627868114913004e-01, + 3.73706088715419561e-01, + 4.43593175238725103e-01, + 5.10867001950827098e-01, + 5.75140446819710315e-01, + 6.36053680726515025e-01, + 6.93237656334751385e-01, + 7.46331906460150793e-01, + 7.95041428837551198e-01, + 8.39116971822218823e-01, + 8.78276811252281976e-01, + 9.12234428251325906e-01, + 9.40822633831754754e-01, + 9.63971927277913791e-01, + 9.81507877450250259e-01, + 9.93128599185094925e-01, + 9.98859031588277664e-01, + }; + static const mp_float_t weights[21] = { + 7.66007119179996564e-02, + 7.63778676720807367e-02, + 7.57044976845566747e-02, + 7.45828754004991890e-02, + 7.30306903327866675e-02, + 7.10544235534440683e-02, + 6.86486729285216193e-02, + 6.58345971336184221e-02, + 6.26532375547811680e-02, + 5.91114008806395724e-02, + 5.51951053482859947e-02, + 5.09445739237286919e-02, + 4.64348218674976747e-02, + 4.16688733279736863e-02, + 3.66001697582007980e-02, + 3.12873067770327990e-02, + 2.58821336049511588e-02, + 2.03883734612665236e-02, + 1.46261692569712530e-02, + 8.60026985564294220e-03, + 3.07358371852053150e-03, + }; + static const mp_float_t gauss_weights[10] = { + 1.52753387130725851e-01, + 1.49172986472603747e-01, + 1.42096109318382051e-01, + 1.31688638449176627e-01, + 1.18194531961518417e-01, + 1.01930119817240435e-01, + 8.32767415767047487e-02, + 6.26720483341090636e-02, + 4.06014298003869413e-02, + 1.76140071391521183e-02, + }; + const mp_obj_type_t *type = mp_obj_get_type(fun); + mp_obj_t fargs[1]; + mp_float_t p = 0; // kronrod quadrature sum + mp_float_t q = 0; // 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*1e-17); // 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)/2; + mp_float_t d = (b-a)/2; + 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 == 0) + tol = t; + if (n > 0 && t < e && tol < e) { + s = qakro(fun, a, c, n-1, t/2, eps, err); + s += qakro(fun, c, b, n-1, t/2, eps, &e); + *err += e; + return s; + } + *err = e; + return s; +} + + +//| def quadgk( +//| fun: Callable[[float], float], +//| a: float, +//| b: float, +//| *, +//| order: int = 5 +//| eps: float = 1e-14, +//| ) -> float: +//| """ +//| :param callable f: The function to integrate +//| :param float a: The left side of the interval +//| :param float b: The right side of the interval +//| :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_quadgk(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 function")); + } + + 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 < 0) { + mp_raise_ValueError(MP_ERROR_TEXT("levels should be > 0")); + } + 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_quadgk_obj, 2, integrate_quadgk); +#endif /* ULAB_INTEGRATE_HAS_QUADGK */ + +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_QUAD + { MP_ROM_QSTR(MP_QSTR_quad), MP_ROM_PTR(&integrate_quad_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_QUADGK + { MP_ROM_QSTR(MP_QSTR_quadgk), MP_ROM_PTR(&integrate_quadgk_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..d2ef4e3d --- /dev/null +++ b/code/scipy/integrate/integrate.h @@ -0,0 +1,52 @@ + +/* + * This file is not 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" + +/* +#ifndef OPTIMIZE_EPSILON +#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT +#define OPTIMIZE_EPSILON MICROPY_FLOAT_CONST(1.2e-7) +#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE +#define OPTIMIZE_EPSILON MICROPY_FLOAT_CONST(2.3e-16) +#endif +#endif + +#define OPTIMIZE_EPS MICROPY_FLOAT_CONST(1.0e-4) +#define OPTIMIZE_NONZDELTA MICROPY_FLOAT_CONST(0.05) +#define OPTIMIZE_ZDELTA MICROPY_FLOAT_CONST(0.00025) +#define OPTIMIZE_ALPHA MICROPY_FLOAT_CONST(1.0) +#define OPTIMIZE_BETA MICROPY_FLOAT_CONST(2.0) +#define OPTIMIZE_GAMMA MICROPY_FLOAT_CONST(0.5) +#define OPTIMIZE_DELTA MICROPY_FLOAT_CONST(0.5) +*/ + +extern const mp_obj_module_t ulab_scipy_integrate_module; + +#if ULAB_INTEGRATE_HAS_QUAD +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_quad_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_QUADGK +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_quadgk_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.h b/code/ulab.h index 78b8a1ba..79e015ba 100644 --- a/code/ulab.h +++ b/code/ulab.h @@ -398,6 +398,29 @@ #define ULAB_NUMPY_HAS_WHERE (1) #endif +// the integrate module +#ifndef ULAB_SCIPY_HAS_INTEGRATE_MODULE +#define ULAB_SCIPY_HAS_INTEGRATE_MODULE (1) +#endif + +#ifndef ULAB_INTEGRATE_HAS_QUAD +#define ULAB_INTEGRATE_HAS_QUAD (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_QUADGK +#define ULAB_INTEGRATE_HAS_QUADGK (1) +#endif + // the linalg module; functions of the linalg module still have // to be defined separately #ifndef ULAB_NUMPY_HAS_LINALG_MODULE From a4c0c0b90fd716c18f9a48b45c3bd662ca60067a Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Thu, 5 Dec 2024 13:21:20 +0100 Subject: [PATCH 02/22] compile unix double-precision, select integrations algos --- build.sh | 4 ++-- code/micropython.mk | 3 +++ code/scipy/integrate/integrate.c | 10 +++++++++- code/scipy/integrate/integrate.h | 2 +- 4 files changed, 15 insertions(+), 4 deletions(-) diff --git a/build.sh b/build.sh index 7927d4ac..47f777d7 100755 --- a/build.sh +++ b/build.sh @@ -59,8 +59,8 @@ fi bash test-common.sh "${dims}" "$PROG" -# Build with single-precision float. -make -C micropython/ports/unix -j${NPROC} USER_C_MODULES="${HERE}" DEBUG=1 STRIP=: MICROPY_PY_FFI=0 MICROPY_PY_BTREE=0 CFLAGS_EXTRA=-DMICROPY_FLOAT_IMPL=MICROPY_FLOAT_IMPL_FLOAT CFLAGS_EXTRA+=-DULAB_MAX_DIMS=$dims CFLAGS_EXTRA+=-DULAB_HASH=$GIT_HASH BUILD=build-nanbox-$dims PROG=micropython-nanbox-$dims +# Build with double-precision float. +make -C micropython/ports/unix -j${NPROC} USER_C_MODULES="${HERE}" DEBUG=1 STRIP=: MICROPY_PY_FFI=0 MICROPY_PY_BTREE=0 CFLAGS_EXTRA=-DMICROPY_FLOAT_IMPL=MICROPY_FLOAT_IMPL_DOUBLE CFLAGS_EXTRA+=-DULAB_MAX_DIMS=$dims CFLAGS_EXTRA+=-DULAB_HASH=$GIT_HASH BUILD=build-nanbox-$dims PROG=micropython-nanbox-$dims # The unix nanbox variant builds as a 32-bit executable and requires gcc-multilib. # macOS doesn't support i386 builds so only build on linux. diff --git a/code/micropython.mk b/code/micropython.mk index 4aa6f615..98c1c999 100644 --- a/code/micropython.mk +++ b/code/micropython.mk @@ -2,6 +2,9 @@ USERMODULES_DIR := $(USERMOD_DIR) # Add all C files to SRC_USERMOD. +#if ULAB_SCIPY_HAS_INTEGRATE_MODULE +SRC_USERMOD += $(USERMODULES_DIR)/scipy/integrate/integrate.c +#endif 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 index f9de7e90..71864692 100644 --- a/code/scipy/integrate/integrate.c +++ b/code/scipy/integrate/integrate.c @@ -1,5 +1,5 @@ /* - * This file is not part of the micropython-ulab project, + * This file is part of the micropython-ulab project, * * https://github.com/v923z/micropython-ulab * @@ -261,9 +261,11 @@ static mp_obj_t integrate_quad(size_t n_args, const mp_obj_t *pos_args, mp_map_t 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 0 if(n < 0) { mp_raise_ValueError(MP_ERROR_TEXT("levels should be > 0")); } +#endif mp_float_t eps = mp_obj_get_float(args[4].u_obj); mp_obj_t res[2]; @@ -355,9 +357,11 @@ static mp_obj_t integrate_romberg(size_t n_args, const mp_obj_t *pos_args, mp_ma 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 0 if(steps < 0) { mp_raise_ValueError(MP_ERROR_TEXT("steps should be > 0")); } +#endif mp_float_t eps = mp_obj_get_float(args[4].u_obj); return mp_obj_new_float(qromb(fun, a, b, steps, eps)); @@ -442,9 +446,11 @@ static mp_obj_t integrate_simpson(size_t n_args, const mp_obj_t *pos_args, mp_ma 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 0 if(steps < 0) { mp_raise_ValueError(MP_ERROR_TEXT("steps should be > 0")); } +#endif mp_float_t eps = mp_obj_get_float(args[4].u_obj); return mp_obj_new_float(qasi(fun, a, b, steps, eps)); @@ -606,9 +612,11 @@ static mp_obj_t integrate_quadgk(size_t n_args, const mp_obj_t *pos_args, mp_map 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 0 if(order < 0) { mp_raise_ValueError(MP_ERROR_TEXT("levels should be > 0")); } +#endif mp_float_t eps = mp_obj_get_float(args[4].u_obj); mp_obj_t res[2]; diff --git a/code/scipy/integrate/integrate.h b/code/scipy/integrate/integrate.h index d2ef4e3d..1dd42cc6 100644 --- a/code/scipy/integrate/integrate.h +++ b/code/scipy/integrate/integrate.h @@ -1,6 +1,6 @@ /* - * This file is not part of the micropython-ulab project, + * This file is part of the micropython-ulab project, * * https://github.com/v923z/micropython-ulab * From 4a17c92aef668d2db76ac6dc025833dadc57831f Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Thu, 5 Dec 2024 14:04:40 +0100 Subject: [PATCH 03/22] bumping ulab version number to 6.7.0 --- code/ulab.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From e50df010d07d70cbfd004989c3f243a4263be264 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Thu, 5 Dec 2024 14:58:59 +0100 Subject: [PATCH 04/22] adding documentation --- docs/scipy-integrate.ipynb | 481 +++++++++++++++++++++++++++++++++++++ 1 file changed, 481 insertions(+) create mode 100644 docs/scipy-integrate.ipynb diff --git a/docs/scipy-integrate.ipynb b/docs/scipy-integrate.ipynb new file mode 100644 index 00000000..903b3bec --- /dev/null +++ b/docs/scipy-integrate.ipynb @@ -0,0 +1,481 @@ +{ + "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 defines four different numeric 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" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## quad\n", + "\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, ported from code found in https://www.genivia.com/qthsh.html. 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). Please check the paper above for a comprehensive discussion.\n", + "\n", + "The function takes three to five positional arguments: \n", + "\n", + "* f, a function defined e.g. with lambda,\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", + "\n", + "The function returns the result and the final error tolerance as a tuple of mp_float_t. " + ] + }, + { + "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 as sc\n", + "\n", + "f = lambda x: x**2 + 2*x + 1\n", + "result = sc.integrate.quad(f, 0, 5)\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", + "Implements the Romberg quadrature algorithm, ported from code found in https://www.genivia.com/qthsh.html. 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). Please check the paper above for a comprehensive discussion.\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 that this function does not support complex numbers. \n", + "\n", + "The function takes three to five positional arguments: \n", + "\n", + "* f, a function defined e.g. with lambda,\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", + "\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 as sc\n", + "\n", + "f = lambda x: x**2 + 2*x + 1\n", + "result = sc.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", + "Implements an Adaptive Simpson quadrature algorithm, ported from code found in https://www.genivia.com/qthsh.html. 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). Please check the paper above for a comprehensive discussion.\n", + "\n", + "Please note that this function does not support complex numbers. \n", + "\n", + "The function takes three to five positional arguments: \n", + "\n", + "* f, a function defined e.g. with lambda,\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", + "\n", + "The function returns the result as a mp_float_t." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "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.simpson(f, 0, 5)\n", + "print (f\"result = {result}\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## quadgk\n", + "\n", + "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quadgk.html \n", + "\n", + "Implements an Adaptive Gauss-Kronrod (G10,K21) quadrature algorithm, ported from code found in https://www.genivia.com/qthsh.html. The Gauss–Kronrod quadrature formula is an adaptive method for numerical integration. It 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). Please check the paper above for a comprehensive discussion.\n", + "\n", + "The function takes three to five positional arguments: \n", + "\n", + "* f, a function defined e.g. with lambda,\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", + "\n", + "The function returns the result and the final error tolerance as a tuple of mp_float_t. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "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" + ] + }, + { + "cell_type": "markdown", + "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 will fail. " + ] + } + ], + "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 +} From a3bd00c6d9f24ab9ae0a1f6cd2bcbff717540e1b Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Thu, 5 Dec 2024 15:02:08 +0100 Subject: [PATCH 05/22] documentation fix --- docs/scipy-integrate.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/scipy-integrate.ipynb b/docs/scipy-integrate.ipynb index 903b3bec..19b4cf99 100644 --- a/docs/scipy-integrate.ipynb +++ b/docs/scipy-integrate.ipynb @@ -278,7 +278,7 @@ "from ulab import scipy as sc\n", "\n", "f = lambda x: x**2 + 2*x + 1\n", - "result = sc.integrate.quad(f, 0, 5)\n", + "result = sc.integrate.quad(f, 0, 5, levels=6, eps=1e-14)\n", "print (f\"result = {result}\")\n" ] }, From ba3ac29069eacd1a88fb83e5103baf79cbc1d41e Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Thu, 5 Dec 2024 15:03:15 +0100 Subject: [PATCH 06/22] documentation fix --- docs/scipy-integrate.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/scipy-integrate.ipynb b/docs/scipy-integrate.ipynb index 19b4cf99..49052c17 100644 --- a/docs/scipy-integrate.ipynb +++ b/docs/scipy-integrate.ipynb @@ -406,7 +406,7 @@ "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 will fail. " + "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 will fail. " ] } ], From 6e878c5f3c779c1df30dba3998de4b9aa70f9b77 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Thu, 5 Dec 2024 15:25:29 +0100 Subject: [PATCH 07/22] documentation fix --- docs/scipy-integrate.ipynb | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/scipy-integrate.ipynb b/docs/scipy-integrate.ipynb index 49052c17..aa755efa 100644 --- a/docs/scipy-integrate.ipynb +++ b/docs/scipy-integrate.ipynb @@ -406,7 +406,9 @@ "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 will fail. " + "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 will 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`. \n" ] } ], From 0a6ad6d5aab1981f40c918886cbca51d8afff5a6 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Sun, 8 Dec 2024 13:08:34 +0100 Subject: [PATCH 08/22] rewritten in some places --- docs/scipy-integrate.ipynb | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/docs/scipy-integrate.ipynb b/docs/scipy-integrate.ipynb index aa755efa..ca0c50af 100644 --- a/docs/scipy-integrate.ipynb +++ b/docs/scipy-integrate.ipynb @@ -234,6 +234,19 @@ "\n" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Motivation\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", + "\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", + "\n", + "Also note that these algorithms do not support complex numbers. " + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -242,7 +255,7 @@ "\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, ported from code found in https://www.genivia.com/qthsh.html. 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). Please check the paper above for a comprehensive discussion.\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", "\n", "The function takes three to five positional arguments: \n", "\n", @@ -290,10 +303,9 @@ "\n", "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romberg.html \n", "\n", - "Implements the Romberg quadrature algorithm, ported from code found in https://www.genivia.com/qthsh.html. 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). Please check the paper above for a comprehensive discussion.\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", "\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 that this function does not support complex numbers. \n", "\n", "The function takes three to five positional arguments: \n", "\n", @@ -336,9 +348,7 @@ "\n", "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.simpson.html \n", "\n", - "Implements an Adaptive Simpson quadrature algorithm, ported from code found in https://www.genivia.com/qthsh.html. 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). Please check the paper above for a comprehensive discussion.\n", - "\n", - "Please note that this function does not support complex numbers. \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", "\n", "The function takes three to five positional arguments: \n", "\n", @@ -373,7 +383,7 @@ "\n", "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quadgk.html \n", "\n", - "Implements an Adaptive Gauss-Kronrod (G10,K21) quadrature algorithm, ported from code found in https://www.genivia.com/qthsh.html. The Gauss–Kronrod quadrature formula is an adaptive method for numerical integration. It 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). Please check the paper above for a comprehensive discussion.\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", "\n", "The function takes three to five positional arguments: \n", "\n", @@ -406,9 +416,9 @@ "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 will fail. \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`. \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" ] } ], From 59c0e716822905c710d0a260cfc4c53022884e72 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Sun, 8 Dec 2024 13:09:02 +0100 Subject: [PATCH 09/22] complex number error handling --- code/scipy/integrate/integrate.c | 68 ++++++++++++++++++++++---------- code/ulab.h | 7 +++- 2 files changed, 53 insertions(+), 22 deletions(-) diff --git a/code/scipy/integrate/integrate.c b/code/scipy/integrate/integrate.c index 71864692..db80a7ac 100644 --- a/code/scipy/integrate/integrate.c +++ b/code/scipy/integrate/integrate.c @@ -258,14 +258,21 @@ static mp_obj_t integrate_quad(size_t n_args, const mp_obj_t *pos_args, mp_map_t 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 n = (uint16_t)args[3].u_int; -#if 0 - if(n < 0) { - mp_raise_ValueError(MP_ERROR_TEXT("levels should be > 0")); - } -#endif + 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]; @@ -354,14 +361,21 @@ static mp_obj_t integrate_romberg(size_t n_args, const mp_obj_t *pos_args, mp_ma 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 0 - if(steps < 0) { - mp_raise_ValueError(MP_ERROR_TEXT("steps should be > 0")); - } -#endif + 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)); @@ -443,14 +457,21 @@ static mp_obj_t integrate_simpson(size_t n_args, const mp_obj_t *pos_args, mp_ma 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 0 - if(steps < 0) { - mp_raise_ValueError(MP_ERROR_TEXT("steps should be > 0")); - } -#endif + 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)); @@ -609,14 +630,21 @@ static mp_obj_t integrate_quadgk(size_t n_args, const mp_obj_t *pos_args, mp_map 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 order = (uint16_t)args[3].u_int; -#if 0 - if(order < 0) { - mp_raise_ValueError(MP_ERROR_TEXT("levels should be > 0")); - } -#endif + 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]; diff --git a/code/ulab.h b/code/ulab.h index 79e015ba..7b5e7f40 100644 --- a/code/ulab.h +++ b/code/ulab.h @@ -398,7 +398,8 @@ #define ULAB_NUMPY_HAS_WHERE (1) #endif -// the integrate module +// 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 @@ -416,10 +417,12 @@ #define ULAB_INTEGRATE_HAS_SIMPSON (1) #endif - +#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE +// we compile quadgk only when we have _DOUBLE #ifndef ULAB_INTEGRATE_HAS_QUADGK #define ULAB_INTEGRATE_HAS_QUADGK (1) #endif +#endif // the linalg module; functions of the linalg module still have // to be defined separately From f446db0ea7925b82d4948b089ae52f2eb3392c4e Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Sun, 8 Dec 2024 13:09:22 +0100 Subject: [PATCH 10/22] added test cases --- tests/2d/scipy/integrate.py | 29 +++++++++++++++++++++++++++++ tests/2d/scipy/integrate.py.exp | 4 ++++ 2 files changed, 33 insertions(+) create mode 100644 tests/2d/scipy/integrate.py create mode 100644 tests/2d/scipy/integrate.py.exp diff --git a/tests/2d/scipy/integrate.py b/tests/2d/scipy/integrate.py new file mode 100644 index 00000000..39b163bd --- /dev/null +++ b/tests/2d/scipy/integrate.py @@ -0,0 +1,29 @@ +import sys +from math import * + +# this test is meaningful only if ulab is compiled with ulab.scipy.integrate +# we're not going to test CPython scipy. +try: + from ulab import scipy +except Exception as e: + print ("could not import ulab.scipy: ", e) + sys.exit(1) + +i = scipy.integrate +try: + if str(type(i)) != "": + print ("scipy.integrate is not available") + sys.exit(1) +except Exception as e: + print ("scipy.integrate is not available") + sys.exit(1) + +f = lambda x: x * sin(x) * exp(x) +a=1 +b=2 +# what if they are not all available? +algorithms = ( i.quad, i.romberg, i.simpson, i.quadgk ) +for quad in algorithms: + print (quad(f, a, b)) + + diff --git a/tests/2d/scipy/integrate.py.exp b/tests/2d/scipy/integrate.py.exp new file mode 100644 index 00000000..6a18b121 --- /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) From 6b3930230d8544d3f91e1ec84956656ff8411c03 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Fri, 13 Dec 2024 10:46:07 +0100 Subject: [PATCH 11/22] resolved importing scipy.integrate --- tests/2d/scipy/integrate.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/tests/2d/scipy/integrate.py b/tests/2d/scipy/integrate.py index 39b163bd..dbd8f1be 100644 --- a/tests/2d/scipy/integrate.py +++ b/tests/2d/scipy/integrate.py @@ -2,16 +2,14 @@ from math import * # this test is meaningful only if ulab is compiled with ulab.scipy.integrate -# we're not going to test CPython scipy. try: - from ulab import scipy + from ulab import scipy except Exception as e: print ("could not import ulab.scipy: ", e) sys.exit(1) -i = scipy.integrate try: - if str(type(i)) != "": + if str(type(scipy.integrate)) != "": print ("scipy.integrate is not available") sys.exit(1) except Exception as e: @@ -21,9 +19,8 @@ f = lambda x: x * sin(x) * exp(x) a=1 b=2 -# what if they are not all available? -algorithms = ( i.quad, i.romberg, i.simpson, i.quadgk ) -for quad in algorithms: - print (quad(f, a, b)) +functions = ( scipy.integrate.quad, scipy.integrate.romberg, scipy.integrate.simpson, scipy.integrate.quadgk ) +for function in functions: + print (function(f, a, b)) From 40397fbc5a0427e48a795f5c23d12c24a0e37664 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Fri, 13 Dec 2024 10:49:25 +0100 Subject: [PATCH 12/22] resolved importing scipy.integrate #2 --- tests/2d/scipy/integrate.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/tests/2d/scipy/integrate.py b/tests/2d/scipy/integrate.py index dbd8f1be..24ff382b 100644 --- a/tests/2d/scipy/integrate.py +++ b/tests/2d/scipy/integrate.py @@ -1,13 +1,8 @@ import sys from math import * +from ulab import scipy # this test is meaningful only if ulab is compiled with ulab.scipy.integrate -try: - from ulab import scipy -except Exception as e: - print ("could not import ulab.scipy: ", e) - sys.exit(1) - try: if str(type(scipy.integrate)) != "": print ("scipy.integrate is not available") From f4fc606ddbce9469ced397d73e8a0e7e0ee658c3 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Fri, 13 Dec 2024 10:53:43 +0100 Subject: [PATCH 13/22] build integrate only when we have MICROPY_FLOAT_IMPL_DOUBLE --- code/ulab.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/code/ulab.h b/code/ulab.h index 7b5e7f40..a859b1f1 100644 --- a/code/ulab.h +++ b/code/ulab.h @@ -400,6 +400,8 @@ // the integrate module; functions of the integrate module still have // to be defined separately +#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE +// we compile integrate only when we have _DOUBLE #ifndef ULAB_SCIPY_HAS_INTEGRATE_MODULE #define ULAB_SCIPY_HAS_INTEGRATE_MODULE (1) #endif @@ -417,8 +419,6 @@ #define ULAB_INTEGRATE_HAS_SIMPSON (1) #endif -#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE -// we compile quadgk only when we have _DOUBLE #ifndef ULAB_INTEGRATE_HAS_QUADGK #define ULAB_INTEGRATE_HAS_QUADGK (1) #endif From 87e669dfad09aac0e0e45e8d9014c8f2c2c113d9 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Fri, 13 Dec 2024 10:56:12 +0100 Subject: [PATCH 14/22] reverting commit a4c0c0b --- build.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/build.sh b/build.sh index 47f777d7..7927d4ac 100755 --- a/build.sh +++ b/build.sh @@ -59,8 +59,8 @@ fi bash test-common.sh "${dims}" "$PROG" -# Build with double-precision float. -make -C micropython/ports/unix -j${NPROC} USER_C_MODULES="${HERE}" DEBUG=1 STRIP=: MICROPY_PY_FFI=0 MICROPY_PY_BTREE=0 CFLAGS_EXTRA=-DMICROPY_FLOAT_IMPL=MICROPY_FLOAT_IMPL_DOUBLE CFLAGS_EXTRA+=-DULAB_MAX_DIMS=$dims CFLAGS_EXTRA+=-DULAB_HASH=$GIT_HASH BUILD=build-nanbox-$dims PROG=micropython-nanbox-$dims +# Build with single-precision float. +make -C micropython/ports/unix -j${NPROC} USER_C_MODULES="${HERE}" DEBUG=1 STRIP=: MICROPY_PY_FFI=0 MICROPY_PY_BTREE=0 CFLAGS_EXTRA=-DMICROPY_FLOAT_IMPL=MICROPY_FLOAT_IMPL_FLOAT CFLAGS_EXTRA+=-DULAB_MAX_DIMS=$dims CFLAGS_EXTRA+=-DULAB_HASH=$GIT_HASH BUILD=build-nanbox-$dims PROG=micropython-nanbox-$dims # The unix nanbox variant builds as a 32-bit executable and requires gcc-multilib. # macOS doesn't support i386 builds so only build on linux. From a10e89fe14ce17306bb35e629509c01b1b6f8094 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Fri, 13 Dec 2024 11:10:39 +0100 Subject: [PATCH 15/22] re-pushing failed commit --- tests/2d/scipy/integrate.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/2d/scipy/integrate.py b/tests/2d/scipy/integrate.py index 24ff382b..52d5cc21 100644 --- a/tests/2d/scipy/integrate.py +++ b/tests/2d/scipy/integrate.py @@ -18,4 +18,3 @@ for function in functions: print (function(f, a, b)) - From af290c155b8e6fb3dbc7c1843228a92f23e8e1d1 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Fri, 13 Dec 2024 11:14:27 +0100 Subject: [PATCH 16/22] Revert "re-pushing failed commit" This reverts commit a10e89fe14ce17306bb35e629509c01b1b6f8094. --- tests/2d/scipy/integrate.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/2d/scipy/integrate.py b/tests/2d/scipy/integrate.py index 52d5cc21..24ff382b 100644 --- a/tests/2d/scipy/integrate.py +++ b/tests/2d/scipy/integrate.py @@ -18,3 +18,4 @@ for function in functions: print (function(f, a, b)) + From 7e5b67d939ea6287702fa268bf9ef78ebf30991d Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Sat, 14 Dec 2024 15:00:55 +0100 Subject: [PATCH 17/22] improve tests using math.isclose() --- tests/2d/scipy/integrate.py | 18 +++++++++++++++--- tests/2d/scipy/integrate.py.exp | 4 ++-- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/tests/2d/scipy/integrate.py b/tests/2d/scipy/integrate.py index 24ff382b..a17110f4 100644 --- a/tests/2d/scipy/integrate.py +++ b/tests/2d/scipy/integrate.py @@ -14,8 +14,20 @@ f = lambda x: x * sin(x) * exp(x) a=1 b=2 -functions = ( scipy.integrate.quad, scipy.integrate.romberg, scipy.integrate.simpson, scipy.integrate.quadgk ) -for function in functions: - print (function(f, a, b)) +(res, err) = scipy.integrate.quad(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.quadgk(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 index 6a18b121..10426e6c 100644 --- a/tests/2d/scipy/integrate.py.exp +++ b/tests/2d/scipy/integrate.py.exp @@ -1,4 +1,4 @@ -(7.11263821415851, 5.438231077315757e-14) +7.11263821415851 5.438231077315757e-14 7.112638214158507 7.112638214158494 -(7.112638214158507, 7.686723611780195e-14) +7.112638214158507 7.686723611780195e-14 From cb7166385451e2128702a8526087193e5139c69d Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Sat, 14 Dec 2024 16:02:09 +0100 Subject: [PATCH 18/22] enabled fp32 builds --- code/scipy/integrate/integrate.c | 223 +++++++++++++++++-------------- code/scipy/integrate/integrate.h | 18 --- code/ulab.h | 3 - 3 files changed, 120 insertions(+), 124 deletions(-) diff --git a/code/scipy/integrate/integrate.c b/code/scipy/integrate/integrate.c index db80a7ac..9a48d746 100644 --- a/code/scipy/integrate/integrate.c +++ b/code/scipy/integrate/integrate.c @@ -39,7 +39,24 @@ #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, ...), @@ -51,7 +68,7 @@ static mp_float_t integrate_python_call(const mp_obj_type_t *type, mp_obj_t fun, // sign helper function int sign(mp_float_t x) { - if (x >= 0) + if (x >= ZERO) return 1; else return -1; @@ -68,7 +85,7 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_ 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 (MICROPY_FLOAT_C_FUN(isfinite)(h2) && MICROPY_FLOAT_C_FUN(fabs)(h2) > 1e-5) { // if |h2| > 2^-16 + 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; @@ -76,8 +93,8 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_ 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 && !MICROPY_FLOAT_C_FUN(isfinite)(h)); - if (j > 1 && MICROPY_FLOAT_C_FUN(isfinite)(h) && sign(h) != sign(h2)) { + } 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 @@ -86,7 +103,7 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_ 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 (MICROPY_FLOAT_C_FUN(isfinite)(h)) { + 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 @@ -101,8 +118,8 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_ 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 != 0) // if last diff != 0, back up r by one step - r /= 2; + 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 @@ -118,8 +135,8 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_ mp_float_t quad(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 = 10*eps; - mp_float_t c = 0, d = 1, s, sign = 1, v, h = 2; + 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; @@ -127,18 +144,18 @@ mp_float_t quad(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint1 a = v; sign = -1; } - if (MICROPY_FLOAT_C_FUN(isfinite)(a) && MICROPY_FLOAT_C_FUN(isfinite)(b)) { - c = (a+b)/2; - d = (b-a)/2; + if (isfinite(a) && isfinite(b)) { + c = (a+b)/TWO; + d = (b-a)/TWO; v = c; } - else if (MICROPY_FLOAT_C_FUN(isfinite)(a)) { + else if (isfinite(a)) { mode = 1; // Exp-Sinh d = exp_sinh_opt_d(fun, a, eps, d); c = a; v = a+d; } - else if (MICROPY_FLOAT_C_FUN(isfinite)(b)) { + else if (isfinite(b)) { mode = 1; // Exp-Sinh // d = -d; d = exp_sinh_opt_d(fun, b, eps, -d); @@ -148,29 +165,29 @@ mp_float_t quad(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint1 } else { mode = 2; // Sinh-Sinh - v = 0; + v = ZERO; } s = integrate_python_call(type, fun, v, fargs, 0); do { - mp_float_t p = 0, q, fp = 0, fm = 0, t, eh; - h /= 2; + mp_float_t p = ZERO, q, fp = ZERO, fm = ZERO, t, eh; + h /= TWO; t = eh = MICROPY_FLOAT_C_FUN(exp)(h); - if (k > 0) + if (k > ZERO) eh *= eh; if (mode == 0) { // Tanh-Sinh do { - mp_float_t u = MICROPY_FLOAT_C_FUN(exp)(1/t-t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2 - mp_float_t r = 2*u/(1+u); // = 1 - tanh(sinh(j*h)) - mp_float_t w = (t+1/t)*r/(1+u); // = cosh(j*h)/cosh(sinh(j*h))^2 + 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 (MICROPY_FLOAT_C_FUN(isfinite)(y)) + 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 (MICROPY_FLOAT_C_FUN(isfinite)(y)) + if (isfinite(y)) fm = y; // if f(x) is finite, add to local sum } q = w*(fp+fm); @@ -179,32 +196,32 @@ mp_float_t quad(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint1 } while (MICROPY_FLOAT_C_FUN(fabs)(q) > eps*MICROPY_FLOAT_C_FUN(fabs)(p)); } else { - t /= 2; + t /= TWO; do { - mp_float_t r = MICROPY_FLOAT_C_FUN(exp)(t-.25/t); // = exp(sinh(j*h)) + 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 = 0; + 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 (MICROPY_FLOAT_C_FUN(isfinite)(y)) // if f(x) is finite, add to local sum + if (isfinite(y)) // if f(x) is finite, add to local sum q += y/w; } else { // Sinh-Sinh - r = (r-1/r)/2; // = sinh(sinh(j*h)) - w = (w+1/w)/2; // = cosh(sinh(j*h)) + 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 (MICROPY_FLOAT_C_FUN(isfinite)(y)) // if f(x) is finite, add to local sum + 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 (MICROPY_FLOAT_C_FUN(isfinite)(y)) // if f(x) is finite, add to local sum + if (isfinite(y)) // if f(x) is finite, add to local sum q += y*w; - q *= t+.25/t; // q *= cosh(j*h) + 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)); @@ -303,12 +320,12 @@ mp_float_t qromb(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint for (i = 1; i < n; ++i) { unsigned long long k = 1UL << i; unsigned long long s = 1; - mp_float_t sum = 0; + mp_float_t sum = ZERO; mp_float_t *Rt; - h /= 2; + 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]/2; + 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); @@ -392,17 +409,17 @@ mp_float_t as(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, mp_floa 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)/2; - mp_float_t f1 = integrate_python_call(type, fun, a + h/2, fargs, 0); - mp_float_t f2 = integrate_python_call(type, fun, b - h/2, fargs, 0); - mp_float_t sl = h*(fa + 4*f1 + fm)/6; - mp_float_t sr = h*(fm + 4*f2 + fb)/6; + 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)/15; + 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 /= 2; + 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); @@ -414,7 +431,7 @@ mp_float_t qasi(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n 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+4*fm+fb)*(b-a)/6; + mp_float_t v = (fa + FOUR * fm + fb) * (b-a) / SIX; return as(fun, a, b, fa, fm, fb, v, eps, n, 0); } @@ -487,67 +504,67 @@ MP_DEFINE_CONST_FUN_OBJ_KW(integrate_simpson_obj, 2, integrate_simpson); 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] = { - 0.00000000000000000e+00, - 7.65265211334973338e-02, - 1.52605465240922676e-01, - 2.27785851141645078e-01, - 3.01627868114913004e-01, - 3.73706088715419561e-01, - 4.43593175238725103e-01, - 5.10867001950827098e-01, - 5.75140446819710315e-01, - 6.36053680726515025e-01, - 6.93237656334751385e-01, - 7.46331906460150793e-01, - 7.95041428837551198e-01, - 8.39116971822218823e-01, - 8.78276811252281976e-01, - 9.12234428251325906e-01, - 9.40822633831754754e-01, - 9.63971927277913791e-01, - 9.81507877450250259e-01, - 9.93128599185094925e-01, - 9.98859031588277664e-01, + 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] = { - 7.66007119179996564e-02, - 7.63778676720807367e-02, - 7.57044976845566747e-02, - 7.45828754004991890e-02, - 7.30306903327866675e-02, - 7.10544235534440683e-02, - 6.86486729285216193e-02, - 6.58345971336184221e-02, - 6.26532375547811680e-02, - 5.91114008806395724e-02, - 5.51951053482859947e-02, - 5.09445739237286919e-02, - 4.64348218674976747e-02, - 4.16688733279736863e-02, - 3.66001697582007980e-02, - 3.12873067770327990e-02, - 2.58821336049511588e-02, - 2.03883734612665236e-02, - 1.46261692569712530e-02, - 8.60026985564294220e-03, - 3.07358371852053150e-03, + 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] = { - 1.52753387130725851e-01, - 1.49172986472603747e-01, - 1.42096109318382051e-01, - 1.31688638449176627e-01, - 1.18194531961518417e-01, - 1.01930119817240435e-01, - 8.32767415767047487e-02, - 6.26720483341090636e-02, - 4.06014298003869413e-02, - 1.76140071391521183e-02, + 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 = 0; // kronrod quadrature sum - mp_float_t q = 0; // gauss quadrature sum + 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; @@ -565,24 +582,24 @@ mp_float_t gk(mp_float_t (*fun)(mp_float_t), mp_float_t c, mp_float_t d, mp_floa p += (fp + fm) * weights[i]; } *err = MICROPY_FLOAT_C_FUN(fabs)(p - q); - e = MICROPY_FLOAT_C_FUN(fabs)(2*p*1e-17); // optional, to take 1e-17 MachEps prec. into account + 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)/2; - mp_float_t d = (b-a)/2; + 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 == 0) + if (tol == ZERO) tol = t; if (n > 0 && t < e && tol < e) { - s = qakro(fun, a, c, n-1, t/2, eps, err); - s += qakro(fun, c, b, n-1, t/2, eps, &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; } diff --git a/code/scipy/integrate/integrate.h b/code/scipy/integrate/integrate.h index 1dd42cc6..0fea59e8 100644 --- a/code/scipy/integrate/integrate.h +++ b/code/scipy/integrate/integrate.h @@ -15,24 +15,6 @@ #include "../../ulab_tools.h" -/* -#ifndef OPTIMIZE_EPSILON -#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT -#define OPTIMIZE_EPSILON MICROPY_FLOAT_CONST(1.2e-7) -#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE -#define OPTIMIZE_EPSILON MICROPY_FLOAT_CONST(2.3e-16) -#endif -#endif - -#define OPTIMIZE_EPS MICROPY_FLOAT_CONST(1.0e-4) -#define OPTIMIZE_NONZDELTA MICROPY_FLOAT_CONST(0.05) -#define OPTIMIZE_ZDELTA MICROPY_FLOAT_CONST(0.00025) -#define OPTIMIZE_ALPHA MICROPY_FLOAT_CONST(1.0) -#define OPTIMIZE_BETA MICROPY_FLOAT_CONST(2.0) -#define OPTIMIZE_GAMMA MICROPY_FLOAT_CONST(0.5) -#define OPTIMIZE_DELTA MICROPY_FLOAT_CONST(0.5) -*/ - extern const mp_obj_module_t ulab_scipy_integrate_module; #if ULAB_INTEGRATE_HAS_QUAD diff --git a/code/ulab.h b/code/ulab.h index a859b1f1..22f32337 100644 --- a/code/ulab.h +++ b/code/ulab.h @@ -400,8 +400,6 @@ // the integrate module; functions of the integrate module still have // to be defined separately -#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE -// we compile integrate only when we have _DOUBLE #ifndef ULAB_SCIPY_HAS_INTEGRATE_MODULE #define ULAB_SCIPY_HAS_INTEGRATE_MODULE (1) #endif @@ -422,7 +420,6 @@ #ifndef ULAB_INTEGRATE_HAS_QUADGK #define ULAB_INTEGRATE_HAS_QUADGK (1) #endif -#endif // the linalg module; functions of the linalg module still have // to be defined separately From 91f96b69ceeac8460934551798a519360b873d92 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Sun, 15 Dec 2024 11:28:26 +0100 Subject: [PATCH 19/22] removed conditional includes --- code/micropython.mk | 2 -- 1 file changed, 2 deletions(-) diff --git a/code/micropython.mk b/code/micropython.mk index 98c1c999..e835d87b 100644 --- a/code/micropython.mk +++ b/code/micropython.mk @@ -2,9 +2,7 @@ USERMODULES_DIR := $(USERMOD_DIR) # Add all C files to SRC_USERMOD. -#if ULAB_SCIPY_HAS_INTEGRATE_MODULE SRC_USERMOD += $(USERMODULES_DIR)/scipy/integrate/integrate.c -#endif SRC_USERMOD += $(USERMODULES_DIR)/scipy/linalg/linalg.c SRC_USERMOD += $(USERMODULES_DIR)/scipy/optimize/optimize.c SRC_USERMOD += $(USERMODULES_DIR)/scipy/signal/signal.c From 97b27d3173ae347cea630202620b744ce6cd93e1 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Sun, 15 Dec 2024 11:48:52 +0100 Subject: [PATCH 20/22] adapted to new function names, corrected importing --- tests/2d/scipy/integrate.py | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/tests/2d/scipy/integrate.py b/tests/2d/scipy/integrate.py index a17110f4..1d0edb7b 100644 --- a/tests/2d/scipy/integrate.py +++ b/tests/2d/scipy/integrate.py @@ -1,21 +1,16 @@ import sys from math import * -from ulab import scipy -# this test is meaningful only if ulab is compiled with ulab.scipy.integrate try: - if str(type(scipy.integrate)) != "": - print ("scipy.integrate is not available") - sys.exit(1) -except Exception as e: - print ("scipy.integrate is not available") - sys.exit(1) - + from ulab import scipy +except ImportError: + import scipy + f = lambda x: x * sin(x) * exp(x) a=1 b=2 -(res, err) = scipy.integrate.quad(f, a, b) +(res, err) = scipy.integrate.tanhsinh(f, a, b) if isclose (res, 7.11263821415851) and isclose (err, 5.438231077315757e-14): print (res, err) @@ -27,7 +22,7 @@ if isclose (res, 7.112638214158494): print (res) -(res, err) = scipy.integrate.quadgk(f, a, b) +(res, err) = scipy.integrate.quad(f, a, b) if isclose (res, 7.112638214158507) and isclose (err, 7.686723611780195e-14): print (res, err) From 2456e16d790f0d39a6df58cacb700d9f59d20260 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Sun, 15 Dec 2024 11:50:13 +0100 Subject: [PATCH 21/22] function names similar to in CPython scipy.integrate, some minor corrections --- code/scipy/integrate/integrate.c | 73 ++++++++++++++++---------------- code/scipy/integrate/integrate.h | 8 ++-- code/ulab.h | 9 ++-- 3 files changed, 44 insertions(+), 46 deletions(-) diff --git a/code/scipy/integrate/integrate.c b/code/scipy/integrate/integrate.c index 9a48d746..80def711 100644 --- a/code/scipy/integrate/integrate.c +++ b/code/scipy/integrate/integrate.c @@ -75,7 +75,7 @@ int sign(mp_float_t x) { } -#if ULAB_INTEGRATE_HAS_QUAD +#if ULAB_INTEGRATE_HAS_TANHSINH // Tanh-Sinh, Sinh-Sinh and Exp-Sinh quadrature // https://www.genivia.com/qthsh.html @@ -132,7 +132,7 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_ // integrate function f, range a..b, max levels n, error tolerance eps -mp_float_t quad(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) { +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; @@ -235,29 +235,28 @@ mp_float_t quad(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint1 return sign*d*s*h; // result with estimated relative error e } -//| def quad( +//| def tanhsinh( //| fun: Callable[[float], float], //| a: float, //| b: float, //| *, //| levels: int = 6 -//| eps: float = 1e-14, +//| eps: float = etolerance //| ) -> float: //| """ //| :param callable f: The function to integrate -//| :param float a: The left side of the interval -//| :param float b: The right side of the interval +//| :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 +//| :param float eps: The error tolerance value //| -//| Find a solution (zero) of the function ``f(x)`` on the interval -//| (``a``..``b``) using the bisection method. The result is accurate to within -//| ``xtol`` unless more than ``maxiter`` steps are required.""" -//| ... +//| 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_quad(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { +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 } }, @@ -272,7 +271,7 @@ static mp_obj_t integrate_quad(size_t n_args, const mp_obj_t *pos_args, mp_map_t 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")); + mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a callable")); } // iterate over args 1, 2, and 4 @@ -294,20 +293,20 @@ static mp_obj_t integrate_quad(size_t n_args, const mp_obj_t *pos_args, mp_map_t mp_obj_t res[2]; mp_float_t e; - res[0] = mp_obj_new_float(quad(fun, a, b, n, eps, &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_quad_obj, 2, integrate_quad); -#endif /* ULAB_INTEGRATE_HAS_QUAD */ +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 it expects an array of function values). +// 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); @@ -345,12 +344,12 @@ mp_float_t qromb(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint //| b: float, //| *, //| steps: int = 100 -//| eps: float = 1e-14, +//| eps: float = etolerance //| ) -> float: //| """ //| :param callable f: The function to integrate -//| :param float a: The left side of the interval -//| :param float b: The right side of the interval +//| :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 //| @@ -375,7 +374,7 @@ static mp_obj_t integrate_romberg(size_t n_args, const mp_obj_t *pos_args, mp_ma 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")); + mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a callable")); } // iterate over args 1, 2, and 4 @@ -441,12 +440,12 @@ mp_float_t qasi(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n //| b: float, //| *, //| steps: int = 100 -//| eps: float = 1e-14, +//| eps: float = etolerance //| ) -> float: //| """ //| :param callable f: The function to integrate -//| :param float a: The left side of the interval -//| :param float b: The right side of the interval +//| :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 //| @@ -497,7 +496,7 @@ static mp_obj_t integrate_simpson(size_t n_args, const mp_obj_t *pos_args, mp_ma MP_DEFINE_CONST_FUN_OBJ_KW(integrate_simpson_obj, 2, integrate_simpson); #endif /* ULAB_INTEGRATE_HAS_SIMPSON */ -#if ULAB_INTEGRATE_HAS_QUADGK +#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 @@ -608,18 +607,18 @@ mp_float_t qakro(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int } -//| def quadgk( +//| def quad( //| fun: Callable[[float], float], //| a: float, //| b: float, //| *, //| order: int = 5 -//| eps: float = 1e-14, +//| eps: float = etolerance //| ) -> float: //| """ //| :param callable f: The function to integrate -//| :param float a: The left side of the interval -//| :param float b: The right side of the interval +//| :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 //| @@ -629,7 +628,7 @@ mp_float_t qakro(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int //| ... //| -static mp_obj_t integrate_quadgk(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { +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 } }, @@ -644,7 +643,7 @@ static mp_obj_t integrate_quadgk(size_t n_args, const mp_obj_t *pos_args, mp_map 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")); + mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a callable")); } // iterate over args 1, 2, and 4 @@ -671,13 +670,13 @@ static mp_obj_t integrate_quadgk(size_t n_args, const mp_obj_t *pos_args, mp_map return mp_obj_new_tuple(2, res); } -MP_DEFINE_CONST_FUN_OBJ_KW(integrate_quadgk_obj, 2, integrate_quadgk); -#endif /* ULAB_INTEGRATE_HAS_QUADGK */ +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_QUAD - { MP_ROM_QSTR(MP_QSTR_quad), MP_ROM_PTR(&integrate_quad_obj) }, +#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) }, @@ -685,8 +684,8 @@ static const mp_rom_map_elem_t ulab_scipy_integrate_globals_table[] = { #if ULAB_INTEGRATE_HAS_SIMPSON { MP_ROM_QSTR(MP_QSTR_simpson), MP_ROM_PTR(&integrate_simpson_obj) }, #endif -#if ULAB_INTEGRATE_HAS_QUADGK - { MP_ROM_QSTR(MP_QSTR_quadgk), MP_ROM_PTR(&integrate_quadgk_obj) }, +#if ULAB_INTEGRATE_HAS_QUAD + { MP_ROM_QSTR(MP_QSTR_quad), MP_ROM_PTR(&integrate_quad_obj) }, #endif }; diff --git a/code/scipy/integrate/integrate.h b/code/scipy/integrate/integrate.h index 0fea59e8..ebfac2ea 100644 --- a/code/scipy/integrate/integrate.h +++ b/code/scipy/integrate/integrate.h @@ -17,8 +17,8 @@ extern const mp_obj_module_t ulab_scipy_integrate_module; -#if ULAB_INTEGRATE_HAS_QUAD -MP_DECLARE_CONST_FUN_OBJ_KW(optimize_quad_obj); +#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); @@ -26,8 +26,8 @@ MP_DECLARE_CONST_FUN_OBJ_KW(optimize_romberg_obj); #if ULAB_INTEGRATE_HAS_SIMPSON MP_DECLARE_CONST_FUN_OBJ_KW(optimize_simpson_obj); #endif -#if ULAB_INTEGRATE_HAS_QUADGK -MP_DECLARE_CONST_FUN_OBJ_KW(optimize_quadgk_obj); +#if ULAB_INTEGRATE_HAS_QUAD +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_quad_obj); #endif #endif /* _SCIPY_INTEGRATE_ */ diff --git a/code/ulab.h b/code/ulab.h index 22f32337..3eb30131 100644 --- a/code/ulab.h +++ b/code/ulab.h @@ -404,21 +404,20 @@ #define ULAB_SCIPY_HAS_INTEGRATE_MODULE (1) #endif -#ifndef ULAB_INTEGRATE_HAS_QUAD -#define ULAB_INTEGRATE_HAS_QUAD (1) +#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_QUADGK -#define ULAB_INTEGRATE_HAS_QUADGK (1) +#ifndef ULAB_INTEGRATE_HAS_QUAD +#define ULAB_INTEGRATE_HAS_QUAD (1) #endif // the linalg module; functions of the linalg module still have From 128defbecdd9995a2dae340623c55f58ef33a049 Mon Sep 17 00:00:00 2001 From: Harald Milz Date: Sun, 15 Dec 2024 14:16:15 +0100 Subject: [PATCH 22/22] major rewrite representing the name changes, mapping to CPython scipy.integrate, more background info --- docs/scipy-integrate.ipynb | 129 +++++++++++++++++++++---------------- 1 file changed, 73 insertions(+), 56 deletions(-) diff --git a/docs/scipy-integrate.ipynb b/docs/scipy-integrate.ipynb index ca0c50af..23220231 100644 --- a/docs/scipy-integrate.ipynb +++ b/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,30 +347,38 @@ "\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" ] }, @@ -379,47 +386,57 @@ "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": {