diff options
Diffstat (limited to 'circuitpython/extmod/ulab/code/scipy')
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/linalg/linalg.c | 280 | ||||
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/linalg/linalg.h | 21 | ||||
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/optimize/optimize.c | 415 | ||||
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/optimize/optimize.h | 41 | ||||
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/scipy.c | 52 | ||||
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/scipy.h | 21 | ||||
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/signal/signal.c | 172 | ||||
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/signal/signal.h | 24 | ||||
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/special/special.c | 43 | ||||
| -rw-r--r-- | circuitpython/extmod/ulab/code/scipy/special/special.h | 21 |
10 files changed, 1090 insertions, 0 deletions
diff --git a/circuitpython/extmod/ulab/code/scipy/linalg/linalg.c b/circuitpython/extmod/ulab/code/scipy/linalg/linalg.c new file mode 100644 index 0000000..d211f72 --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/linalg/linalg.c @@ -0,0 +1,280 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2021 Vikas Udupa + * +*/ + +#include <stdlib.h> +#include <string.h> +#include <math.h> +#include "py/obj.h" +#include "py/runtime.h" +#include "py/misc.h" + +#include "../../ulab.h" +#include "../../ulab_tools.h" +#include "../../numpy/linalg/linalg_tools.h" +#include "linalg.h" + +#if ULAB_SCIPY_HAS_LINALG_MODULE +//| +//| import ulab.scipy +//| import ulab.numpy +//| +//| """Linear algebra functions""" +//| + +#if ULAB_MAX_DIMS > 1 + +//| def solve_triangular(A: ulab.numpy.ndarray, b: ulab.numpy.ndarray, lower: bool) -> ulab.numpy.ndarray: +//| """ +//| :param ~ulab.numpy.ndarray A: a matrix +//| :param ~ulab.numpy.ndarray b: a vector +//| :param ~bool lower: if true, use only data contained in lower triangle of A, else use upper triangle of A +//| :return: solution to the system A x = b. Shape of return matches b +//| :raises TypeError: if A and b are not of type ndarray and are not dense +//| :raises ValueError: if A is a singular matrix +//| +//| Solve the equation A x = b for x, assuming A is a triangular matrix""" +//| ... +//| + +static mp_obj_t solve_triangular(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + + size_t i, j; + + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = mp_const_none} } , + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = mp_const_none} } , + { MP_QSTR_lower, MP_ARG_OBJ, { .u_rom_obj = mp_const_false } }, + }; + + 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); + + if(!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type) || !mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) { + mp_raise_TypeError(translate("first two arguments must be ndarrays")); + } + + ndarray_obj_t *A = MP_OBJ_TO_PTR(args[0].u_obj); + ndarray_obj_t *b = MP_OBJ_TO_PTR(args[1].u_obj); + + if(!ndarray_is_dense(A) || !ndarray_is_dense(b)) { + mp_raise_TypeError(translate("input must be a dense ndarray")); + } + + size_t A_rows = A->shape[ULAB_MAX_DIMS - 2]; + size_t A_cols = A->shape[ULAB_MAX_DIMS - 1]; + + uint8_t *A_arr = (uint8_t *)A->array; + uint8_t *b_arr = (uint8_t *)b->array; + + mp_float_t (*get_A_ele)(void *) = ndarray_get_float_function(A->dtype); + mp_float_t (*get_b_ele)(void *) = ndarray_get_float_function(b->dtype); + + uint8_t *temp_A = A_arr; + + // check if input matrix A is singular + for (i = 0; i < A_rows; i++) { + if (MICROPY_FLOAT_C_FUN(fabs)(get_A_ele(A_arr)) < LINALG_EPSILON) + mp_raise_ValueError(translate("input matrix is singular")); + A_arr += A->strides[ULAB_MAX_DIMS - 2]; + A_arr += A->strides[ULAB_MAX_DIMS - 1]; + } + + A_arr = temp_A; + + ndarray_obj_t *x = ndarray_new_dense_ndarray(b->ndim, b->shape, NDARRAY_FLOAT); + mp_float_t *x_arr = (mp_float_t *)x->array; + + if (mp_obj_is_true(args[2].u_obj)) { + // Solve the lower triangular matrix by iterating each row of A. + // Start by finding the first unknown using the first row. + // On finding this unknown, find the second unknown using the second row. + // Continue the same till the last unknown is found using the last row. + + for (i = 0; i < A_rows; i++) { + mp_float_t sum = 0.0; + for (j = 0; j < i; j++) { + sum += (get_A_ele(A_arr) * (*x_arr++)); + A_arr += A->strides[ULAB_MAX_DIMS - 1]; + } + + sum = (get_b_ele(b_arr) - sum) / (get_A_ele(A_arr)); + *x_arr = sum; + + x_arr -= j; + A_arr -= A->strides[ULAB_MAX_DIMS - 1] * j; + A_arr += A->strides[ULAB_MAX_DIMS - 2]; + b_arr += b->strides[ULAB_MAX_DIMS - 1]; + } + } else { + // Solve the upper triangular matrix by iterating each row of A. + // Start by finding the last unknown using the last row. + // On finding this unknown, find the last-but-one unknown using the last-but-one row. + // Continue the same till the first unknown is found using the first row. + + A_arr += (A->strides[ULAB_MAX_DIMS - 2] * A_rows); + b_arr += (b->strides[ULAB_MAX_DIMS - 1] * A_cols); + x_arr += A_cols; + + for (i = A_rows - 1; i < A_rows; i--) { + mp_float_t sum = 0.0; + for (j = i + 1; j < A_cols; j++) { + sum += (get_A_ele(A_arr) * (*x_arr++)); + A_arr += A->strides[ULAB_MAX_DIMS - 1]; + } + + x_arr -= (j - i); + A_arr -= (A->strides[ULAB_MAX_DIMS - 1] * (j - i)); + b_arr -= b->strides[ULAB_MAX_DIMS - 1]; + + sum = (get_b_ele(b_arr) - sum) / get_A_ele(A_arr); + *x_arr = sum; + + A_arr -= A->strides[ULAB_MAX_DIMS - 2]; + } + } + + return MP_OBJ_FROM_PTR(x); +} + +MP_DEFINE_CONST_FUN_OBJ_KW(linalg_solve_triangular_obj, 2, solve_triangular); + +//| def cho_solve(L: ulab.numpy.ndarray, b: ulab.numpy.ndarray) -> ulab.numpy.ndarray: +//| """ +//| :param ~ulab.numpy.ndarray L: the lower triangular, Cholesky factorization of A +//| :param ~ulab.numpy.ndarray b: right-hand-side vector b +//| :return: solution to the system A x = b. Shape of return matches b +//| :raises TypeError: if L and b are not of type ndarray and are not dense +//| +//| Solve the linear equations A x = b, given the Cholesky factorization of A as input""" +//| ... +//| + +static mp_obj_t cho_solve(mp_obj_t _L, mp_obj_t _b) { + + if(!mp_obj_is_type(_L, &ulab_ndarray_type) || !mp_obj_is_type(_b, &ulab_ndarray_type)) { + mp_raise_TypeError(translate("first two arguments must be ndarrays")); + } + + ndarray_obj_t *L = MP_OBJ_TO_PTR(_L); + ndarray_obj_t *b = MP_OBJ_TO_PTR(_b); + + if(!ndarray_is_dense(L) || !ndarray_is_dense(b)) { + mp_raise_TypeError(translate("input must be a dense ndarray")); + } + + mp_float_t (*get_L_ele)(void *) = ndarray_get_float_function(L->dtype); + mp_float_t (*get_b_ele)(void *) = ndarray_get_float_function(b->dtype); + void (*set_L_ele)(void *, mp_float_t) = ndarray_set_float_function(L->dtype); + + size_t L_rows = L->shape[ULAB_MAX_DIMS - 2]; + size_t L_cols = L->shape[ULAB_MAX_DIMS - 1]; + + // Obtain transpose of the input matrix L in L_t + size_t L_t_shape[ULAB_MAX_DIMS]; + size_t L_t_rows = L_t_shape[ULAB_MAX_DIMS - 2] = L_cols; + size_t L_t_cols = L_t_shape[ULAB_MAX_DIMS - 1] = L_rows; + ndarray_obj_t *L_t = ndarray_new_dense_ndarray(L->ndim, L_t_shape, L->dtype); + + uint8_t *L_arr = (uint8_t *)L->array; + uint8_t *L_t_arr = (uint8_t *)L_t->array; + uint8_t *b_arr = (uint8_t *)b->array; + + size_t i, j; + + uint8_t *L_ptr = L_arr; + uint8_t *L_t_ptr = L_t_arr; + for (i = 0; i < L_rows; i++) { + for (j = 0; j < L_cols; j++) { + set_L_ele(L_t_ptr, get_L_ele(L_ptr)); + L_t_ptr += L_t->strides[ULAB_MAX_DIMS - 2]; + L_ptr += L->strides[ULAB_MAX_DIMS - 1]; + } + + L_t_ptr -= j * L_t->strides[ULAB_MAX_DIMS - 2]; + L_t_ptr += L_t->strides[ULAB_MAX_DIMS - 1]; + L_ptr -= j * L->strides[ULAB_MAX_DIMS - 1]; + L_ptr += L->strides[ULAB_MAX_DIMS - 2]; + } + + ndarray_obj_t *x = ndarray_new_dense_ndarray(b->ndim, b->shape, NDARRAY_FLOAT); + mp_float_t *x_arr = (mp_float_t *)x->array; + + ndarray_obj_t *y = ndarray_new_dense_ndarray(b->ndim, b->shape, NDARRAY_FLOAT); + mp_float_t *y_arr = (mp_float_t *)y->array; + + // solve L y = b to obtain y, where L_t x = y + for (i = 0; i < L_rows; i++) { + mp_float_t sum = 0.0; + for (j = 0; j < i; j++) { + sum += (get_L_ele(L_arr) * (*y_arr++)); + L_arr += L->strides[ULAB_MAX_DIMS - 1]; + } + + sum = (get_b_ele(b_arr) - sum) / (get_L_ele(L_arr)); + *y_arr = sum; + + y_arr -= j; + L_arr -= L->strides[ULAB_MAX_DIMS - 1] * j; + L_arr += L->strides[ULAB_MAX_DIMS - 2]; + b_arr += b->strides[ULAB_MAX_DIMS - 1]; + } + + // using y, solve L_t x = y to obtain x + L_t_arr += (L_t->strides[ULAB_MAX_DIMS - 2] * L_t_rows); + y_arr += L_t_cols; + x_arr += L_t_cols; + + for (i = L_t_rows - 1; i < L_t_rows; i--) { + mp_float_t sum = 0.0; + for (j = i + 1; j < L_t_cols; j++) { + sum += (get_L_ele(L_t_arr) * (*x_arr++)); + L_t_arr += L_t->strides[ULAB_MAX_DIMS - 1]; + } + + x_arr -= (j - i); + L_t_arr -= (L_t->strides[ULAB_MAX_DIMS - 1] * (j - i)); + y_arr--; + + sum = ((*y_arr) - sum) / get_L_ele(L_t_arr); + *x_arr = sum; + + L_t_arr -= L_t->strides[ULAB_MAX_DIMS - 2]; + } + + return MP_OBJ_FROM_PTR(x); +} + +MP_DEFINE_CONST_FUN_OBJ_2(linalg_cho_solve_obj, cho_solve); + +#endif + +static const mp_rom_map_elem_t ulab_scipy_linalg_globals_table[] = { + { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_linalg) }, + #if ULAB_MAX_DIMS > 1 + #if ULAB_SCIPY_LINALG_HAS_SOLVE_TRIANGULAR + { MP_ROM_QSTR(MP_QSTR_solve_triangular), (mp_obj_t)&linalg_solve_triangular_obj }, + #endif + #if ULAB_SCIPY_LINALG_HAS_CHO_SOLVE + { MP_ROM_QSTR(MP_QSTR_cho_solve), (mp_obj_t)&linalg_cho_solve_obj }, + #endif + #endif +}; + +static MP_DEFINE_CONST_DICT(mp_module_ulab_scipy_linalg_globals, ulab_scipy_linalg_globals_table); + +const mp_obj_module_t ulab_scipy_linalg_module = { + .base = { &mp_type_module }, + .globals = (mp_obj_dict_t*)&mp_module_ulab_scipy_linalg_globals, +}; +MP_REGISTER_MODULE(MP_QSTR_ulab_dot_scipy_dot_linalg, ulab_scipy_linalg_module, MODULE_ULAB_ENABLED && CIRCUITPY_ULAB); + +#endif diff --git a/circuitpython/extmod/ulab/code/scipy/linalg/linalg.h b/circuitpython/extmod/ulab/code/scipy/linalg/linalg.h new file mode 100644 index 0000000..628051f --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/linalg/linalg.h @@ -0,0 +1,21 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2021 Vikas Udupa + * +*/ + +#ifndef _SCIPY_LINALG_ +#define _SCIPY_LINALG_ + +extern const mp_obj_module_t ulab_scipy_linalg_module; + +MP_DECLARE_CONST_FUN_OBJ_KW(linalg_solve_triangular_obj); +MP_DECLARE_CONST_FUN_OBJ_2(linalg_cho_solve_obj); + +#endif /* _SCIPY_LINALG_ */ diff --git a/circuitpython/extmod/ulab/code/scipy/optimize/optimize.c b/circuitpython/extmod/ulab/code/scipy/optimize/optimize.c new file mode 100644 index 0000000..f1c746a --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/optimize/optimize.c @@ -0,0 +1,415 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2020 Jeff Epler for Adafruit Industries + * 2020 Scott Shawcroft for Adafruit Industries + * 2020-2021 Zoltán Vörös + * 2020 Taku Fukada +*/ + +#include <math.h> +#include "py/obj.h" +#include "py/runtime.h" +#include "py/misc.h" + +#include "../../ndarray.h" +#include "../../ulab.h" +#include "../../ulab_tools.h" +#include "optimize.h" + +const mp_obj_float_t xtolerance = {{&mp_type_float}, MICROPY_FLOAT_CONST(2.4e-7)}; +const mp_obj_float_t rtolerance = {{&mp_type_float}, MICROPY_FLOAT_CONST(0.0)}; + +static mp_float_t optimize_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(type->MP_TYPE_CALL(fun, nparams+1, 0, fargs)); +} + +#if ULAB_SCIPY_OPTIMIZE_HAS_BISECT +//| def bisect( +//| fun: Callable[[float], float], +//| a: float, +//| b: float, +//| *, +//| xtol: float = 2.4e-7, +//| maxiter: int = 100 +//| ) -> float: +//| """ +//| :param callable f: The function to bisect +//| :param float a: The left side of the interval +//| :param float b: The right side of the interval +//| :param float xtol: The tolerance value +//| :param float maxiter: The maximum number of iterations to perform +//| +//| 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 optimize_bisect(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + // Simple bisection routine + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_xtol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} }, + { MP_QSTR_maxiter, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 100} }, + }; + + 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_type_get_call_slot(type) == NULL) { + mp_raise_TypeError(translate("first argument must be a function")); + } + mp_float_t xtol = mp_obj_get_float(args[3].u_obj); + mp_obj_t *fargs = m_new(mp_obj_t, 1); + mp_float_t left, right; + mp_float_t x_mid; + mp_float_t a = mp_obj_get_float(args[1].u_obj); + mp_float_t b = mp_obj_get_float(args[2].u_obj); + left = optimize_python_call(type, fun, a, fargs, 0); + right = optimize_python_call(type, fun, b, fargs, 0); + if(left * right > 0) { + mp_raise_ValueError(translate("function has the same sign at the ends of interval")); + } + mp_float_t rtb = left < MICROPY_FLOAT_CONST(0.0) ? a : b; + mp_float_t dx = left < MICROPY_FLOAT_CONST(0.0) ? b - a : a - b; + if(args[4].u_int < 0) { + mp_raise_ValueError(translate("maxiter should be > 0")); + } + for(uint16_t i=0; i < args[4].u_int; i++) { + dx *= MICROPY_FLOAT_CONST(0.5); + x_mid = rtb + dx; + if(optimize_python_call(type, fun, x_mid, fargs, 0) < MICROPY_FLOAT_CONST(0.0)) { + rtb = x_mid; + } + if(MICROPY_FLOAT_C_FUN(fabs)(dx) < xtol) break; + } + return mp_obj_new_float(rtb); +} + +MP_DEFINE_CONST_FUN_OBJ_KW(optimize_bisect_obj, 3, optimize_bisect); +#endif + +#if ULAB_SCIPY_OPTIMIZE_HAS_FMIN +//| def fmin( +//| fun: Callable[[float], float], +//| x0: float, +//| *, +//| xatol: float = 2.4e-7, +//| fatol: float = 2.4e-7, +//| maxiter: int = 200 +//| ) -> float: +//| """ +//| :param callable f: The function to bisect +//| :param float x0: The initial x value +//| :param float xatol: The absolute tolerance value +//| :param float fatol: The relative tolerance value +//| +//| Find a minimum of the function ``f(x)`` using the downhill simplex method. +//| The located ``x`` is within ``fxtol`` of the actual minimum, and ``f(x)`` +//| is within ``fatol`` of the actual minimum unless more than ``maxiter`` +//| steps are requried.""" +//| ... +//| + +STATIC mp_obj_t optimize_fmin(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + // downhill simplex method in 1D + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_xatol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} }, + { MP_QSTR_fatol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} }, + { MP_QSTR_maxiter, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 200} }, + }; + + 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_type_get_call_slot(type) == NULL) { + mp_raise_TypeError(translate("first argument must be a function")); + } + + // parameters controlling convergence conditions + mp_float_t xatol = mp_obj_get_float(args[2].u_obj); + mp_float_t fatol = mp_obj_get_float(args[3].u_obj); + if(args[4].u_int <= 0) { + mp_raise_ValueError(translate("maxiter must be > 0")); + } + uint16_t maxiter = (uint16_t)args[4].u_int; + + mp_float_t x0 = mp_obj_get_float(args[1].u_obj); + mp_float_t x1 = MICROPY_FLOAT_C_FUN(fabs)(x0) > OPTIMIZE_EPSILON ? (MICROPY_FLOAT_CONST(1.0) + OPTIMIZE_NONZDELTA) * x0 : OPTIMIZE_ZDELTA; + mp_obj_t *fargs = m_new(mp_obj_t, 1); + mp_float_t f0 = optimize_python_call(type, fun, x0, fargs, 0); + mp_float_t f1 = optimize_python_call(type, fun, x1, fargs, 0); + if(f1 < f0) { + SWAP(mp_float_t, x0, x1); + SWAP(mp_float_t, f0, f1); + } + for(uint16_t i=0; i < maxiter; i++) { + uint8_t shrink = 0; + f0 = optimize_python_call(type, fun, x0, fargs, 0); + f1 = optimize_python_call(type, fun, x1, fargs, 0); + + // reflection + mp_float_t xr = (MICROPY_FLOAT_CONST(1.0) + OPTIMIZE_ALPHA) * x0 - OPTIMIZE_ALPHA * x1; + mp_float_t fr = optimize_python_call(type, fun, xr, fargs, 0); + if(fr < f0) { // expansion + mp_float_t xe = (1 + OPTIMIZE_ALPHA * OPTIMIZE_BETA) * x0 - OPTIMIZE_ALPHA * OPTIMIZE_BETA * x1; + mp_float_t fe = optimize_python_call(type, fun, xe, fargs, 0); + if(fe < fr) { + x1 = xe; + f1 = fe; + } else { + x1 = xr; + f1 = fr; + } + } else { + if(fr < f1) { // contraction + mp_float_t xc = (1 + OPTIMIZE_GAMMA * OPTIMIZE_ALPHA) * x0 - OPTIMIZE_GAMMA * OPTIMIZE_ALPHA * x1; + mp_float_t fc = optimize_python_call(type, fun, xc, fargs, 0); + if(fc < fr) { + x1 = xc; + f1 = fc; + } else { + shrink = 1; + } + } else { // inside contraction + mp_float_t xc = (MICROPY_FLOAT_CONST(1.0) - OPTIMIZE_GAMMA) * x0 + OPTIMIZE_GAMMA * x1; + mp_float_t fc = optimize_python_call(type, fun, xc, fargs, 0); + if(fc < f1) { + x1 = xc; + f1 = fc; + } else { + shrink = 1; + } + } + if(shrink == 1) { + x1 = x0 + OPTIMIZE_DELTA * (x1 - x0); + f1 = optimize_python_call(type, fun, x1, fargs, 0); + } + if((MICROPY_FLOAT_C_FUN(fabs)(f1 - f0) < fatol) || + (MICROPY_FLOAT_C_FUN(fabs)(x1 - x0) < xatol)) { + break; + } + if(f1 < f0) { + SWAP(mp_float_t, x0, x1); + SWAP(mp_float_t, f0, f1); + } + } + } + return mp_obj_new_float(x0); +} + +MP_DEFINE_CONST_FUN_OBJ_KW(optimize_fmin_obj, 2, optimize_fmin); +#endif + +#if ULAB_SCIPY_OPTIMIZE_HAS_CURVE_FIT +static void optimize_jacobi(const mp_obj_type_t *type, mp_obj_t fun, mp_float_t *x, mp_float_t *y, uint16_t len, mp_float_t *params, uint8_t nparams, mp_float_t *jacobi, mp_float_t *grad) { + /* Calculates the Jacobian and the gradient of the cost function + * + * The entries in the Jacobian are + * J(m, n) = de_m/da_n, + * + * where + * + * e_m = (f(x_m, a1, a2, ...) - y_m)/sigma_m is the error at x_m, + * + * and + * + * a1, a2, ..., a_n are the free parameters + */ + mp_obj_t *fargs0 = m_new(mp_obj_t, lenp+1); + mp_obj_t *fargs1 = m_new(mp_obj_t, lenp+1); + for(uint8_t p=0; p < nparams; p++) { + fargs0[p+1] = mp_obj_new_float(params[p]); + fargs1[p+1] = mp_obj_new_float(params[p]); + } + for(uint8_t p=0; p < nparams; p++) { + mp_float_t da = params[p] != MICROPY_FLOAT_CONST(0.0) ? (MICROPY_FLOAT_CONST(1.0) + APPROX_NONZDELTA) * params[p] : APPROX_ZDELTA; + fargs1[p+1] = mp_obj_new_float(params[p] + da); + grad[p] = MICROPY_FLOAT_CONST(0.0); + for(uint16_t i=0; i < len; i++) { + mp_float_t f0 = optimize_python_call(type, fun, x[i], fargs0, nparams); + mp_float_t f1 = optimize_python_call(type, fun, x[i], fargs1, nparams); + jacobi[i*nparamp+p] = (f1 - f0) / da; + grad[p] += (f0 - y[i]) * jacobi[i*nparamp+p]; + } + fargs1[p+1] = fargs0[p+1]; // set back to the original value + } +} + +static void optimize_delta(mp_float_t *jacobi, mp_float_t *grad, uint16_t len, uint8_t nparams, mp_float_t lambda) { + // +} + +mp_obj_t optimize_curve_fit(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + // Levenberg-Marquardt non-linear fit + // The implementation follows the introductory discussion in Mark Tanstrum's paper, https://arxiv.org/abs/1201.5885 + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_p0, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_xatol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} }, + { MP_QSTR_fatol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} }, + { MP_QSTR_maxiter, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} }, + }; + + 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_type_get_call_slot(type) == NULL) { + mp_raise_TypeError(translate("first argument must be a function")); + } + + mp_obj_t x_obj = args[1].u_obj; + mp_obj_t y_obj = args[2].u_obj; + mp_obj_t p0_obj = args[3].u_obj; + if(!ndarray_object_is_array_like(x_obj) || !ndarray_object_is_array_like(y_obj)) { + mp_raise_TypeError(translate("data must be iterable")); + } + if(!ndarray_object_is_nditerable(p0_obj)) { + mp_raise_TypeError(translate("initial values must be iterable")); + } + size_t len = (size_t)mp_obj_get_int(mp_obj_len_maybe(x_obj)); + uint8_t lenp = (uint8_t)mp_obj_get_int(mp_obj_len_maybe(p0_obj)); + if(len != (uint16_t)mp_obj_get_int(mp_obj_len_maybe(y_obj))) { + mp_raise_ValueError(translate("data must be of equal length")); + } + + mp_float_t *x = m_new(mp_float_t, len); + fill_array_iterable(x, x_obj); + mp_float_t *y = m_new(mp_float_t, len); + fill_array_iterable(y, y_obj); + mp_float_t *p0 = m_new(mp_float_t, lenp); + fill_array_iterable(p0, p0_obj); + mp_float_t *grad = m_new(mp_float_t, len); + mp_float_t *jacobi = m_new(mp_float_t, len*len); + mp_obj_t *fargs = m_new(mp_obj_t, lenp+1); + + m_del(mp_float_t, p0, lenp); + // parameters controlling convergence conditions + //mp_float_t xatol = mp_obj_get_float(args[2].u_obj); + //mp_float_t fatol = mp_obj_get_float(args[3].u_obj); + + // this has finite binary representation; we will multiply/divide by 4 + //mp_float_t lambda = 0.0078125; + + //linalg_invert_matrix(mp_float_t *data, size_t N) + + m_del(mp_float_t, x, len); + m_del(mp_float_t, y, len); + m_del(mp_float_t, grad, len); + m_del(mp_float_t, jacobi, len*len); + m_del(mp_obj_t, fargs, lenp+1); + return mp_const_none; +} + +MP_DEFINE_CONST_FUN_OBJ_KW(optimize_curve_fit_obj, 2, optimize_curve_fit); +#endif + +#if ULAB_SCIPY_OPTIMIZE_HAS_NEWTON +//| def newton( +//| fun: Callable[[float], float], +//| x0: float, +//| *, +//| xtol: float = 2.4e-7, +//| rtol: float = 0.0, +//| maxiter: int = 50 +//| ) -> float: +//| """ +//| :param callable f: The function to bisect +//| :param float x0: The initial x value +//| :param float xtol: The absolute tolerance value +//| :param float rtol: The relative tolerance value +//| :param float maxiter: The maximum number of iterations to perform +//| +//| Find a solution (zero) of the function ``f(x)`` using Newton's Method. +//| The result is accurate to within ``xtol * rtol * |f(x)|`` unless more than +//| ``maxiter`` steps are requried.""" +//| ... +//| + +static mp_obj_t optimize_newton(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + // this is actually the secant method, as the first derivative of the function + // is not accepted as an argument. The function whose root we want to solve for + // must depend on a single variable without parameters, i.e., f(x) + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } }, + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } }, + { MP_QSTR_tol, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_PTR(&xtolerance) } }, + { MP_QSTR_rtol, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_PTR(&rtolerance) } }, + { MP_QSTR_maxiter, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = 50 } }, + }; + + 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_type_get_call_slot(type) == NULL) { + mp_raise_TypeError(translate("first argument must be a function")); + } + mp_float_t x = mp_obj_get_float(args[1].u_obj); + mp_float_t tol = mp_obj_get_float(args[2].u_obj); + mp_float_t rtol = mp_obj_get_float(args[3].u_obj); + mp_float_t dx, df, fx; + dx = x > MICROPY_FLOAT_CONST(0.0) ? OPTIMIZE_EPS * x : -OPTIMIZE_EPS * x; + mp_obj_t *fargs = m_new(mp_obj_t, 1); + if(args[4].u_int <= 0) { + mp_raise_ValueError(translate("maxiter must be > 0")); + } + for(uint16_t i=0; i < args[4].u_int; i++) { + fx = optimize_python_call(type, fun, x, fargs, 0); + df = (optimize_python_call(type, fun, x + dx, fargs, 0) - fx) / dx; + dx = fx / df; + x -= dx; + if(MICROPY_FLOAT_C_FUN(fabs)(dx) < (tol + rtol * MICROPY_FLOAT_C_FUN(fabs)(x))) break; + } + return mp_obj_new_float(x); +} + +MP_DEFINE_CONST_FUN_OBJ_KW(optimize_newton_obj, 2, optimize_newton); +#endif + +static const mp_rom_map_elem_t ulab_scipy_optimize_globals_table[] = { + { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_optimize) }, + #if ULAB_SCIPY_OPTIMIZE_HAS_BISECT + { MP_OBJ_NEW_QSTR(MP_QSTR_bisect), (mp_obj_t)&optimize_bisect_obj }, + #endif + #if ULAB_SCIPY_OPTIMIZE_HAS_CURVE_FIT + { MP_OBJ_NEW_QSTR(MP_QSTR_curve_fit), (mp_obj_t)&optimize_curve_fit_obj }, + #endif + #if ULAB_SCIPY_OPTIMIZE_HAS_FMIN + { MP_OBJ_NEW_QSTR(MP_QSTR_fmin), (mp_obj_t)&optimize_fmin_obj }, + #endif + #if ULAB_SCIPY_OPTIMIZE_HAS_NEWTON + { MP_OBJ_NEW_QSTR(MP_QSTR_newton), (mp_obj_t)&optimize_newton_obj }, + #endif +}; + +static MP_DEFINE_CONST_DICT(mp_module_ulab_scipy_optimize_globals, ulab_scipy_optimize_globals_table); + +const mp_obj_module_t ulab_scipy_optimize_module = { + .base = { &mp_type_module }, + .globals = (mp_obj_dict_t*)&mp_module_ulab_scipy_optimize_globals, +}; +MP_REGISTER_MODULE(MP_QSTR_ulab_dot_scipy_dot_optimize, ulab_scipy_optimize_module, MODULE_ULAB_ENABLED && CIRCUITPY_ULAB); diff --git a/circuitpython/extmod/ulab/code/scipy/optimize/optimize.h b/circuitpython/extmod/ulab/code/scipy/optimize/optimize.h new file mode 100644 index 0000000..174b386 --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/optimize/optimize.h @@ -0,0 +1,41 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2020-2021 Zoltán Vörös + * +*/ + +#ifndef _SCIPY_OPTIMIZE_ +#define _SCIPY_OPTIMIZE_ + +#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_optimize_module; + +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_bisect_obj); +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_curve_fit_obj); +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_fmin_obj); +MP_DECLARE_CONST_FUN_OBJ_KW(optimize_newton_obj); + +#endif /* _SCIPY_OPTIMIZE_ */ diff --git a/circuitpython/extmod/ulab/code/scipy/scipy.c b/circuitpython/extmod/ulab/code/scipy/scipy.c new file mode 100644 index 0000000..ba37dde --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/scipy.c @@ -0,0 +1,52 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2020 Jeff Epler for Adafruit Industries + * 2020 Scott Shawcroft for Adafruit Industries + * 2020-2021 Zoltán Vörös + * 2020 Taku Fukada +*/ + +#include <math.h> +#include "py/runtime.h" + +#include "../ulab.h" +#include "optimize/optimize.h" +#include "signal/signal.h" +#include "special/special.h" +#include "linalg/linalg.h" + +#if ULAB_HAS_SCIPY + +//| """Compatibility layer for scipy""" +//| + +static const mp_rom_map_elem_t ulab_scipy_globals_table[] = { + { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_scipy) }, + #if ULAB_SCIPY_HAS_LINALG_MODULE + { MP_ROM_QSTR(MP_QSTR_linalg), MP_ROM_PTR(&ulab_scipy_linalg_module) }, + #endif + #if ULAB_SCIPY_HAS_OPTIMIZE_MODULE + { MP_ROM_QSTR(MP_QSTR_optimize), MP_ROM_PTR(&ulab_scipy_optimize_module) }, + #endif + #if ULAB_SCIPY_HAS_SIGNAL_MODULE + { MP_ROM_QSTR(MP_QSTR_signal), MP_ROM_PTR(&ulab_scipy_signal_module) }, + #endif + #if ULAB_SCIPY_HAS_SPECIAL_MODULE + { MP_ROM_QSTR(MP_QSTR_special), MP_ROM_PTR(&ulab_scipy_special_module) }, + #endif +}; + +static MP_DEFINE_CONST_DICT(mp_module_ulab_scipy_globals, ulab_scipy_globals_table); + +const mp_obj_module_t ulab_scipy_module = { + .base = { &mp_type_module }, + .globals = (mp_obj_dict_t*)&mp_module_ulab_scipy_globals, +}; +MP_REGISTER_MODULE(MP_QSTR_ulab_dot_scipy, ulab_scipy_module, MODULE_ULAB_ENABLED && CIRCUITPY_ULAB); +#endif diff --git a/circuitpython/extmod/ulab/code/scipy/scipy.h b/circuitpython/extmod/ulab/code/scipy/scipy.h new file mode 100644 index 0000000..ec8c804 --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/scipy.h @@ -0,0 +1,21 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2020-2021 Zoltán Vörös + * +*/ + +#ifndef _SCIPY_ +#define _SCIPY_ + +#include "../ulab.h" +#include "../ndarray.h" + +extern const mp_obj_module_t ulab_scipy_module; + +#endif /* _SCIPY_ */ diff --git a/circuitpython/extmod/ulab/code/scipy/signal/signal.c b/circuitpython/extmod/ulab/code/scipy/signal/signal.c new file mode 100644 index 0000000..69d5609 --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/signal/signal.c @@ -0,0 +1,172 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2020 Jeff Epler for Adafruit Industries + * 2020 Scott Shawcroft for Adafruit Industries + * 2020-2021 Zoltán Vörös + * 2020 Taku Fukada +*/ + +#include <math.h> +#include <string.h> +#include "py/runtime.h" + +#include "../../ulab.h" +#include "../../ndarray.h" +#include "../../numpy/carray/carray_tools.h" +#include "../../numpy/fft/fft_tools.h" + +#if ULAB_SCIPY_SIGNAL_HAS_SPECTROGRAM +//| import ulab.numpy +//| +//| def spectrogram(r: ulab.numpy.ndarray) -> ulab.numpy.ndarray: +//| """ +//| :param ulab.numpy.ndarray r: A 1-dimension array of values whose size is a power of 2 +//| +//| Computes the spectrum of the input signal. This is the absolute value of the (complex-valued) fft of the signal. +//| This function is similar to scipy's ``scipy.signal.spectrogram``.""" +//| ... +//| + +mp_obj_t signal_spectrogram(size_t n_args, const mp_obj_t *args) { + #if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE + return fft_fft_ifft_spectrogram(args[0], FFT_SPECTROGRAM); + #else + if(n_args == 2) { + return fft_fft_ifft_spectrogram(n_args, args[0], args[1], FFT_SPECTROGRAM); + } else { + return fft_fft_ifft_spectrogram(n_args, args[0], mp_const_none, FFT_SPECTROGRAM); + } + #endif +} + +#if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE +MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(signal_spectrogram_obj, 1, 1, signal_spectrogram); +#else +MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(signal_spectrogram_obj, 1, 2, signal_spectrogram); +#endif + +#endif /* ULAB_SCIPY_SIGNAL_HAS_SPECTROGRAM */ + +#if ULAB_SCIPY_SIGNAL_HAS_SOSFILT +static void signal_sosfilt_array(mp_float_t *x, const mp_float_t *coeffs, mp_float_t *zf, const size_t len) { + for(size_t i=0; i < len; i++) { + mp_float_t xn = *x; + *x = coeffs[0] * xn + zf[0]; + zf[0] = zf[1] + coeffs[1] * xn - coeffs[4] * *x; + zf[1] = coeffs[2] * xn - coeffs[5] * *x; + x++; + } + x -= len; +} + +mp_obj_t signal_sosfilt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + static const mp_arg_t allowed_args[] = { + { MP_QSTR_sos, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_x, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + { MP_QSTR_zi, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } }, + }; + + 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); + + if(!ndarray_object_is_array_like(args[0].u_obj) || !ndarray_object_is_array_like(args[1].u_obj)) { + mp_raise_TypeError(translate("sosfilt requires iterable arguments")); + } + #if ULAB_SUPPORTS_COMPLEX + if(mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) { + ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[1].u_obj); + COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype) + } + #endif + size_t lenx = (size_t)mp_obj_get_int(mp_obj_len_maybe(args[1].u_obj)); + ndarray_obj_t *y = ndarray_new_linear_array(lenx, NDARRAY_FLOAT); + mp_float_t *yarray = (mp_float_t *)y->array; + mp_float_t coeffs[6]; + if(mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) { + ndarray_obj_t *inarray = MP_OBJ_TO_PTR(args[1].u_obj); + #if ULAB_MAX_DIMS > 1 + if(inarray->ndim > 1) { + mp_raise_ValueError(translate("input must be one-dimensional")); + } + #endif + uint8_t *iarray = (uint8_t *)inarray->array; + for(size_t i=0; i < lenx; i++) { + *yarray++ = ndarray_get_float_value(iarray, inarray->dtype); + iarray += inarray->strides[ULAB_MAX_DIMS - 1]; + } + yarray -= lenx; + } else { + fill_array_iterable(yarray, args[1].u_obj); + } + + mp_obj_iter_buf_t iter_buf; + mp_obj_t item, iterable = mp_getiter(args[0].u_obj, &iter_buf); + size_t lensos = (size_t)mp_obj_get_int(mp_obj_len_maybe(args[0].u_obj)); + + size_t *shape = ndarray_shape_vector(0, 0, lensos, 2); + ndarray_obj_t *zf = ndarray_new_dense_ndarray(2, shape, NDARRAY_FLOAT); + mp_float_t *zf_array = (mp_float_t *)zf->array; + + if(args[2].u_obj != mp_const_none) { + if(!mp_obj_is_type(args[2].u_obj, &ulab_ndarray_type)) { + mp_raise_TypeError(translate("zi must be an ndarray")); + } else { + ndarray_obj_t *zi = MP_OBJ_TO_PTR(args[2].u_obj); + if((zi->shape[ULAB_MAX_DIMS - 1] != lensos) || (zi->shape[ULAB_MAX_DIMS - 1] != 2)) { + mp_raise_ValueError(translate("zi must be of shape (n_section, 2)")); + } + if(zi->dtype != NDARRAY_FLOAT) { + mp_raise_ValueError(translate("zi must be of float type")); + } + // TODO: this won't work with sparse arrays + memcpy(zf_array, zi->array, 2*lensos*sizeof(mp_float_t)); + } + } + while((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) { + if(mp_obj_get_int(mp_obj_len_maybe(item)) != 6) { + mp_raise_ValueError(translate("sos array must be of shape (n_section, 6)")); + } else { + fill_array_iterable(coeffs, item); + if(coeffs[3] != MICROPY_FLOAT_CONST(1.0)) { + mp_raise_ValueError(translate("sos[:, 3] should be all ones")); + } + signal_sosfilt_array(yarray, coeffs, zf_array, lenx); + zf_array += 2; + } + } + if(args[2].u_obj == mp_const_none) { + return MP_OBJ_FROM_PTR(y); + } else { + mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(mp_obj_new_tuple(2, NULL)); + tuple->items[0] = MP_OBJ_FROM_PTR(y); + tuple->items[1] = MP_OBJ_FROM_PTR(zf); + return tuple; + } +} + +MP_DEFINE_CONST_FUN_OBJ_KW(signal_sosfilt_obj, 2, signal_sosfilt); +#endif /* ULAB_SCIPY_SIGNAL_HAS_SOSFILT */ + +static const mp_rom_map_elem_t ulab_scipy_signal_globals_table[] = { + { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_signal) }, + #if ULAB_SCIPY_SIGNAL_HAS_SPECTROGRAM + { MP_OBJ_NEW_QSTR(MP_QSTR_spectrogram), (mp_obj_t)&signal_spectrogram_obj }, + #endif + #if ULAB_SCIPY_SIGNAL_HAS_SOSFILT + { MP_OBJ_NEW_QSTR(MP_QSTR_sosfilt), (mp_obj_t)&signal_sosfilt_obj }, + #endif +}; + +static MP_DEFINE_CONST_DICT(mp_module_ulab_scipy_signal_globals, ulab_scipy_signal_globals_table); + +const mp_obj_module_t ulab_scipy_signal_module = { + .base = { &mp_type_module }, + .globals = (mp_obj_dict_t*)&mp_module_ulab_scipy_signal_globals, +}; +MP_REGISTER_MODULE(MP_QSTR_ulab_dot_scipy_dot_signal, ulab_scipy_signal_module, MODULE_ULAB_ENABLED && CIRCUITPY_ULAB); diff --git a/circuitpython/extmod/ulab/code/scipy/signal/signal.h b/circuitpython/extmod/ulab/code/scipy/signal/signal.h new file mode 100644 index 0000000..21299a6 --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/signal/signal.h @@ -0,0 +1,24 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2020-2021 Zoltán Vörös + * +*/ + +#ifndef _SCIPY_SIGNAL_ +#define _SCIPY_SIGNAL_ + +#include "../../ulab.h" +#include "../../ndarray.h" + +extern const mp_obj_module_t ulab_scipy_signal_module; + +MP_DECLARE_CONST_FUN_OBJ_VAR_BETWEEN(signal_spectrogram_obj); +MP_DECLARE_CONST_FUN_OBJ_KW(signal_sosfilt_obj); + +#endif /* _SCIPY_SIGNAL_ */ diff --git a/circuitpython/extmod/ulab/code/scipy/special/special.c b/circuitpython/extmod/ulab/code/scipy/special/special.c new file mode 100644 index 0000000..decfde0 --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/special/special.c @@ -0,0 +1,43 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2020 Jeff Epler for Adafruit Industries + * 2020 Scott Shawcroft for Adafruit Industries + * 2020-2021 Zoltán Vörös + * 2020 Taku Fukada +*/ + +#include <math.h> +#include "py/runtime.h" + +#include "../../ulab.h" +#include "../../numpy/vector.h" + +static const mp_rom_map_elem_t ulab_scipy_special_globals_table[] = { + { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_special) }, + #if ULAB_SCIPY_SPECIAL_HAS_ERF + { MP_OBJ_NEW_QSTR(MP_QSTR_erf), (mp_obj_t)&vector_erf_obj }, + #endif + #if ULAB_SCIPY_SPECIAL_HAS_ERFC + { MP_OBJ_NEW_QSTR(MP_QSTR_erfc), (mp_obj_t)&vector_erfc_obj }, + #endif + #if ULAB_SCIPY_SPECIAL_HAS_GAMMA + { MP_OBJ_NEW_QSTR(MP_QSTR_gamma), (mp_obj_t)&vector_gamma_obj }, + #endif + #if ULAB_SCIPY_SPECIAL_HAS_GAMMALN + { MP_OBJ_NEW_QSTR(MP_QSTR_gammaln), (mp_obj_t)&vector_lgamma_obj }, + #endif +}; + +static MP_DEFINE_CONST_DICT(mp_module_ulab_scipy_special_globals, ulab_scipy_special_globals_table); + +const mp_obj_module_t ulab_scipy_special_module = { + .base = { &mp_type_module }, + .globals = (mp_obj_dict_t*)&mp_module_ulab_scipy_special_globals, +}; +MP_REGISTER_MODULE(MP_QSTR_ulab_dot_scipy_dot_special, ulab_scipy_special_module, MODULE_ULAB_ENABLED && CIRCUITPY_ULAB); diff --git a/circuitpython/extmod/ulab/code/scipy/special/special.h b/circuitpython/extmod/ulab/code/scipy/special/special.h new file mode 100644 index 0000000..bb34e27 --- /dev/null +++ b/circuitpython/extmod/ulab/code/scipy/special/special.h @@ -0,0 +1,21 @@ + +/* + * This file is part of the micropython-ulab project, + * + * https://github.com/v923z/micropython-ulab + * + * The MIT License (MIT) + * + * Copyright (c) 2020-2021 Zoltán Vörös + * +*/ + +#ifndef _SCIPY_SPECIAL_ +#define _SCIPY_SPECIAL_ + +#include "../../ulab.h" +#include "../../ndarray.h" + +extern const mp_obj_module_t ulab_scipy_special_module; + +#endif /* _SCIPY_SPECIAL_ */ |