UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-03-18 21:30:38 +01:00
Code Issues Packages Projects Releases Wiki Activity

[python] Upgrade MicroPython to version 1.19.1, Ulab and add unit tests for Ulab (#259)

Browse Source
This commit is contained in:
Yaya-Cout
2022-06-26 13:03:22 +02:00
committed by GitHub
parent 25d4793441
commit 92c653f2f2
177 changed files with 8809 additions and 3819 deletions
Download Patch File Download Diff File Expand all files Collapse all files

115
python/port/genhdr/moduledefs.h
View File

@@ -1,13 +1,112 @@
// Automatically generated by makemoduledefs.py.
/* In the standard MicroPython build system, this file is autogenerated from the
* reset of the sources. We manually include it here some modules are not included
* by the build system, so we need to manually update the MicroPython part
*
* How to update this file with a new MicroPython release
* - Get a clean copy of MicroPython
* - Copy our mpconfigport.h over the "bare-arm" port of MicroPython
* - "make" the bare-arm port of MicroPython (don't worry if it doesn't finish)
* - "cat build/genhdr/moduledefs.h".
* - Insert the result below in the MicroPython section,
* until the definition of MICROPY_REGISTERED_MODULES
* - copy the MICROPY_REGISTERED_MODULES section at the end of this file,
* /!\ this section is present twice in the file, so you need to copy it twice
* Keep the Upsilon part when copying the MICROPY_REGISTERED_MODULES section
*/
#if (MICROPY_PY_ARRAY)
extern const struct _mp_obj_module_t mp_module_uarray;
#define MODULE_DEF_MP_QSTR_UARRAY { MP_ROM_QSTR(MP_QSTR_uarray), MP_ROM_PTR(&mp_module_uarray) },
#else
#define MODULE_DEF_MP_QSTR_UARRAY
#endif
// MicroPython part
extern const struct _mp_obj_module_t mp_module___main__;
#undef MODULE_DEF_MP_QSTR___MAIN__
#define MODULE_DEF_MP_QSTR___MAIN__ { MP_ROM_QSTR(MP_QSTR___main__), MP_ROM_PTR(&mp_module___main__) },
extern const struct _mp_obj_module_t mp_module_builtins;
#undef MODULE_DEF_MP_QSTR_BUILTINS
#define MODULE_DEF_MP_QSTR_BUILTINS { MP_ROM_QSTR(MP_QSTR_builtins), MP_ROM_PTR(&mp_module_builtins) },
extern const struct _mp_obj_module_t mp_module_cmath;
#undef MODULE_DEF_MP_QSTR_CMATH
#define MODULE_DEF_MP_QSTR_CMATH { MP_ROM_QSTR(MP_QSTR_cmath), MP_ROM_PTR(&mp_module_cmath) },
extern const struct _mp_obj_module_t mp_module_math;
#undef MODULE_DEF_MP_QSTR_MATH
#define MODULE_DEF_MP_QSTR_MATH { MP_ROM_QSTR(MP_QSTR_math), MP_ROM_PTR(&mp_module_math) },
extern const struct _mp_obj_module_t mp_module_micropython;
#undef MODULE_DEF_MP_QSTR_MICROPYTHON
#define MODULE_DEF_MP_QSTR_MICROPYTHON { MP_ROM_QSTR(MP_QSTR_micropython), MP_ROM_PTR(&mp_module_micropython) },
extern const struct _mp_obj_module_t mp_module_urandom;
#undef MODULE_DEF_MP_QSTR_URANDOM
#define MODULE_DEF_MP_QSTR_URANDOM { MP_ROM_QSTR(MP_QSTR_urandom), MP_ROM_PTR(&mp_module_urandom) },
// Upsilon's modules part
extern const struct _mp_obj_module_t modion_module;
#undef MODULE_DEF_MP_QSTR_ION
#define MODULE_DEF_MP_QSTR_ION { MP_ROM_QSTR(MP_QSTR_ion), MP_ROM_PTR(&modion_module) },
extern const struct _mp_obj_module_t modkandinsky_module;
#undef MODULE_DEF_MP_QSTR_KANDINSKY
#define MODULE_DEF_MP_QSTR_KANDINSKY { MP_ROM_QSTR(MP_QSTR_kandinsky), MP_ROM_PTR(&modkandinsky_module) },
extern const struct _mp_obj_module_t modmatplotlib_module;
#undef MODULE_DEF_MP_QSTR_MATPLOTLIB
#define MODULE_DEF_MP_QSTR_MATPLOTLIB { MP_ROM_QSTR(MP_QSTR_matplotlib), MP_ROM_PTR(&modmatplotlib_module) },
extern const struct _mp_obj_module_t modpyplot_module;
#undef MODULE_DEF_MP_QSTR_PYPLOT
#define MODULE_DEF_MP_QSTR_PYPLOT { MP_ROM_QSTR(MP_QSTR_matplotlib_dot_pyplot), MP_ROM_PTR(&modpyplot_module) },
extern const struct _mp_obj_module_t modtime_module;
#undef MODULE_DEF_MP_QSTR_TIME
#define MODULE_DEF_MP_QSTR_TIME { MP_ROM_QSTR(MP_QSTR_time), MP_ROM_PTR(&modtime_module) },
extern const struct _mp_obj_module_t modos_module;
#undef MODULE_DEF_MP_QSTR_OS
#define MODULE_DEF_MP_QSTR_OS { MP_ROM_QSTR(MP_QSTR_os), MP_ROM_PTR(&modos_module) },
extern const struct _mp_obj_module_t modturtle_module;
#undef MODULE_DEF_MP_QSTR_TURTLE
#define MODULE_DEF_MP_QSTR_TURTLE { MP_ROM_QSTR(MP_QSTR_turtle), MP_ROM_PTR(&modturtle_module) },
#if !defined(INCLUDE_ULAB)
#define MICROPY_REGISTERED_MODULES \
MODULE_DEF_MP_QSTR_UARRAY \
MODULE_DEF_MP_QSTR_BUILTINS \
MODULE_DEF_MP_QSTR_CMATH \
MODULE_DEF_MP_QSTR_MATH \
MODULE_DEF_MP_QSTR_MICROPYTHON \
MODULE_DEF_MP_QSTR_URANDOM \
MODULE_DEF_MP_QSTR___MAIN__ \
/* Upsilon's modules part */ \
MODULE_DEF_MP_QSTR_ION \
MODULE_DEF_MP_QSTR_KANDINSKY \
MODULE_DEF_MP_QSTR_MATPLOTLIB \
MODULE_DEF_MP_QSTR_PYPLOT \
MODULE_DEF_MP_QSTR_TIME \
MODULE_DEF_MP_QSTR_OS \
MODULE_DEF_MP_QSTR_TURTLE
#else
extern const struct _mp_obj_module_t ulab_user_cmodule;
#undef MODULE_DEF_MP_QSTR_ULAB
#define MODULE_DEF_MP_QSTR_ULAB { MP_ROM_QSTR(MP_QSTR_ulab), MP_ROM_PTR(&ulab_user_cmodule) },
#define MICROPY_REGISTERED_MODULES \
MODULE_DEF_MP_QSTR_BUILTINS \
MODULE_DEF_MP_QSTR_CMATH \
MODULE_DEF_MP_QSTR_MATH \
MODULE_DEF_MP_QSTR_MICROPYTHON \
MODULE_DEF_MP_QSTR_URANDOM \
MODULE_DEF_MP_QSTR___MAIN__ \
/* Upsilon's modules part */ \
MODULE_DEF_MP_QSTR_ION \
MODULE_DEF_MP_QSTR_KANDINSKY \
MODULE_DEF_MP_QSTR_MATPLOTLIB \
MODULE_DEF_MP_QSTR_PYPLOT \
MODULE_DEF_MP_QSTR_TIME \
MODULE_DEF_MP_QSTR_OS \
MODULE_DEF_MP_QSTR_TURTLE \
MODULE_DEF_MP_QSTR_ULAB
#endif
// MICROPY_REGISTERED_MODULES

4
python/port/genhdr/mpversion.h Normal file
View File

@@ -0,0 +1,4 @@
// This file was generated by py/makeversionhdr.py
#define MICROPY_GIT_TAG "v1.19.1"
#define MICROPY_GIT_HASH "9b486340d"
#define MICROPY_BUILD_DATE "2022-06-22"

24
python/port/genhdr/qstrdefs.in.h
View File

@@ -74,6 +74,7 @@ Q(__call__)
Q(__class__)
Q(__contains__)
Q(__delitem__)
Q(__dir__)
Q(__divmod__)
Q(__enter__)
Q(__eq__)
@@ -101,9 +102,7 @@ Q(__mod__)
Q(__module__)
Q(__mul__)
Q(__name__)
#if __EMSCRIPTEN__
Q(__ne__)
#endif
Q(__neg__)
Q(__new__)
Q(__next__)
@@ -142,7 +141,6 @@ Q(all)
Q(any)
Q(append)
Q(args)
Q(argv)
Q(asin)
Q(asinh)
Q(atan)
@@ -154,7 +152,6 @@ Q(bound_method)
Q(builtins)
Q(bytearray)
Q(bytecode)
Q(byteorder)
Q(bytes)
Q(callable)
Q(ceil)
@@ -192,7 +189,6 @@ Q(erfc)
Q(errno)
Q(eval)
Q(exec)
Q(exit)
Q(exp)
Q(expm1)
Q(extend)
@@ -223,14 +219,12 @@ Q(heap_unlock)
Q(hex)
Q(id)
Q(imag)
Q(implementation)
Q(index)
Q(input)
Q(insert)
Q(int)
Q(intersection)
Q(intersection_update)
Q(ion)
Q(isalpha)
Q(isdigit)
Q(isdisjoint)
@@ -248,7 +242,6 @@ Q(items)
Q(iter)
Q(iterator)
Q(join)
Q(kandinsky)
Q(kbd_intr)
Q(key)
Q(keys)
@@ -265,23 +258,18 @@ Q(lower)
Q(lstrip)
Q(map)
Q(math)
Q(matplotlib)
Q(matplotlib.pyplot)
Q(max)
Q(maximum_space_recursion_space_depth_space_exceeded)
Q(micropython)
Q(min)
Q(modf)
Q(module)
Q(modules)
Q(next)
Q(object)
Q(oct)
Q(open)
Q(opt_level)
Q(ord)
Q(os)
Q(path)
Q(pend_throw)
Q(phase)
Q(pi)
@@ -290,7 +278,6 @@ Q(pop)
Q(popitem)
Q(pow)
Q(print)
Q(print_exception)
Q(property)
Q(radians)
Q(randint)
@@ -334,26 +321,20 @@ Q(sum)
Q(super)
Q(symmetric_difference)
Q(symmetric_difference_update)
Q(sys)
Q(tan)
Q(tanh)
Q(throw)
Q(time)
Q(to_bytes)
Q(trunc)
Q(tuple)
Q(turtle)
Q(type)
Q(uniform)
Q(union)
Q(update)
Q(upper)
Q(random)
Q(sys)
Q(value)
Q(values)
Q(version)
Q(version_info)
Q(zip)
// Ion QSTR
@@ -643,7 +624,9 @@ Q(set_printoptions)
Q(get_printoptions)
Q(ndinfo)
Q(arange)
Q(compress)
Q(concatenate)
Q(delete)
Q(diag)
Q(empty)
Q(eye)
@@ -706,6 +689,7 @@ Q(byteswap)
Q(flatten)
Q(k)
Q(tobytes)
Q(tolist)
Q(M)
Q(ulab)
Q(num)

17
python/port/mod/matplotlib/modmatplotlib_table.c
View File

@@ -4,13 +4,26 @@ extern const mp_obj_module_t modpyplot_module;
STATIC MP_DEFINE_CONST_FUN_OBJ_0(modmatplotlib___init___obj, modmatplotlib___init__);
STATIC const mp_rom_map_elem_t modmatplotlib_module_globals_table[] = {
// Define the module table as non-const, because MicroPython needs to be able to modify it.
STATIC mp_rom_map_elem_t modmatplotlib_module_globals_table[] = {
{ MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_matplotlib) },
{ MP_ROM_QSTR(MP_QSTR___init__), MP_ROM_PTR(&modmatplotlib___init___obj) },
{ MP_ROM_QSTR(MP_QSTR_pyplot), MP_ROM_PTR(&modpyplot_module) }
};
STATIC MP_DEFINE_CONST_DICT(modmatplotlib_module_globals, modmatplotlib_module_globals_table);
// Define the module object, not as a constant, because MicroPython needs to be able to dynamically add attributes to it.
mp_obj_dict_t modmatplotlib_module_globals = { \
.base = {&mp_type_dict}, \
.map = { \
.all_keys_are_qstrs = 1, \
.is_fixed = 0, \
.is_ordered = 1, \
.used = MP_ARRAY_SIZE(modmatplotlib_module_globals_table), \
.alloc = MP_ARRAY_SIZE(modmatplotlib_module_globals_table), \
.table = (mp_map_elem_t *)(mp_rom_map_elem_t *)modmatplotlib_module_globals_table, \
}, \
};
const mp_obj_module_t modmatplotlib_module = {
.base = { &mp_type_module },

18
python/port/mod/ulab/micropython.cmake Normal file
View File

@@ -0,0 +1,18 @@
add_library(usermod_ulab INTERFACE)
file(GLOB_RECURSE ULAB_SOURCES ${CMAKE_CURRENT_LIST_DIR}/*.c)
target_sources(usermod_ulab INTERFACE
${ULAB_SOURCES}
)
target_include_directories(usermod_ulab INTERFACE
${CMAKE_CURRENT_LIST_DIR}
)
target_compile_definitions(usermod_ulab INTERFACE
MODULE_ULAB_ENABLED=1
)
target_link_libraries(usermod INTERFACE usermod_ulab)

39
python/port/mod/ulab/micropython.mk Normal file
View File

@@ -0,0 +1,39 @@
USERMODULES_DIR := $(USERMOD_DIR)
# Add all C files to SRC_USERMOD.
SRC_USERMOD += $(USERMODULES_DIR)/scipy/linalg/linalg.c
SRC_USERMOD += $(USERMODULES_DIR)/scipy/optimize/optimize.c
SRC_USERMOD += $(USERMODULES_DIR)/scipy/signal/signal.c
SRC_USERMOD += $(USERMODULES_DIR)/scipy/special/special.c
SRC_USERMOD += $(USERMODULES_DIR)/ndarray_operators.c
SRC_USERMOD += $(USERMODULES_DIR)/ulab_tools.c
SRC_USERMOD += $(USERMODULES_DIR)/ndarray.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/ndarray/ndarray_iter.c
SRC_USERMOD += $(USERMODULES_DIR)/ndarray_properties.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/approx.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/compare.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/carray/carray.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/carray/carray_tools.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/create.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/fft/fft.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/fft/fft_tools.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/filter.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/io/io.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/linalg/linalg.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/linalg/linalg_tools.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/numerical.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/poly.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/stats.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/transform.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/vector.c
SRC_USERMOD += $(USERMODULES_DIR)/numpy/numpy.c
SRC_USERMOD += $(USERMODULES_DIR)/scipy/scipy.c
SRC_USERMOD += $(USERMODULES_DIR)/user/user.c
SRC_USERMOD += $(USERMODULES_DIR)/utils/utils.c
SRC_USERMOD += $(USERMODULES_DIR)/ulab.c
CFLAGS_USERMOD += -I$(USERMODULES_DIR)
override CFLAGS_EXTRA += -DMODULE_ULAB_ENABLED=1

632
python/port/mod/ulab/ndarray.c
View File

File diff suppressed because it is too large Load Diff

121
python/port/mod/ulab/ndarray.h
View File

@@ -63,6 +63,8 @@ typedef struct _mp_obj_slice_t {
void ndarray_set_value(char , void *, size_t , mp_obj_t );
#endif
void ndarray_set_complex_value(void *, size_t , mp_obj_t );
#define NDARRAY_NUMERIC 0
#define NDARRAY_BOOLEAN 1
@@ -77,6 +79,9 @@ enum NDARRAY_TYPE {
NDARRAY_INT8 = 'b',
NDARRAY_UINT16 = 'H',
NDARRAY_INT16 = 'h',
#if ULAB_SUPPORTS_COMPLEX
NDARRAY_COMPLEX = 'c',
#endif
NDARRAY_FLOAT = FLOAT_TYPECODE,
};
@@ -131,6 +136,7 @@ void ndarray_assign_elements(ndarray_obj_t *, mp_obj_t , uint8_t , size_t *);
size_t *ndarray_contract_shape(ndarray_obj_t *, uint8_t );
int32_t *ndarray_contract_strides(ndarray_obj_t *, uint8_t );
ndarray_obj_t *ndarray_from_iterable(mp_obj_t , uint8_t );
ndarray_obj_t *ndarray_new_dense_ndarray(uint8_t , size_t *, uint8_t );
ndarray_obj_t *ndarray_new_ndarray_from_tuple(mp_obj_tuple_t *, uint8_t );
ndarray_obj_t *ndarray_new_ndarray(uint8_t , size_t *, int32_t *, uint8_t );
@@ -138,7 +144,8 @@ ndarray_obj_t *ndarray_new_linear_array(size_t , uint8_t );
ndarray_obj_t *ndarray_new_view(ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t );
bool ndarray_is_dense(ndarray_obj_t *);
ndarray_obj_t *ndarray_copy_view(ndarray_obj_t *);
void ndarray_copy_array(ndarray_obj_t *, ndarray_obj_t *);
ndarray_obj_t *ndarray_copy_view_convert_type(ndarray_obj_t *, uint8_t );
void ndarray_copy_array(ndarray_obj_t *, ndarray_obj_t *, uint8_t );
MP_DECLARE_CONST_FUN_OBJ_KW(ndarray_array_constructor_obj);
mp_obj_t ndarray_make_new(const mp_obj_type_t *, size_t , size_t , const mp_obj_t *);
@@ -185,6 +192,11 @@ mp_obj_t ndarray_tobytes(mp_obj_t );
MP_DECLARE_CONST_FUN_OBJ_1(ndarray_tobytes_obj);
#endif
#if NDARRAY_HAS_TOBYTES
mp_obj_t ndarray_tolist(mp_obj_t );
MP_DECLARE_CONST_FUN_OBJ_1(ndarray_tolist_obj);
#endif
#if NDARRAY_HAS_TRANSPOSE
mp_obj_t ndarray_transpose(mp_obj_t );
MP_DECLARE_CONST_FUN_OBJ_1(ndarray_transpose_obj);
@@ -201,15 +213,15 @@ mp_int_t ndarray_get_buffer(mp_obj_t , mp_buffer_info_t *, mp_uint_t );
ndarray_obj_t *ndarray_from_mp_obj(mp_obj_t , uint8_t );
#define BOOLEAN_ASSIGNMENT_LOOP(type_left, type_right, ndarray, iarray, istride, varray, vstride)\
#define BOOLEAN_ASSIGNMENT_LOOP(type_left, type_right, ndarray, lstrides, iarray, istride, varray, vstride)\
type_left *array = (type_left *)(ndarray)->array;\
for(size_t i=0; i < (ndarray)->len; i++) {\
if(*(iarray)) {\
*array = (type_left)(*((type_right *)(varray)));\
(varray) += (vstride);\
}\
array += (ndarray)->strides[ULAB_MAX_DIMS - 1] / (ndarray)->itemsize;\
array += (lstrides);\
(iarray) += (istride);\
(varray) += (vstride);\
} while(0)
#if ULAB_HAS_FUNCTION_ITERATOR
@@ -634,105 +646,4 @@ ndarray_obj_t *ndarray_from_mp_obj(mp_obj_t , uint8_t );
#endif /* ULAB_MAX_DIMS == 4 */
#endif /* ULAB_HAS_FUNCTION_ITERATOR */
#if ULAB_MAX_DIMS == 1
#define ASSIGNMENT_LOOP(results, type_left, type_right, lstrides, rarray, rstrides)\
type_left *larray = (type_left *)(results)->array;\
size_t l = 0;\
do {\
*larray = (type_left)(*((type_right *)(rarray)));\
(larray) += (lstrides)[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 1]);\
#endif /* ULAB_MAX_DIMS == 1 */
#if ULAB_MAX_DIMS == 2
#define ASSIGNMENT_LOOP(results, type_left, type_right, lstrides, rarray, rstrides)\
type_left *larray = (type_left *)(results)->array;\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
*larray = (type_left)(*((type_right *)(rarray)));\
(larray) += (lstrides)[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 1]);\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS-1];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS-1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 2]);\
#endif /* ULAB_MAX_DIMS == 2 */
#if ULAB_MAX_DIMS == 3
#define ASSIGNMENT_LOOP(results, type_left, type_right, lstrides, rarray, rstrides)\
type_left *larray = (type_left *)(results)->array;\
size_t j = 0;\
do {\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
*larray = (type_left)(*((type_right *)(rarray)));\
(larray) += (lstrides)[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 1]);\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS-1];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS-1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 2]);\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 3];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 3];\
j++;\
} while(j < (results)->shape[ULAB_MAX_DIMS - 3]);\
#endif /* ULAB_MAX_DIMS == 3 */
#if ULAB_MAX_DIMS == 4
#define ASSIGNMENT_LOOP(results, type_left, type_right, lstrides, rarray, rstrides)\
type_left *larray = (type_left *)(results)->array;\
size_t i = 0;\
do {\
size_t j = 0;\
do {\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
*larray = (type_left)(*((type_right *)(rarray)));\
(larray) += (lstrides)[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 1]);\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS-1];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS-1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 2]);\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 3];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 3];\
j++;\
} while(j < (results)->shape[ULAB_MAX_DIMS - 3]);\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS-3];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 4];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS-3];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 4];\
i++;\
} while(i < (results)->shape[ULAB_MAX_DIMS - 4]);\
#endif /* ULAB_MAX_DIMS == 4 */
#endif

34
python/port/mod/ulab/ndarray_operators.c
View File

@@ -17,6 +17,7 @@
#include "ndarray_operators.h"
#include "ulab.h"
#include "ulab_tools.h"
#include "numpy/carray/carray.h"
/*
This file contains the actual implementations of the various
@@ -24,7 +25,8 @@
These are the upcasting rules of the binary operators
- if one of the operarands is a float, the result is always float
- if complex is supported, and if one of the operarands is a complex, the result is always complex
- if both operarands are real one of them is a float, then the result is also a float
- operation on identical types preserves type
uint8 + int8 => int16
@@ -39,6 +41,12 @@
mp_obj_t ndarray_binary_equality(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides, mp_binary_op_t op) {
#if ULAB_SUPPORTS_COMPLEX
if((lhs->dtype == NDARRAY_COMPLEX) || (rhs->dtype == NDARRAY_COMPLEX)) {
return carray_binary_equal_not_equal(lhs, rhs, ndim, shape, lstrides, rstrides, op);
}
#endif
ndarray_obj_t *results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_UINT8);
results->boolean = 1;
uint8_t *array = (uint8_t *)results->array;
@@ -161,6 +169,12 @@ mp_obj_t ndarray_binary_equality(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
mp_obj_t ndarray_binary_add(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {
#if ULAB_SUPPORTS_COMPLEX
if((lhs->dtype == NDARRAY_COMPLEX) || (rhs->dtype == NDARRAY_COMPLEX)) {
return carray_binary_add(lhs, rhs, ndim, shape, lstrides, rstrides);
}
#endif
ndarray_obj_t *results = NULL;
uint8_t *larray = (uint8_t *)lhs->array;
uint8_t *rarray = (uint8_t *)rhs->array;
@@ -238,6 +252,12 @@ mp_obj_t ndarray_binary_add(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
mp_obj_t ndarray_binary_multiply(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {
#if ULAB_SUPPORTS_COMPLEX
if((lhs->dtype == NDARRAY_COMPLEX) || (rhs->dtype == NDARRAY_COMPLEX)) {
return carray_binary_multiply(lhs, rhs, ndim, shape, lstrides, rstrides);
}
#endif
ndarray_obj_t *results = NULL;
uint8_t *larray = (uint8_t *)lhs->array;
uint8_t *rarray = (uint8_t *)rhs->array;
@@ -460,6 +480,12 @@ mp_obj_t ndarray_binary_more(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
mp_obj_t ndarray_binary_subtract(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {
#if ULAB_SUPPORTS_COMPLEX
if((lhs->dtype == NDARRAY_COMPLEX) || (rhs->dtype == NDARRAY_COMPLEX)) {
return carray_binary_subtract(lhs, rhs, ndim, shape, lstrides, rstrides);
}
#endif
ndarray_obj_t *results = NULL;
uint8_t *larray = (uint8_t *)lhs->array;
uint8_t *rarray = (uint8_t *)rhs->array;
@@ -559,6 +585,12 @@ mp_obj_t ndarray_binary_subtract(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
mp_obj_t ndarray_binary_true_divide(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {
#if ULAB_SUPPORTS_COMPLEX
if((lhs->dtype == NDARRAY_COMPLEX) || (rhs->dtype == NDARRAY_COMPLEX)) {
return carray_binary_divide(lhs, rhs, ndim, shape, lstrides, rstrides);
}
#endif
ndarray_obj_t *results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
uint8_t *larray = (uint8_t *)lhs->array;
uint8_t *rarray = (uint8_t *)rhs->array;

15
python/port/mod/ulab/ndarray_properties.c
View File

@@ -20,6 +20,9 @@
#include "ulab.h"
#include "ndarray.h"
#include "numpy/ndarray/ndarray_iter.h"
#if ULAB_SUPPORTS_COMPLEX
#include "numpy/carray/carray.h"
#endif
#ifndef CIRCUITPY
@@ -82,6 +85,18 @@ void ndarray_properties_attr(mp_obj_t self_in, qstr attr, mp_obj_t *dest) {
dest[0] = ndarray_transpose(self_in);
break;
#endif
#if ULAB_SUPPORTS_COMPLEX
#if ULAB_NUMPY_HAS_IMAG
case MP_QSTR_imag:
dest[0] = carray_imag(self_in);
break;
#endif
#if ULAB_NUMPY_HAS_IMAG
case MP_QSTR_real:
dest[0] = carray_real(self_in);
break;
#endif
#endif /* ULAB_SUPPORTS_COMPLEX */
default:
call_local_method(self_in, attr, dest);
break;

6
python/port/mod/ulab/numpy/approx.c
View File

@@ -19,6 +19,7 @@
#include "../ulab.h"
#include "../ulab_tools.h"
#include "carray/carray_tools.h"
#include "approx.h"
//| """Numerical approximation methods"""
@@ -60,6 +61,9 @@ STATIC mp_obj_t approx_interp(size_t n_args, const mp_obj_t *pos_args, mp_map_t
ndarray_obj_t *x = ndarray_from_mp_obj(args[0].u_obj, 0);
ndarray_obj_t *xp = ndarray_from_mp_obj(args[1].u_obj, 0); // xp must hold an increasing sequence of independent values
ndarray_obj_t *fp = ndarray_from_mp_obj(args[2].u_obj, 0);
COMPLEX_DTYPE_NOT_IMPLEMENTED(x->dtype)
COMPLEX_DTYPE_NOT_IMPLEMENTED(xp->dtype)
COMPLEX_DTYPE_NOT_IMPLEMENTED(fp->dtype)
if((xp->ndim != 1) || (fp->ndim != 1) || (xp->len < 2) || (fp->len < 2) || (xp->len != fp->len)) {
mp_raise_ValueError(translate("interp is defined for 1D iterables of equal length"));
}
@@ -157,6 +161,7 @@ STATIC mp_obj_t approx_trapz(size_t n_args, const mp_obj_t *pos_args, mp_map_t *
mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
ndarray_obj_t *y = ndarray_from_mp_obj(args[0].u_obj, 0);
COMPLEX_DTYPE_NOT_IMPLEMENTED(y->dtype)
ndarray_obj_t *x;
mp_float_t mean = MICROPY_FLOAT_CONST(0.0);
if(y->len < 2) {
@@ -174,6 +179,7 @@ STATIC mp_obj_t approx_trapz(size_t n_args, const mp_obj_t *pos_args, mp_map_t *
if(args[1].u_obj != mp_const_none) {
x = ndarray_from_mp_obj(args[1].u_obj, 0); // x must hold an increasing sequence of independent values
COMPLEX_DTYPE_NOT_IMPLEMENTED(x->dtype)
if((x->ndim != 1) || (y->len != x->len)) {
mp_raise_ValueError(translate("trapz is defined for 1D arrays of equal length"));
}

826
python/port/mod/ulab/numpy/carray/carray.c Normal file
View File

@@ -0,0 +1,826 @@
/*
* This file is part of the micropython-ulab project,
*
* https://github.com/v923z/micropython-ulab
*
* The MIT License (MIT)
*
* Copyright (c) 2021-2022 Zoltán Vörös
*/
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "py/obj.h"
#include "py/objint.h"
#include "py/runtime.h"
#include "py/builtin.h"
#include "py/misc.h"
#include "../../ulab.h"
#include "../../ndarray.h"
#include "../../ulab_tools.h"
#include "carray.h"
#if ULAB_SUPPORTS_COMPLEX
//| import ulab.numpy
//| def real(val):
//| """
//| Return the real part of the complex argument, which can be
//| either an ndarray, or a scalar."""
//| ...
//|
mp_obj_t carray_real(mp_obj_t _source) {
if(mp_obj_is_type(_source, &ulab_ndarray_type)) {
ndarray_obj_t *source = MP_OBJ_TO_PTR(_source);
if(source->dtype != NDARRAY_COMPLEX) {
ndarray_obj_t *target = ndarray_new_dense_ndarray(source->ndim, source->shape, source->dtype);
ndarray_copy_array(source, target, 0);
return MP_OBJ_FROM_PTR(target);
} else { // the input is most definitely a complex array
ndarray_obj_t *target = ndarray_new_dense_ndarray(source->ndim, source->shape, NDARRAY_FLOAT);
ndarray_copy_array(source, target, 0);
return MP_OBJ_FROM_PTR(target);
}
} else {
mp_raise_NotImplementedError(translate("function is implemented for ndarrays only"));
}
return mp_const_none;
}
MP_DEFINE_CONST_FUN_OBJ_1(carray_real_obj, carray_real);
//| def imag(val):
//| """
//| Return the imaginary part of the complex argument, which can be
//| either an ndarray, or a scalar."""
//| ...
//|
mp_obj_t carray_imag(mp_obj_t _source) {
if(mp_obj_is_type(_source, &ulab_ndarray_type)) {
ndarray_obj_t *source = MP_OBJ_TO_PTR(_source);
if(source->dtype != NDARRAY_COMPLEX) { // if not complex, then the imaginary part is zero
ndarray_obj_t *target = ndarray_new_dense_ndarray(source->ndim, source->shape, source->dtype);
return MP_OBJ_FROM_PTR(target);
} else { // the input is most definitely a complex array
ndarray_obj_t *target = ndarray_new_dense_ndarray(source->ndim, source->shape, NDARRAY_FLOAT);
ndarray_copy_array(source, target, source->itemsize / 2);
return MP_OBJ_FROM_PTR(target);
}
} else {
mp_raise_NotImplementedError(translate("function is implemented for ndarrays only"));
}
return mp_const_none;
}
MP_DEFINE_CONST_FUN_OBJ_1(carray_imag_obj, carray_imag);
#if ULAB_NUMPY_HAS_CONJUGATE
//| def conjugate(val):
//| """
//| Return the conjugate of the complex argument, which can be
//| either an ndarray, or a scalar."""
//| ...
//|
mp_obj_t carray_conjugate(mp_obj_t _source) {
if(mp_obj_is_type(_source, &ulab_ndarray_type)) {
ndarray_obj_t *source = MP_OBJ_TO_PTR(_source);
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(source->ndim, source->shape, source->dtype);
ndarray_copy_array(source, ndarray, 0);
if(source->dtype == NDARRAY_COMPLEX) {
mp_float_t *array = (mp_float_t *)ndarray->array;
array++;
for(size_t i = 0; i < ndarray->len; i++) {
*array *= MICROPY_FLOAT_CONST(-1.0);
array += 2;
}
}
return MP_OBJ_FROM_PTR(ndarray);
} else {
if(mp_obj_is_type(_source, &mp_type_complex)) {
mp_float_t real, imag;
mp_obj_get_complex(_source, &real, &imag);
imag = imag * MICROPY_FLOAT_CONST(-1.0);
return mp_obj_new_complex(real, imag);
} else if(mp_obj_is_int(_source) || mp_obj_is_float(_source)) {
return _source;
} else {
mp_raise_TypeError(translate("input must be an ndarray, or a scalar"));
}
}
// this should never happen
return mp_const_none;
}
MP_DEFINE_CONST_FUN_OBJ_1(carray_conjugate_obj, carray_conjugate);
#endif
#if ULAB_NUMPY_HAS_SORT_COMPLEX
//| def sort_complex(a: ulab.numpy.ndarray) -> ulab.numpy.ndarray:
//| """
//| .. param: a
//| a one-dimensional ndarray
//|
//| Sort a complex array using the real part first, then the imaginary part.
//| Always returns a sorted complex array, even if the input was real."""
//| ...
//|
static void carray_sort_complex_(mp_float_t *array, size_t len) {
// array is assumed to be a floating vector containing the real and imaginary parts
// of a complex array at alternating positions as
// array[0] = real[0]
// array[1] = imag[0]
// array[2] = real[1]
// array[3] = imag[1]
mp_float_t real, imag;
size_t c, q = len, p, r = len >> 1;
for (;;) {
if (r > 0) {
r--;
real = array[2 * r];
imag = array[2 * r + 1];
} else {
q--;
if(q == 0) {
break;
}
real = array[2 * q];
imag = array[2 * q + 1];
array[2 * q] = array[0];
array[2 * q + 1] = array[1];
}
p = r;
c = r + r + 1;
while (c < q) {
if(c + 1 < q) {
if((array[2 * (c+1)] > array[2 * c]) ||
((array[2 * (c+1)] == array[2 * c]) && (array[2 * (c+1) + 1] > array[2 * c + 1]))) {
c++;
}
}
if((array[2 * c] > real) ||
((array[2 * c] == real) && (array[2 * c + 1] > imag))) {
array[2 * p] = array[2 * c]; // real part
array[2 * p + 1] = array[2 * c + 1]; // imag part
p = c;
c = p + p + 1;
} else {
break;
}
}
array[2 * p] = real;
array[2 * p + 1] = imag;
}
}
mp_obj_t carray_sort_complex(mp_obj_t _source) {
if(!mp_obj_is_type(_source, &ulab_ndarray_type)) {
mp_raise_TypeError(translate("input must be a 1D ndarray"));
}
ndarray_obj_t *source = MP_OBJ_TO_PTR(_source);
if(source->ndim != 1) {
mp_raise_TypeError(translate("input must be a 1D ndarray"));
}
ndarray_obj_t *ndarray = ndarray_copy_view_convert_type(source, NDARRAY_COMPLEX);
mp_float_t *array = (mp_float_t *)ndarray->array;
carray_sort_complex_(array, ndarray->len);
return MP_OBJ_FROM_PTR(ndarray);
}
MP_DEFINE_CONST_FUN_OBJ_1(carray_sort_complex_obj, carray_sort_complex);
#endif
//| def abs(a: ulab.numpy.ndarray) -> ulab.numpy.ndarray:
//| """
//| .. param: a
//| a one-dimensional ndarray
//|
//| Return the absolute value of complex ndarray."""
//| ...
//|
mp_obj_t carray_abs(ndarray_obj_t *source, ndarray_obj_t *target) {
// calculates the absolute value of a complex array and returns a dense array
uint8_t *sarray = (uint8_t *)source->array;
mp_float_t *tarray = (mp_float_t *)target->array;
uint8_t itemsize = mp_binary_get_size('@', NDARRAY_FLOAT, NULL);
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
mp_float_t rvalue = *(mp_float_t *)sarray;
mp_float_t ivalue = *(mp_float_t *)(sarray + itemsize);
*tarray++ = MICROPY_FLOAT_C_FUN(sqrt)(rvalue * rvalue + ivalue * ivalue);
sarray += source->strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < source->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
sarray -= source->strides[ULAB_MAX_DIMS - 1] * source->shape[ULAB_MAX_DIMS-1];
sarray += source->strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < source->shape[ULAB_MAX_DIMS - 2]);
#endif
#if ULAB_MAX_DIMS > 2
sarray -= source->strides[ULAB_MAX_DIMS - 2] * source->shape[ULAB_MAX_DIMS-2];
sarray += source->strides[ULAB_MAX_DIMS - 3];
j++;
} while(j < source->shape[ULAB_MAX_DIMS - 3]);
#endif
#if ULAB_MAX_DIMS > 3
sarray -= source->strides[ULAB_MAX_DIMS - 3] * source->shape[ULAB_MAX_DIMS-3];
sarray += source->strides[ULAB_MAX_DIMS - 4];
i++;
} while(i < source->shape[ULAB_MAX_DIMS - 4]);
#endif
return MP_OBJ_FROM_PTR(target);
}
static void carray_copy_part(uint8_t *tarray, uint8_t *sarray, size_t *shape, int32_t *strides) {
// copies the real or imaginary part of an array
// into the respective part of a dense complex array
uint8_t sz = sizeof(mp_float_t);
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
memcpy(tarray, sarray, sz);
tarray += 2 * sz;
sarray += strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
sarray -= strides[ULAB_MAX_DIMS - 1] * shape[ULAB_MAX_DIMS-1];
sarray += strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_MAX_DIMS > 2
sarray -= strides[ULAB_MAX_DIMS - 2] * shape[ULAB_MAX_DIMS-2];
sarray += strides[ULAB_MAX_DIMS - 3];
j++;
} while(j < shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS > 2 */
#if ULAB_MAX_DIMS > 3
sarray -= strides[ULAB_MAX_DIMS - 3] * shape[ULAB_MAX_DIMS-3];
sarray += strides[ULAB_MAX_DIMS - 4];
i++;
} while(i < shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS > 3 */
}
mp_obj_t carray_binary_equal_not_equal(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides, mp_binary_op_t op) {
ndarray_obj_t *results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_UINT8);
results->boolean = 1;
uint8_t *array = (uint8_t *)results->array;
if(op == MP_BINARY_OP_NOT_EQUAL) {
memset(array, 1, results->len);
}
if((lhs->dtype == NDARRAY_COMPLEX) && (rhs->dtype == NDARRAY_COMPLEX)) {
mp_float_t *larray = (mp_float_t *)lhs->array;
mp_float_t *rarray = (mp_float_t *)rhs->array;
ulab_rescale_float_strides(lstrides);
ulab_rescale_float_strides(rstrides);
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
if((larray[0] == rarray[0]) && (larray[1] == rarray[1])) {
*array ^= 0x01;
}
array++;
larray += lstrides[ULAB_MAX_DIMS - 1];
rarray += rstrides[ULAB_MAX_DIMS - 1];
l++;
} while(l < results->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
larray -= lstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
larray += lstrides[ULAB_MAX_DIMS - 2];
rarray -= rstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
rarray += rstrides[ULAB_MAX_DIMS - 2];
k++;
} while(k < results->shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_MAX_DIMS > 2
larray -= lstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
larray += lstrides[ULAB_MAX_DIMS - 3];
rarray -= rstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
rarray += rstrides[ULAB_MAX_DIMS - 3];
j++;
} while(j < results->shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS > 2 */
#if ULAB_MAX_DIMS > 3
larray -= lstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
larray += lstrides[ULAB_MAX_DIMS - 4];
rarray -= rstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
rarray += rstrides[ULAB_MAX_DIMS - 4];
i++;
} while(i < results->shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS > 3 */
} else { // only one of the operands is complex
mp_float_t *larray = (mp_float_t *)lhs->array;
uint8_t *rarray = (uint8_t *)rhs->array;
// align the complex array to the left
uint8_t rdtype = rhs->dtype;
int32_t *lstrides_ = lstrides;
int32_t *rstrides_ = rstrides;
if(rhs->dtype == NDARRAY_COMPLEX) {
larray = (mp_float_t *)rhs->array;
rarray = (uint8_t *)lhs->array;
lstrides_ = rstrides;
rstrides_ = lstrides;
rdtype = lhs->dtype;
}
ulab_rescale_float_strides(lstrides_);
if(rdtype == NDARRAY_UINT8) {
BINARY_LOOP_COMPLEX_EQUAL(results, array, uint8_t, larray, lstrides_, rarray, rstrides_);
} else if(rdtype == NDARRAY_INT8) {
BINARY_LOOP_COMPLEX_EQUAL(results, array, int8_t, larray, lstrides_, rarray, rstrides_);
} else if(rdtype == NDARRAY_UINT16) {
BINARY_LOOP_COMPLEX_EQUAL(results, array, uint16_t, larray, lstrides_, rarray, rstrides_);
} else if(rdtype == NDARRAY_INT16) {
BINARY_LOOP_COMPLEX_EQUAL(results, array, int16_t, larray, lstrides_, rarray, rstrides_);
} else if(rdtype == NDARRAY_FLOAT) {
BINARY_LOOP_COMPLEX_EQUAL(results, array, mp_float_t, larray, lstrides_, rarray, rstrides_);
}
}
return MP_OBJ_FROM_PTR(results);
}
mp_obj_t carray_binary_add(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {
ndarray_obj_t *results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_COMPLEX);
mp_float_t *resarray = (mp_float_t *)results->array;
if((lhs->dtype == NDARRAY_COMPLEX) && (rhs->dtype == NDARRAY_COMPLEX)) {
mp_float_t *larray = (mp_float_t *)lhs->array;
mp_float_t *rarray = (mp_float_t *)rhs->array;
ulab_rescale_float_strides(lstrides);
ulab_rescale_float_strides(rstrides);
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
// real part
*resarray++ = larray[0] + rarray[0];
// imaginary part
*resarray++ = larray[1] + rarray[1];
larray += lstrides[ULAB_MAX_DIMS - 1];
rarray += rstrides[ULAB_MAX_DIMS - 1];
l++;
} while(l < results->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
larray -= lstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
larray += lstrides[ULAB_MAX_DIMS - 2];
rarray -= rstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
rarray += rstrides[ULAB_MAX_DIMS - 2];
k++;
} while(k < results->shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_MAX_DIMS > 2
larray -= lstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
larray += lstrides[ULAB_MAX_DIMS - 3];
rarray -= rstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
rarray += rstrides[ULAB_MAX_DIMS - 3];
j++;
} while(j < results->shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS > 2 */
#if ULAB_MAX_DIMS > 3
larray -= lstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
larray += lstrides[ULAB_MAX_DIMS - 4];
rarray -= rstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
rarray += rstrides[ULAB_MAX_DIMS - 4];
i++;
} while(i < results->shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS > 3 */
} else { // only one of the operands is complex
uint8_t *larray = (uint8_t *)lhs->array;
uint8_t *rarray = (uint8_t *)rhs->array;
// align the complex array to the left
uint8_t rdtype = rhs->dtype;
int32_t *lstrides_ = lstrides;
int32_t *rstrides_ = rstrides;
if(rhs->dtype == NDARRAY_COMPLEX) {
larray = (uint8_t *)rhs->array;
rarray = (uint8_t *)lhs->array;
lstrides_ = rstrides;
rstrides_ = lstrides;
rdtype = lhs->dtype;
}
if(rdtype == NDARRAY_UINT8) {
BINARY_LOOP_COMPLEX(results, resarray, uint8_t, larray, lstrides_, rarray, rstrides_, +);
} else if(rdtype == NDARRAY_INT8) {
BINARY_LOOP_COMPLEX(results, resarray, int8_t, larray, lstrides_, rarray, rstrides_, +);
} else if(rdtype == NDARRAY_UINT16) {
BINARY_LOOP_COMPLEX(results, resarray, uint16_t, larray, lstrides_, rarray, rstrides_, +);
} else if(rdtype == NDARRAY_INT16) {
BINARY_LOOP_COMPLEX(results, resarray, int16_t, larray, lstrides_, rarray, rstrides_, +);
} else if(rdtype == NDARRAY_FLOAT) {
BINARY_LOOP_COMPLEX(results, resarray, mp_float_t, larray, lstrides_, rarray, rstrides_, +);
}
// simply copy the imaginary part
uint8_t *tarray = (uint8_t *)results->array;
tarray += sizeof(mp_float_t);
if(lhs->dtype == NDARRAY_COMPLEX) {
rarray = (uint8_t *)lhs->array;
rstrides = lstrides;
} else {
rarray = (uint8_t *)rhs->array;
}
rarray += sizeof(mp_float_t);
carray_copy_part(tarray, rarray, results->shape, rstrides);
}
return MP_OBJ_FROM_PTR(results);
}
static void carray_binary_multiply_(ndarray_obj_t *results, mp_float_t *resarray, uint8_t *larray, uint8_t *rarray,
int32_t *lstrides, int32_t *rstrides, uint8_t rdtype) {
if(rdtype == NDARRAY_UINT8) {
BINARY_LOOP_COMPLEX(results, resarray, uint8_t, larray, lstrides, rarray, rstrides, *);
} else if(rdtype == NDARRAY_INT8) {
BINARY_LOOP_COMPLEX(results, resarray, int8_t, larray, lstrides, rarray, rstrides, *);
} else if(rdtype == NDARRAY_UINT16) {
BINARY_LOOP_COMPLEX(results, resarray, uint16_t, larray, lstrides, rarray, rstrides, *);
} else if(rdtype == NDARRAY_INT16) {
BINARY_LOOP_COMPLEX(results, resarray, int16_t, larray, lstrides, rarray, rstrides, *);
} else if(rdtype == NDARRAY_FLOAT) {
BINARY_LOOP_COMPLEX(results, resarray, mp_float_t, larray, lstrides, rarray, rstrides, *);
}
}
mp_obj_t carray_binary_multiply(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {
ndarray_obj_t *results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_COMPLEX);
mp_float_t *resarray = (mp_float_t *)results->array;
if((lhs->dtype == NDARRAY_COMPLEX) && (rhs->dtype == NDARRAY_COMPLEX)) {
mp_float_t *larray = (mp_float_t *)lhs->array;
mp_float_t *rarray = (mp_float_t *)rhs->array;
ulab_rescale_float_strides(lstrides);
ulab_rescale_float_strides(rstrides);
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
// real part
*resarray++ = larray[0] * rarray[0] - larray[1] * rarray[1];
// imaginary part
*resarray++ = larray[0] * rarray[1] + larray[1] * rarray[0];
larray += lstrides[ULAB_MAX_DIMS - 1];
rarray += rstrides[ULAB_MAX_DIMS - 1];
l++;
} while(l < results->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
larray -= lstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
larray += lstrides[ULAB_MAX_DIMS - 2];
rarray -= rstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
rarray += rstrides[ULAB_MAX_DIMS - 2];
k++;
} while(k < results->shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_MAX_DIMS > 2
larray -= lstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
larray += lstrides[ULAB_MAX_DIMS - 3];
rarray -= rstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
rarray += rstrides[ULAB_MAX_DIMS - 3];
j++;
} while(j < results->shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS > 2 */
#if ULAB_MAX_DIMS > 3
larray -= lstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
larray += lstrides[ULAB_MAX_DIMS - 4];
rarray -= rstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
rarray += rstrides[ULAB_MAX_DIMS - 4];
i++;
} while(i < results->shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS > 3 */
} else { // only one of the operands is complex
uint8_t *larray = (uint8_t *)lhs->array;
uint8_t *rarray = (uint8_t *)rhs->array;
uint8_t *lo = larray, *ro = rarray;
int32_t *left_strides = lstrides;
int32_t *right_strides = rstrides;
uint8_t rdtype = rhs->dtype;
// align the complex array to the left
if(rhs->dtype == NDARRAY_COMPLEX) {
lo = (uint8_t *)rhs->array;
ro = (uint8_t *)lhs->array;
rdtype = lhs->dtype;
left_strides = rstrides;
right_strides = lstrides;
}
larray = lo;
rarray = ro;
// real part
carray_binary_multiply_(results, resarray, larray, rarray, left_strides, right_strides, rdtype);
larray = lo + sizeof(mp_float_t);
rarray = ro;
resarray = (mp_float_t *)results->array;
resarray++;
// imaginary part
carray_binary_multiply_(results, resarray, larray, rarray, left_strides, right_strides, rdtype);
}
return MP_OBJ_FROM_PTR(results);
}
mp_obj_t carray_binary_subtract(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {
ndarray_obj_t *results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_COMPLEX);
mp_float_t *resarray = (mp_float_t *)results->array;
if((lhs->dtype == NDARRAY_COMPLEX) && (rhs->dtype == NDARRAY_COMPLEX)) {
mp_float_t *larray = (mp_float_t *)lhs->array;
mp_float_t *rarray = (mp_float_t *)rhs->array;
ulab_rescale_float_strides(lstrides);
ulab_rescale_float_strides(rstrides);
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
// real part
*resarray++ = larray[0] - rarray[0];
// imaginary part
*resarray++ = larray[1] - rarray[1];
larray += lstrides[ULAB_MAX_DIMS - 1];
rarray += rstrides[ULAB_MAX_DIMS - 1];
l++;
} while(l < results->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
larray -= lstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
larray += lstrides[ULAB_MAX_DIMS - 2];
rarray -= rstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
rarray += rstrides[ULAB_MAX_DIMS - 2];
k++;
} while(k < results->shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_MAX_DIMS > 2
larray -= lstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
larray += lstrides[ULAB_MAX_DIMS - 3];
rarray -= rstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
rarray += rstrides[ULAB_MAX_DIMS - 3];
j++;
} while(j < results->shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS > 2 */
#if ULAB_MAX_DIMS > 3
larray -= lstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
larray += lstrides[ULAB_MAX_DIMS - 4];
rarray -= rstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
rarray += rstrides[ULAB_MAX_DIMS - 4];
i++;
} while(i < results->shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS > 3 */
} else {
uint8_t *larray = (uint8_t *)lhs->array;
if(lhs->dtype == NDARRAY_COMPLEX) {
uint8_t *rarray = (uint8_t *)rhs->array;
if(rhs->dtype == NDARRAY_UINT8) {
BINARY_LOOP_COMPLEX(results, resarray, uint8_t, larray, lstrides, rarray, rstrides, -);
} else if(rhs->dtype == NDARRAY_INT8) {
BINARY_LOOP_COMPLEX(results, resarray, int8_t, larray, lstrides, rarray, rstrides, -);
} else if(rhs->dtype == NDARRAY_UINT16) {
BINARY_LOOP_COMPLEX(results, resarray, uint16_t, larray, lstrides, rarray, rstrides, -);
} else if(rhs->dtype == NDARRAY_INT16) {
BINARY_LOOP_COMPLEX(results, resarray, int16_t, larray, lstrides, rarray, rstrides, -);
} else if(rhs->dtype == NDARRAY_FLOAT) {
BINARY_LOOP_COMPLEX(results, resarray, mp_float_t, larray, lstrides, rarray, rstrides, -);
}
// copy the imaginary part
uint8_t *tarray = (uint8_t *)results->array;
tarray += sizeof(mp_float_t);
larray = (uint8_t *)lhs->array;
larray += sizeof(mp_float_t);
carray_copy_part(tarray, larray, results->shape, lstrides);
} else if(rhs->dtype == NDARRAY_COMPLEX) {
mp_float_t *rarray = (mp_float_t *)rhs->array;
ulab_rescale_float_strides(rstrides);
if(lhs->dtype == NDARRAY_UINT8) {
BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT(results, resarray, uint8_t, larray, lstrides, rarray, rstrides);
} else if(lhs->dtype == NDARRAY_INT8) {
BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT(results, resarray, int8_t, larray, lstrides, rarray, rstrides);
} else if(lhs->dtype == NDARRAY_UINT16) {
BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT(results, resarray, uint16_t, larray, lstrides, rarray, rstrides);
} else if(lhs->dtype == NDARRAY_INT16) {
BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT(results, resarray, int16_t, larray, lstrides, rarray, rstrides);
} else if(lhs->dtype == NDARRAY_FLOAT) {
BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT(results, resarray, mp_float_t, larray, lstrides, rarray, rstrides);
}
}
}
return MP_OBJ_FROM_PTR(results);
}
static void carray_binary_left_divide_(ndarray_obj_t *results, mp_float_t *resarray, uint8_t *larray, uint8_t *rarray,
int32_t *lstrides, int32_t *rstrides, uint8_t rdtype) {
if(rdtype == NDARRAY_UINT8) {
BINARY_LOOP_COMPLEX(results, resarray, uint8_t, larray, lstrides, rarray, rstrides, /);
} else if(rdtype == NDARRAY_INT8) {
BINARY_LOOP_COMPLEX(results, resarray, int8_t, larray, lstrides, rarray, rstrides, /);
} else if(rdtype == NDARRAY_UINT16) {
BINARY_LOOP_COMPLEX(results, resarray, uint16_t, larray, lstrides, rarray, rstrides, /);
} else if(rdtype == NDARRAY_INT16) {
BINARY_LOOP_COMPLEX(results, resarray, int16_t, larray, lstrides, rarray, rstrides, /);
} else if(rdtype == NDARRAY_FLOAT) {
BINARY_LOOP_COMPLEX(results, resarray, mp_float_t, larray, lstrides, rarray, rstrides, /);
}
}
mp_obj_t carray_binary_divide(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {
ndarray_obj_t *results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_COMPLEX);
mp_float_t *resarray = (mp_float_t *)results->array;
if((lhs->dtype == NDARRAY_COMPLEX) && (rhs->dtype == NDARRAY_COMPLEX)) {
mp_float_t *larray = (mp_float_t *)lhs->array;
mp_float_t *rarray = (mp_float_t *)rhs->array;
ulab_rescale_float_strides(lstrides);
ulab_rescale_float_strides(rstrides);
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
// (a + bi) / (c + di) =
// (ac + bd) / (c^2 + d^2) + i (bc - ad) / (c^2 + d^2)
// denominator
mp_float_t denom = rarray[0] * rarray[0] + rarray[1] * rarray[1];
// real part
*resarray++ = (larray[0] * rarray[0] + larray[1] * rarray[1]) / denom;
// imaginary part
*resarray++ = (larray[1] * rarray[0] - larray[0] * rarray[1]) / denom;
larray += lstrides[ULAB_MAX_DIMS - 1];
rarray += rstrides[ULAB_MAX_DIMS - 1];
l++;
} while(l < results->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
larray -= lstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
larray += lstrides[ULAB_MAX_DIMS - 2];
rarray -= rstrides[ULAB_MAX_DIMS - 1] * results->shape[ULAB_MAX_DIMS-1];
rarray += rstrides[ULAB_MAX_DIMS - 2];
k++;
} while(k < results->shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_MAX_DIMS > 2
larray -= lstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
larray += lstrides[ULAB_MAX_DIMS - 3];
rarray -= rstrides[ULAB_MAX_DIMS - 2] * results->shape[ULAB_MAX_DIMS-2];
rarray += rstrides[ULAB_MAX_DIMS - 3];
j++;
} while(j < results->shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS > 2 */
#if ULAB_MAX_DIMS > 3
larray -= lstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
larray += lstrides[ULAB_MAX_DIMS - 4];
rarray -= rstrides[ULAB_MAX_DIMS - 3] * results->shape[ULAB_MAX_DIMS-3];
rarray += rstrides[ULAB_MAX_DIMS - 4];
i++;
} while(i < results->shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS > 3 */
} else {
uint8_t *larray = (uint8_t *)lhs->array;
uint8_t *rarray = (uint8_t *)rhs->array;
if(lhs->dtype == NDARRAY_COMPLEX) {
// real part
carray_binary_left_divide_(results, resarray, larray, rarray, lstrides, rstrides, rhs->dtype);
// imaginary part
resarray = (mp_float_t *)results->array;
resarray++;
larray = (uint8_t *)lhs->array;
larray += sizeof(mp_float_t);
rarray = (uint8_t *)rhs->array;
carray_binary_left_divide_(results, resarray, larray, rarray, lstrides, rstrides, rhs->dtype);
} else {
if(lhs->dtype == NDARRAY_UINT8) {
BINARY_LOOP_COMPLEX_RIGHT_DIVIDE(results, resarray, uint8_t, larray, lstrides, rarray, rstrides);
} else if(lhs->dtype == NDARRAY_INT8) {
BINARY_LOOP_COMPLEX_RIGHT_DIVIDE(results, resarray, int8_t, larray, lstrides, rarray, rstrides);
} else if(lhs->dtype == NDARRAY_UINT16) {
BINARY_LOOP_COMPLEX_RIGHT_DIVIDE(results, resarray, uint16_t, larray, lstrides, rarray, rstrides);
} else if(lhs->dtype == NDARRAY_INT16) {
BINARY_LOOP_COMPLEX_RIGHT_DIVIDE(results, resarray, int16_t, larray, lstrides, rarray, rstrides);
} else if(lhs->dtype == NDARRAY_FLOAT) {
BINARY_LOOP_COMPLEX_RIGHT_DIVIDE(results, resarray, mp_float_t, larray, lstrides, rarray, rstrides);
}
}
}
return MP_OBJ_FROM_PTR(results);
}
#endif

237
python/port/mod/ulab/numpy/carray/carray.h Normal file
View File

@@ -0,0 +1,237 @@
/*
* This file is part of the micropython-ulab project,
*
* https://github.com/v923z/micropython-ulab
*
* The MIT License (MIT)
*
* Copyright (c) 2021-2022 Zoltán Vörös
*/
#ifndef _CARRAY_
#define _CARRAY_
MP_DECLARE_CONST_FUN_OBJ_1(carray_real_obj);
MP_DECLARE_CONST_FUN_OBJ_1(carray_imag_obj);
MP_DECLARE_CONST_FUN_OBJ_1(carray_conjugate_obj);
MP_DECLARE_CONST_FUN_OBJ_1(carray_sort_complex_obj);
mp_obj_t carray_imag(mp_obj_t );
mp_obj_t carray_real(mp_obj_t );
mp_obj_t carray_abs(ndarray_obj_t *, ndarray_obj_t *);
mp_obj_t carray_binary_add(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *);
mp_obj_t carray_binary_multiply(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *);
mp_obj_t carray_binary_subtract(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *);
mp_obj_t carray_binary_divide(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *);
mp_obj_t carray_binary_equal_not_equal(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *, mp_binary_op_t );
#define BINARY_LOOP_COMPLEX1(results, resarray, type_right, larray, lstrides, rarray, rstrides, OPERATOR)\
size_t l = 0;\
do {\
*(resarray) = *((mp_float_t *)(larray)) OPERATOR *((type_right *)(rarray));\
(resarray) += 2;\
(larray) += (lstrides)[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 1]);\
#define BINARY_LOOP_COMPLEX2(results, resarray, type_right, larray, lstrides, rarray, rstrides, OPERATOR)\
size_t k = 0;\
do {\
BINARY_LOOP_COMPLEX1((results), (resarray), type_right, (larray), (lstrides), (rarray), (rstrides), OPERATOR);\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 2]);\
#define BINARY_LOOP_COMPLEX3(results, resarray, type_right, larray, lstrides, rarray, rstrides, OPERATOR)\
size_t j = 0;\
do {\
BINARY_LOOP_COMPLEX2((results), (resarray), type_right, (larray), (lstrides), (rarray), (rstrides), OPERATOR);\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 3];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 3];\
j++;\
} while(j < (results)->shape[ULAB_MAX_DIMS - 3]);\
#define BINARY_LOOP_COMPLEX4(results, resarray, type_right, larray, lstrides, rarray, rstrides, OPERATOR)\
size_t i = 0;\
do {\
BINARY_LOOP_COMPLEX3((results), (resarray), type_right, (larray), (lstrides), (rarray), (rstrides), OPERATOR);\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 4];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 4];\
i++;\
} while(i < (results)->shape[ULAB_MAX_DIMS - 4]);\
#define BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT1(results, resarray, type_left, larray, lstrides, rarray, rstrides)\
size_t l = 0;\
do {\
*(resarray)++ = *((type_left *)(larray)) - (rarray)[0];\
*(resarray)++ = -(rarray)[1];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 1]);\
#define BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT2(results, resarray, type_left, larray, lstrides, rarray, rstrides)\
size_t k = 0;\
do {\
BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT1((results), (resarray), type_left, (larray), (lstrides), (rarray), (rstrides));\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS-1];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS-1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 2]);\
#define BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT3(results, resarray, type_left, larray, lstrides, rarray, rstrides)\
size_t j = 0;\
do {\
BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT2((results), (resarray), type_left, (larray), (lstrides), (rarray), (rstrides));\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 3];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 3];\
j++;\
} while(j < (results)->shape[ULAB_MAX_DIMS - 3]);\
#define BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT4(results, resarray, type_left, larray, lstrides, rarray, rstrides)\
size_t i = 0;\
do {\
BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT3((results), (resarray), type_left, (larray), (lstrides), (rarray), (rstrides));\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 4];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 4];\
i++;\
} while(i < (results)->shape[ULAB_MAX_DIMS - 4]);\
#define BINARY_LOOP_COMPLEX_RIGHT_DIVIDE1(results, resarray, type_left, larray, lstrides, rarray, rstrides)\
size_t l = 0;\
do {\
mp_float_t *c = (mp_float_t *)(rarray);\
mp_float_t denom = c[0] * c[0] + c[1] * c[1];\
mp_float_t a = *((type_left *)(larray)) / denom;\
*(resarray)++ = a * c[0];\
*(resarray)++ = -a * c[1];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 1]);\
#define BINARY_LOOP_COMPLEX_RIGHT_DIVIDE2(results, resarray, type_left, larray, lstrides, rarray, rstrides)\
size_t k = 0;\
do {\
BINARY_LOOP_COMPLEX_RIGHT_DIVIDE1((results), (resarray), type_left, (larray), (lstrides), (rarray), (rstrides));\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 2]);\
#define BINARY_LOOP_COMPLEX_RIGHT_DIVIDE3(results, resarray, type_left, larray, lstrides, rarray, rstrides)\
size_t j = 0;\
do {\
BINARY_LOOP_COMPLEX_RIGHT_DIVIDE2((results), (resarray), type_left, (larray), (lstrides), (rarray), (rstrides));\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 3];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 3];\
j++;\
} while(j < (results)->shape[ULAB_MAX_DIMS - 3]);\
#define BINARY_LOOP_COMPLEX_RIGHT_DIVIDE4(results, resarray, type_left, larray, lstrides, rarray, rstrides)\
size_t i = 0;\
do {\
BINARY_LOOP_COMPLEX_RIGHT_DIVIDE3((results), (resarray), type_left, (larray), (lstrides), (rarray), (rstrides));\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 4];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 4];\
i++;\
} while(i < (results)->shape[ULAB_MAX_DIMS - 4]);\
#define BINARY_LOOP_COMPLEX_EQUAL1(results, array, type_right, larray, lstrides, rarray, rstrides)\
size_t l = 0;\
do {\
if((*(larray) == *((type_right *)(rarray))) && ((larray)[1] == MICROPY_FLOAT_CONST(0.0))) {\
*(array) ^= 0x01;\
}\
(array)++;\
(larray) += (lstrides)[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 1]);\
#define BINARY_LOOP_COMPLEX_EQUAL2(results, array, type_right, larray, lstrides, rarray, rstrides)\
size_t k = 0;\
do {\
BINARY_LOOP_COMPLEX_EQUAL1((results), (array), type_right, (larray), (lstrides), (rarray), (rstrides));\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 2]);\
#define BINARY_LOOP_COMPLEX_EQUAL3(results, array, type_right, larray, lstrides, rarray, rstrides)\
size_t j = 0;\
do {\
BINARY_LOOP_COMPLEX_EQUAL2((results), (array), type_right, (larray), (lstrides), (rarray), (rstrides));\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 3];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 3];\
j++;\
} while(j < (results)->shape[ULAB_MAX_DIMS - 3]);\
#define BINARY_LOOP_COMPLEX_EQUAL4(results, array, type_right, larray, lstrides, rarray, rstrides)\
size_t i = 0;\
do {\
BINARY_LOOP_COMPLEX_EQUAL3((results), (array), type_right, (larray), (lstrides), (rarray), (rstrides));\
(larray) -= (lstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(larray) += (lstrides)[ULAB_MAX_DIMS - 4];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 4];\
i++;\
} while(i < (results)->shape[ULAB_MAX_DIMS - 4]);\
#if ULAB_MAX_DIMS == 1
#define BINARY_LOOP_COMPLEX BINARY_LOOP_COMPLEX1
#define BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT1
#define BINARY_LOOP_COMPLEX_RIGHT_DIVIDE BINARY_LOOP_COMPLEX_RIGHT_DIVIDE1
#define BINARY_LOOP_COMPLEX_EQUAL BINARY_LOOP_COMPLEX_EQUAL1
#endif /* ULAB_MAX_DIMS == 1 */
#if ULAB_MAX_DIMS == 2
#define BINARY_LOOP_COMPLEX BINARY_LOOP_COMPLEX2
#define BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT2
#define BINARY_LOOP_COMPLEX_RIGHT_DIVIDE BINARY_LOOP_COMPLEX_RIGHT_DIVIDE2
#define BINARY_LOOP_COMPLEX_EQUAL BINARY_LOOP_COMPLEX_EQUAL2
#endif /* ULAB_MAX_DIMS == 2 */
#if ULAB_MAX_DIMS == 3
#define BINARY_LOOP_COMPLEX BINARY_LOOP_COMPLEX3
#define BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT3
#define BINARY_LOOP_COMPLEX_RIGHT_DIVIDE BINARY_LOOP_COMPLEX_RIGHT_DIVIDE3
#define BINARY_LOOP_COMPLEX_EQUAL BINARY_LOOP_COMPLEX_EQUAL3
#endif /* ULAB_MAX_DIMS == 3 */
#if ULAB_MAX_DIMS == 4
#define BINARY_LOOP_COMPLEX BINARY_LOOP_COMPLEX4
#define BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT BINARY_LOOP_COMPLEX_REVERSED_SUBTRACT4
#define BINARY_LOOP_COMPLEX_RIGHT_DIVIDE BINARY_LOOP_COMPLEX_RIGHT_DIVIDE4
#define BINARY_LOOP_COMPLEX_EQUAL BINARY_LOOP_COMPLEX_EQUAL4
#endif /* ULAB_MAX_DIMS == 4 */
#endif

28
python/port/mod/ulab/numpy/carray/carray_tools.c Normal file
View File

@@ -0,0 +1,28 @@
/*
* This file is part of the micropython-ulab project,
*
* https://github.com/v923z/micropython-ulab
*
* The MIT License (MIT)
*
* Copyright (c) 2022 Zoltán Vörös
*/
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "py/obj.h"
#include "py/runtime.h"
#include "py/misc.h"
#include "../../ulab.h"
#include "../../ndarray.h"
#if ULAB_SUPPORTS_COMPLEX
void raise_complex_NotImplementedError(void) {
mp_raise_NotImplementedError(translate("not implemented for complex dtype"));
}
#endif

25
python/port/mod/ulab/numpy/carray/carray_tools.h Normal file
View File

@@ -0,0 +1,25 @@
/*
* This file is part of the micropython-ulab project,
*
* https://github.com/v923z/micropython-ulab
*
* The MIT License (MIT)
*
* Copyright (c) 2022 Zoltán Vörös
*/
#ifndef _CARRAY_TOOLS_
#define _CARRAY_TOOLS_
void raise_complex_NotImplementedError(void);
#if ULAB_SUPPORTS_COMPLEX
#define NOT_IMPLEMENTED_FOR_COMPLEX() raise_complex_NotImplementedError();
#define COMPLEX_DTYPE_NOT_IMPLEMENTED(dtype) if((dtype) == NDARRAY_COMPLEX) raise_complex_NotImplementedError();
#else
#define NOT_IMPLEMENTED_FOR_COMPLEX() // do nothing
#define COMPLEX_DTYPE_NOT_IMPLEMENTED(dtype) // do nothing
#endif
#endif

11
python/port/mod/ulab/numpy/compare.c
View File

@@ -20,11 +20,17 @@
#include "../ulab.h"
#include "../ndarray_operators.h"
#include "../ulab_tools.h"
#include "carray/carray_tools.h"
#include "compare.h"
static mp_obj_t compare_function(mp_obj_t x1, mp_obj_t x2, uint8_t op) {
ndarray_obj_t *lhs = ndarray_from_mp_obj(x1, 0);
ndarray_obj_t *rhs = ndarray_from_mp_obj(x2, 0);
#if ULAB_SUPPORTS_COMPLEX
if((lhs->dtype == NDARRAY_COMPLEX) || (rhs->dtype == NDARRAY_COMPLEX)) {
NOT_IMPLEMENTED_FOR_COMPLEX()
}
#endif
uint8_t ndim = 0;
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
int32_t *lstrides = m_new(int32_t, ULAB_MAX_DIMS);
@@ -197,6 +203,7 @@ static mp_obj_t compare_isinf_isfinite(mp_obj_t _x, uint8_t mask) {
}
} else if(mp_obj_is_type(_x, &ulab_ndarray_type)) {
ndarray_obj_t *x = MP_OBJ_TO_PTR(_x);
COMPLEX_DTYPE_NOT_IMPLEMENTED(x->dtype)
ndarray_obj_t *results = ndarray_new_dense_ndarray(x->ndim, x->shape, NDARRAY_BOOL);
// At this point, results is all False
uint8_t *rarray = (uint8_t *)results->array;
@@ -313,6 +320,10 @@ mp_obj_t compare_where(mp_obj_t _condition, mp_obj_t _x, mp_obj_t _y) {
ndarray_obj_t *x = ndarray_from_mp_obj(_x, 0);
ndarray_obj_t *y = ndarray_from_mp_obj(_y, 0);
COMPLEX_DTYPE_NOT_IMPLEMENTED(c->dtype)
COMPLEX_DTYPE_NOT_IMPLEMENTED(x->dtype)
COMPLEX_DTYPE_NOT_IMPLEMENTED(y->dtype)
int32_t *cstrides = m_new(int32_t, ULAB_MAX_DIMS);
int32_t *xstrides = m_new(int32_t, ULAB_MAX_DIMS);
int32_t *ystrides = m_new(int32_t, ULAB_MAX_DIMS);

843
python/port/mod/ulab/numpy/create.c Normal file
View File

@@ -0,0 +1,843 @@
/*
* 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
* 2019-2021 Zoltán Vörös
* 2020 Taku Fukada
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "py/obj.h"
#include "py/runtime.h"
#include "create.h"
#include "../ulab.h"
#include "../ulab_tools.h"
#if ULAB_NUMPY_HAS_ONES | ULAB_NUMPY_HAS_ZEROS | ULAB_NUMPY_HAS_FULL | ULAB_NUMPY_HAS_EMPTY
static mp_obj_t create_zeros_ones_full(mp_obj_t oshape, uint8_t dtype, mp_obj_t value) {
if(!mp_obj_is_int(oshape) && !mp_obj_is_type(oshape, &mp_type_tuple) && !mp_obj_is_type(oshape, &mp_type_list)) {
mp_raise_TypeError(translate("input argument must be an integer, a tuple, or a list"));
}
ndarray_obj_t *ndarray = NULL;
if(mp_obj_is_int(oshape)) {
size_t n = mp_obj_get_int(oshape);
ndarray = ndarray_new_linear_array(n, dtype);
} else if(mp_obj_is_type(oshape, &mp_type_tuple) || mp_obj_is_type(oshape, &mp_type_list)) {
uint8_t len = (uint8_t)mp_obj_get_int(mp_obj_len_maybe(oshape));
if(len > ULAB_MAX_DIMS) {
mp_raise_TypeError(translate("too many dimensions"));
}
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
size_t i = 0;
mp_obj_iter_buf_t iter_buf;
mp_obj_t item, iterable = mp_getiter(oshape, &iter_buf);
while((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION){
shape[ULAB_MAX_DIMS - len + i] = (size_t)mp_obj_get_int(item);
i++;
}
ndarray = ndarray_new_dense_ndarray(len, shape, dtype);
}
if(value != mp_const_none) {
if(dtype == NDARRAY_BOOL) {
dtype = NDARRAY_UINT8;
if(mp_obj_is_true(value)) {
value = mp_obj_new_int(1);
} else {
value = mp_obj_new_int(0);
}
}
for(size_t i=0; i < ndarray->len; i++) {
#if ULAB_SUPPORTS_COMPLEX
if(dtype == NDARRAY_COMPLEX) {
ndarray_set_complex_value(ndarray->array, i, value);
} else {
ndarray_set_value(dtype, ndarray->array, i, value);
}
#else
ndarray_set_value(dtype, ndarray->array, i, value);
#endif
}
}
// if zeros calls the function, we don't have to do anything
return MP_OBJ_FROM_PTR(ndarray);
}
#endif
#if ULAB_NUMPY_HAS_ARANGE | ULAB_NUMPY_HAS_LINSPACE
static ndarray_obj_t *create_linspace_arange(mp_float_t start, mp_float_t step, mp_float_t stop, size_t len, uint8_t dtype) {
mp_float_t value = start;
ndarray_obj_t *ndarray = ndarray_new_linear_array(len, dtype);
if(ndarray->boolean == NDARRAY_BOOLEAN) {
uint8_t *array = (uint8_t *)ndarray->array;
for(size_t i=0; i < len; i++, value += step) {
*array++ = value == MICROPY_FLOAT_CONST(0.0) ? 0 : 1;
}
} else if(dtype == NDARRAY_UINT8) {
ARANGE_LOOP(uint8_t, ndarray, len, step, stop);
} else if(dtype == NDARRAY_INT8) {
ARANGE_LOOP(int8_t, ndarray, len, step, stop);
} else if(dtype == NDARRAY_UINT16) {
ARANGE_LOOP(uint16_t, ndarray, len, step, stop);
} else if(dtype == NDARRAY_INT16) {
ARANGE_LOOP(int16_t, ndarray, len, step, stop);
} else {
ARANGE_LOOP(mp_float_t, ndarray, len, step, stop);
}
return ndarray;
}
#endif
#if ULAB_NUMPY_HAS_ARANGE
//| @overload
//| def arange(stop: _float, step: _float = 1, *, dtype: _DType = ulab.numpy.float) -> ulab.numpy.ndarray: ...
//| @overload
//| def arange(start: _float, stop: _float, step: _float = 1, *, dtype: _DType = ulab.numpy.float) -> ulab.numpy.ndarray:
//| """
//| .. param: start
//| First value in the array, optional, defaults to 0
//| .. param: stop
//| Final value in the array
//| .. param: step
//| Difference between consecutive elements, optional, defaults to 1.0
//| .. param: dtype
//| Type of values in the array
//|
//| Return a new 1-D array with elements ranging from ``start`` to ``stop``, with step size ``step``."""
//| ...
//|
mp_obj_t create_arange(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_const_none } },
{ MP_QSTR_, MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_, MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_dtype, 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);
uint8_t dtype = NDARRAY_FLOAT;
mp_float_t start, stop, step;
if(n_args == 1) {
start = MICROPY_FLOAT_CONST(0.0);
stop = mp_obj_get_float(args[0].u_obj);
step = MICROPY_FLOAT_CONST(1.0);
if(mp_obj_is_int(args[0].u_obj)) dtype = NDARRAY_INT16;
} else if(n_args == 2) {
start = mp_obj_get_float(args[0].u_obj);
stop = mp_obj_get_float(args[1].u_obj);
step = MICROPY_FLOAT_CONST(1.0);
if(mp_obj_is_int(args[0].u_obj) && mp_obj_is_int(args[1].u_obj)) dtype = NDARRAY_INT16;
} else if(n_args == 3) {
start = mp_obj_get_float(args[0].u_obj);
stop = mp_obj_get_float(args[1].u_obj);
step = mp_obj_get_float(args[2].u_obj);
if(mp_obj_is_int(args[0].u_obj) && mp_obj_is_int(args[1].u_obj) && mp_obj_is_int(args[2].u_obj)) dtype = NDARRAY_INT16;
} else {
mp_raise_TypeError(translate("wrong number of arguments"));
}
if((MICROPY_FLOAT_C_FUN(fabs)(stop) > 32768) || (MICROPY_FLOAT_C_FUN(fabs)(start) > 32768) || (MICROPY_FLOAT_C_FUN(fabs)(step) > 32768)) {
dtype = NDARRAY_FLOAT;
}
if(args[3].u_obj != mp_const_none) {
dtype = (uint8_t)mp_obj_get_int(args[3].u_obj);
}
ndarray_obj_t *ndarray;
if((stop - start)/step < 0) {
ndarray = ndarray_new_linear_array(0, dtype);
} else {
size_t len = (size_t)(MICROPY_FLOAT_C_FUN(ceil)((stop - start) / step));
stop = start + (len - 1) * step;
ndarray = create_linspace_arange(start, step, stop, len, dtype);
}
return MP_OBJ_FROM_PTR(ndarray);
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_arange_obj, 1, create_arange);
#endif
#if ULAB_NUMPY_HAS_ASARRAY
mp_obj_t create_asarray(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_const_none } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_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);
uint8_t _dtype;
#if ULAB_HAS_DTYPE_OBJECT
if(mp_obj_is_type(args[1].u_obj, &ulab_dtype_type)) {
dtype_obj_t *dtype = MP_OBJ_TO_PTR(args[1].u_obj);
_dtype = dtype->dtype;
} else { // this must be an integer defined as a class constant (ulab.numpy.uint8 etc.)
if(args[1].u_obj == mp_const_none) {
_dtype = 0;
} else {
_dtype = mp_obj_get_int(args[1].u_obj);
}
}
#else
if(args[1].u_obj == mp_const_none) {
_dtype = 0;
} else {
_dtype = mp_obj_get_int(args[1].u_obj);
}
#endif
if(ulab_tools_mp_obj_is_scalar(args[0].u_obj)) {
return args[0].u_obj;
} else if(mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type)) {
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
if((_dtype == ndarray->dtype) || (_dtype == 0)) {
return args[0].u_obj;
} else {
return MP_OBJ_FROM_PTR(ndarray_copy_view_convert_type(ndarray, _dtype));
}
} else if(ndarray_object_is_array_like(args[0].u_obj)) {
if(_dtype == 0) {
_dtype = NDARRAY_FLOAT;
}
return MP_OBJ_FROM_PTR(ndarray_from_iterable(args[0].u_obj, _dtype));
} else {
mp_raise_TypeError(translate("wrong input type"));
}
return mp_const_none; // this should never happen
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_asarray_obj, 1, create_asarray);
#endif
#if ULAB_NUMPY_HAS_CONCATENATE
//| def concatenate(arrays: Tuple[ulab.numpy.ndarray], *, axis: int = 0) -> ulab.numpy.ndarray:
//| """
//| .. param: arrays
//| tuple of ndarrays
//| .. param: axis
//| axis along which the arrays will be joined
//|
//| Join a sequence of arrays along an existing axis."""
//| ...
//|
mp_obj_t create_concatenate(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_const_none } },
{ MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = 0 } },
};
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, &mp_type_tuple)) {
mp_raise_TypeError(translate("first argument must be a tuple of ndarrays"));
}
int8_t axis = (int8_t)args[1].u_int;
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
mp_obj_tuple_t *ndarrays = MP_OBJ_TO_PTR(args[0].u_obj);
// first check, whether the arrays are compatible
ndarray_obj_t *_ndarray = MP_OBJ_TO_PTR(ndarrays->items[0]);
uint8_t dtype = _ndarray->dtype;
uint8_t ndim = _ndarray->ndim;
if(axis < 0) {
axis += ndim;
}
if((axis < 0) || (axis >= ndim)) {
mp_raise_ValueError(translate("wrong axis specified"));
}
// shift axis
axis = ULAB_MAX_DIMS - ndim + axis;
for(uint8_t j=0; j < ULAB_MAX_DIMS; j++) {
shape[j] = _ndarray->shape[j];
}
for(uint8_t i=1; i < ndarrays->len; i++) {
_ndarray = MP_OBJ_TO_PTR(ndarrays->items[i]);
// check, whether the arrays are compatible
if((dtype != _ndarray->dtype) || (ndim != _ndarray->ndim)) {
mp_raise_ValueError(translate("input arrays are not compatible"));
}
for(uint8_t j=0; j < ULAB_MAX_DIMS; j++) {
if(j == axis) {
shape[j] += _ndarray->shape[j];
} else {
if(shape[j] != _ndarray->shape[j]) {
mp_raise_ValueError(translate("input arrays are not compatible"));
}
}
}
}
ndarray_obj_t *target = ndarray_new_dense_ndarray(ndim, shape, dtype);
uint8_t *tpos = (uint8_t *)target->array;
uint8_t *tarray;
for(uint8_t p=0; p < ndarrays->len; p++) {
// reset the pointer along the axis
ndarray_obj_t *source = MP_OBJ_TO_PTR(ndarrays->items[p]);
uint8_t *sarray = (uint8_t *)source->array;
tarray = tpos;
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
memcpy(tarray, sarray, source->itemsize);
tarray += target->strides[ULAB_MAX_DIMS - 1];
sarray += source->strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < source->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
tarray -= target->strides[ULAB_MAX_DIMS - 1] * source->shape[ULAB_MAX_DIMS-1];
tarray += target->strides[ULAB_MAX_DIMS - 2];
sarray -= source->strides[ULAB_MAX_DIMS - 1] * source->shape[ULAB_MAX_DIMS-1];
sarray += source->strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < source->shape[ULAB_MAX_DIMS - 2]);
#endif
#if ULAB_MAX_DIMS > 2
tarray -= target->strides[ULAB_MAX_DIMS - 2] * source->shape[ULAB_MAX_DIMS-2];
tarray += target->strides[ULAB_MAX_DIMS - 3];
sarray -= source->strides[ULAB_MAX_DIMS - 2] * source->shape[ULAB_MAX_DIMS-2];
sarray += source->strides[ULAB_MAX_DIMS - 3];
j++;
} while(j < source->shape[ULAB_MAX_DIMS - 3]);
#endif
#if ULAB_MAX_DIMS > 3
tarray -= target->strides[ULAB_MAX_DIMS - 3] * source->shape[ULAB_MAX_DIMS-3];
tarray += target->strides[ULAB_MAX_DIMS - 4];
sarray -= source->strides[ULAB_MAX_DIMS - 3] * source->shape[ULAB_MAX_DIMS-3];
sarray += source->strides[ULAB_MAX_DIMS - 4];
i++;
} while(i < source->shape[ULAB_MAX_DIMS - 4]);
#endif
if(p < ndarrays->len - 1) {
tpos += target->strides[axis] * source->shape[axis];
}
}
return MP_OBJ_FROM_PTR(target);
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_concatenate_obj, 1, create_concatenate);
#endif
#if ULAB_MAX_DIMS > 1
#if ULAB_NUMPY_HAS_DIAG
//| def diag(a: ulab.numpy.ndarray, *, k: int = 0) -> ulab.numpy.ndarray:
//| """
//| .. param: a
//| an ndarray
//| .. param: k
//| Offset of the diagonal from the main diagonal. Can be positive or negative.
//|
//| Return specified diagonals."""
//| ...
//|
mp_obj_t create_diag(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_const_none } },
{ MP_QSTR_k, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = 0 } },
};
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);
ndarray_obj_t *source = ndarray_from_iterable(args[0].u_obj, NDARRAY_FLOAT);
ndarray_obj_t *target = NULL;
int32_t k = args[1].u_int;
size_t k_abs = k >= 0 ? (size_t)k : (size_t)(-k);
if(source->ndim == 2) { // return the diagonal
size_t len;
if(k >= 0) {
len = (k_abs <= source->shape[ULAB_MAX_DIMS - 1]) ? source->shape[ULAB_MAX_DIMS - 1] - k_abs : 0;
} else {
len = (k_abs <= source->shape[ULAB_MAX_DIMS - 2]) ? source->shape[ULAB_MAX_DIMS - 2] - k_abs : 0;
}
target = ndarray_new_linear_array(len, source->dtype);
if(len == 0) {
return MP_OBJ_FROM_PTR(target);
}
uint8_t *sarray = (uint8_t *)source->array;
uint8_t *tarray = (uint8_t *)target->array;
if(k >= 0) {
sarray += source->strides[ULAB_MAX_DIMS - 1] * k;
} else {
sarray += source->strides[ULAB_MAX_DIMS - 2] * k_abs;
}
for(size_t i=0; i < len; i++) {
memcpy(tarray, sarray, source->itemsize);
sarray += (source->strides[ULAB_MAX_DIMS - 1] + source->strides[ULAB_MAX_DIMS - 2]);
tarray += target->itemsize;
}
} else if(source->ndim == 1) { // return a rank-2 tensor with the prescribed diagonal
size_t len = source->len + k_abs;
target = ndarray_new_dense_ndarray(2, ndarray_shape_vector(0, 0, len, len), source->dtype);
uint8_t *sarray = (uint8_t *)source->array;
uint8_t *tarray = (uint8_t *)target->array;
if(k < 0) {
tarray += len * k_abs * target->itemsize;
} else {
tarray += k_abs * target->itemsize;
}
for(size_t i = 0; i < source->len; i++) {
memcpy(tarray, sarray, source->itemsize);
sarray += source->strides[ULAB_MAX_DIMS - 1];
tarray += (len + 1) * target->itemsize;
}
}
#if ULAB_MAX_DIMS > 2
else {
mp_raise_ValueError(translate("input must be 1- or 2-d"));
}
#endif
return MP_OBJ_FROM_PTR(target);
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_diag_obj, 1, create_diag);
#endif /* ULAB_NUMPY_HAS_DIAG */
#if ULAB_NUMPY_HAS_EMPTY
// This function is bound in numpy.c to numpy.zeros(), and is simply an alias for that
//| def empty(shape: Union[int, Tuple[int, ...]], *, dtype: _DType = ulab.numpy.float) -> ulab.numpy.ndarray:
//| """
//| .. param: shape
//| Shape of the array, either an integer (for a 1-D array) or a tuple of 2 integers (for a 2-D array)
//| .. param: dtype
//| Type of values in the array
//|
//| Return a new array of the given shape with all elements set to 0. An alias for numpy.zeros."""
//| ...
//|
#endif
#if ULAB_NUMPY_HAS_EYE
//| def eye(size: int, *, M: Optional[int] = None, k: int = 0, dtype: _DType = ulab.numpy.float) -> ulab.numpy.ndarray:
//| """Return a new square array of size, with the diagonal elements set to 1
//| and the other elements set to 0. If k is given, the diagonal is shifted by the specified amount."""
//| ...
//|
mp_obj_t create_eye(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_INT, { .u_int = 0 } },
{ MP_QSTR_M, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_k, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = 0 } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = NDARRAY_FLOAT } },
};
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);
size_t n = args[0].u_int, m;
size_t k = args[2].u_int > 0 ? (size_t)args[2].u_int : (size_t)(-args[2].u_int);
uint8_t dtype = args[3].u_int;
if(args[1].u_rom_obj == mp_const_none) {
m = n;
} else {
m = mp_obj_get_int(args[1].u_rom_obj);
}
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(2, ndarray_shape_vector(0, 0, n, m), dtype);
if(dtype == NDARRAY_BOOL) {
dtype = NDARRAY_UINT8;
}
mp_obj_t one = mp_obj_new_int(1);
size_t i = 0;
if((args[2].u_int >= 0)) {
while(k < m) {
ndarray_set_value(dtype, ndarray->array, i*m+k, one);
k++;
i++;
}
} else {
while(k < n) {
ndarray_set_value(dtype, ndarray->array, k*m+i, one);
k++;
i++;
}
}
return MP_OBJ_FROM_PTR(ndarray);
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_eye_obj, 1, create_eye);
#endif /* ULAB_NUMPY_HAS_EYE */
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_NUMPY_HAS_FULL
//| def full(shape: Union[int, Tuple[int, ...]], fill_value: Union[_float, _bool], *, dtype: _DType = ulab.numpy.float) -> ulab.numpy.ndarray:
//| """
//| .. param: shape
//| Shape of the array, either an integer (for a 1-D array) or a tuple of integers (for tensors of higher rank)
//| .. param: fill_value
//| scalar, the value with which the array is filled
//| .. param: dtype
//| Type of values in the array
//|
//| Return a new array of the given shape with all elements set to 0."""
//| ...
//|
mp_obj_t create_full(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_obj = MP_OBJ_NULL } },
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_obj = MP_OBJ_NULL } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = NDARRAY_FLOAT } },
};
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);
uint8_t dtype = args[2].u_int;
return create_zeros_ones_full(args[0].u_obj, dtype, args[1].u_obj);
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_full_obj, 0, create_full);
#endif
#if ULAB_NUMPY_HAS_LINSPACE
//| def linspace(
//| start: _float,
//| stop: _float,
//| *,
//| dtype: _DType = ulab.numpy.float,
//| num: int = 50,
//| endpoint: _bool = True,
//| retstep: _bool = False
//| ) -> ulab.numpy.ndarray:
//| """
//| .. param: start
//| First value in the array
//| .. param: stop
//| Final value in the array
//| .. param int: num
//| Count of values in the array.
//| .. param: dtype
//| Type of values in the array
//| .. param bool: endpoint
//| Whether the ``stop`` value is included. Note that even when
//| endpoint=True, the exact ``stop`` value may not be included due to the
//| inaccuracy of floating point arithmetic.
//| .. param bool: retstep,
//| If True, return (`samples`, `step`), where `step` is the spacing between samples.
//|
//| Return a new 1-D array with ``num`` elements ranging from ``start`` to ``stop`` linearly."""
//| ...
//|
mp_obj_t create_linspace(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_const_none } },
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_num, MP_ARG_INT, { .u_int = 50 } },
{ MP_QSTR_endpoint, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_true } },
{ MP_QSTR_retstep, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_false } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = NDARRAY_FLOAT } },
};
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(args[2].u_int < 2) {
mp_raise_ValueError(translate("number of points must be at least 2"));
}
size_t len = (size_t)args[2].u_int;
mp_float_t start, step, stop;
ndarray_obj_t *ndarray = NULL;
#if ULAB_SUPPORTS_COMPLEX
mp_float_t step_real, step_imag;
bool complex_out = false;
if(mp_obj_is_type(args[0].u_obj, &mp_type_complex) || mp_obj_is_type(args[1].u_obj, &mp_type_complex)) {
complex_out = true;
ndarray = ndarray_new_linear_array(len, NDARRAY_COMPLEX);
mp_float_t *array = (mp_float_t *)ndarray->array;
mp_float_t start_real, start_imag;
mp_float_t stop_real, stop_imag;
mp_obj_get_complex(args[0].u_obj, &start_real, &start_imag);
mp_obj_get_complex(args[1].u_obj, &stop_real, &stop_imag);
if(args[3].u_obj == mp_const_true) {
step_real = (stop_real - start_real) / (len - 1);
step_imag = (stop_imag - start_imag) / (len - 1);
} else {
step_real = (stop_real - start_real) / len;
step_imag = (stop_imag - start_imag) / len;
}
for(size_t i = 0; i < len; i++) {
*array++ = start_real;
*array++ = start_imag;
start_real += step_real;
start_imag += step_imag;
}
} else {
#endif
start = mp_obj_get_float(args[0].u_obj);
stop = mp_obj_get_float(args[1].u_obj);
uint8_t typecode = args[5].u_int;
if(args[3].u_obj == mp_const_true) {
step = (stop - start) / (len - 1);
} else {
step = (stop - start) / len;
stop = start + step * (len - 1);
}
ndarray = create_linspace_arange(start, step, stop, len, typecode);
#if ULAB_SUPPORTS_COMPLEX
}
#endif
if(args[4].u_obj == mp_const_false) {
return MP_OBJ_FROM_PTR(ndarray);
} else {
mp_obj_t tuple[2];
tuple[0] = ndarray;
#if ULAB_SUPPORTS_COMPLEX
if(complex_out) {
tuple[1] = mp_obj_new_complex(step_real, step_imag);
} else {
tuple[1] = mp_obj_new_float(step);
}
#else /* ULAB_SUPPORTS_COMPLEX */
tuple[1] = mp_obj_new_float(step);
#endif
return mp_obj_new_tuple(2, tuple);
}
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_linspace_obj, 2, create_linspace);
#endif
#if ULAB_NUMPY_HAS_LOGSPACE
//| def logspace(
//| start: _float,
//| stop: _float,
//| *,
//| dtype: _DType = ulab.numpy.float,
//| num: int = 50,
//| endpoint: _bool = True,
//| base: _float = 10.0
//| ) -> ulab.numpy.ndarray:
//| """
//| .. param: start
//| First value in the array
//| .. param: stop
//| Final value in the array
//| .. param int: num
//| Count of values in the array. Defaults to 50.
//| .. param: base
//| The base of the log space. The step size between the elements in
//| ``ln(samples) / ln(base)`` (or ``log_base(samples)``) is uniform. Defaults to 10.0.
//| .. param: dtype
//| Type of values in the array
//| .. param bool: endpoint
//| Whether the ``stop`` value is included. Note that even when
//| endpoint=True, the exact ``stop`` value may not be included due to the
//| inaccuracy of floating point arithmetic. Defaults to True.
//|
//| Return a new 1-D array with ``num`` evenly spaced elements on a log scale.
//| The sequence starts at ``base ** start``, and ends with ``base ** stop``."""
//| ...
//|
const mp_obj_float_t create_float_const_ten = {{&mp_type_float}, MICROPY_FLOAT_CONST(10.0)};
mp_obj_t create_logspace(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_const_none } },
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_num, MP_ARG_INT, { .u_int = 50 } },
{ MP_QSTR_base, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_PTR(&create_float_const_ten) } },
{ MP_QSTR_endpoint, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_true } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = NDARRAY_FLOAT } },
};
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(args[2].u_int < 2) {
mp_raise_ValueError(translate("number of points must be at least 2"));
}
size_t len = (size_t)args[2].u_int;
mp_float_t start, step, quotient;
start = mp_obj_get_float(args[0].u_obj);
uint8_t dtype = args[5].u_int;
mp_float_t base = mp_obj_get_float(args[3].u_obj);
if(args[4].u_obj == mp_const_true) step = (mp_obj_get_float(args[1].u_obj) - start)/(len - 1);
else step = (mp_obj_get_float(args[1].u_obj) - start) / len;
quotient = MICROPY_FLOAT_C_FUN(pow)(base, step);
ndarray_obj_t *ndarray = ndarray_new_linear_array(len, dtype);
mp_float_t value = MICROPY_FLOAT_C_FUN(pow)(base, start);
if(ndarray->dtype == NDARRAY_UINT8) {
uint8_t *array = (uint8_t *)ndarray->array;
if(ndarray->boolean) {
memset(array, 1, len);
} else {
for(size_t i=0; i < len; i++, value *= quotient) *array++ = (uint8_t)value;
}
} else if(ndarray->dtype == NDARRAY_INT8) {
int8_t *array = (int8_t *)ndarray->array;
for(size_t i=0; i < len; i++, value *= quotient) *array++ = (int8_t)value;
} else if(ndarray->dtype == NDARRAY_UINT16) {
uint16_t *array = (uint16_t *)ndarray->array;
for(size_t i=0; i < len; i++, value *= quotient) *array++ = (uint16_t)value;
} else if(ndarray->dtype == NDARRAY_INT16) {
int16_t *array = (int16_t *)ndarray->array;
for(size_t i=0; i < len; i++, value *= quotient) *array++ = (int16_t)value;
} else {
mp_float_t *array = (mp_float_t *)ndarray->array;
for(size_t i=0; i < len; i++, value *= quotient) *array++ = value;
}
return MP_OBJ_FROM_PTR(ndarray);
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_logspace_obj, 2, create_logspace);
#endif
#if ULAB_NUMPY_HAS_ONES
//| def ones(shape: Union[int, Tuple[int, ...]], *, dtype: _DType = ulab.numpy.float) -> ulab.numpy.ndarray:
//| """
//| .. param: shape
//| Shape of the array, either an integer (for a 1-D array) or a tuple of 2 integers (for a 2-D array)
//| .. param: dtype
//| Type of values in the array
//|
//| Return a new array of the given shape with all elements set to 1."""
//| ...
//|
mp_obj_t create_ones(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_obj = MP_OBJ_NULL } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = NDARRAY_FLOAT } },
};
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);
uint8_t dtype = args[1].u_int;
mp_obj_t one = mp_obj_new_int(1);
return create_zeros_ones_full(args[0].u_obj, dtype, one);
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_ones_obj, 0, create_ones);
#endif
#if ULAB_NUMPY_HAS_ZEROS
//| def zeros(shape: Union[int, Tuple[int, ...]], *, dtype: _DType = ulab.numpy.float) -> ulab.numpy.ndarray:
//| """
//| .. param: shape
//| Shape of the array, either an integer (for a 1-D array) or a tuple of 2 integers (for a 2-D array)
//| .. param: dtype
//| Type of values in the array
//|
//| Return a new array of the given shape with all elements set to 0."""
//| ...
//|
mp_obj_t create_zeros(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_obj = MP_OBJ_NULL } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = NDARRAY_FLOAT } },
};
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);
uint8_t dtype = args[1].u_int;
return create_zeros_ones_full(args[0].u_obj, dtype, mp_const_none);
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_zeros_obj, 0, create_zeros);
#endif
#if ULAB_NUMPY_HAS_FROMBUFFER
mp_obj_t create_frombuffer(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_const_none } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_INT(NDARRAY_FLOAT) } },
{ MP_QSTR_count, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_INT(-1) } },
{ MP_QSTR_offset, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_INT(0) } },
};
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);
uint8_t dtype = mp_obj_get_int(args[1].u_obj);
size_t offset = mp_obj_get_int(args[3].u_obj);
mp_buffer_info_t bufinfo;
if(mp_get_buffer(args[0].u_obj, &bufinfo, MP_BUFFER_READ)) {
size_t sz = ulab_binary_get_size(dtype);
if(bufinfo.len < offset) {
mp_raise_ValueError(translate("offset must be non-negative and no greater than buffer length"));
}
size_t len = (bufinfo.len - offset) / sz;
if((len * sz) != (bufinfo.len - offset)) {
mp_raise_ValueError(translate("buffer size must be a multiple of element size"));
}
if(mp_obj_get_int(args[2].u_obj) > 0) {
size_t count = mp_obj_get_int(args[2].u_obj);
if(len < count) {
mp_raise_ValueError(translate("buffer is smaller than requested size"));
} else {
len = count;
}
}
ndarray_obj_t *ndarray = m_new_obj(ndarray_obj_t);
ndarray->base.type = &ulab_ndarray_type;
ndarray->dtype = dtype == NDARRAY_BOOL ? NDARRAY_UINT8 : dtype;
ndarray->boolean = dtype == NDARRAY_BOOL ? NDARRAY_BOOLEAN : NDARRAY_NUMERIC;
ndarray->ndim = 1;
ndarray->len = len;
ndarray->itemsize = sz;
ndarray->shape[ULAB_MAX_DIMS - 1] = len;
ndarray->strides[ULAB_MAX_DIMS - 1] = sz;
uint8_t *buffer = bufinfo.buf;
ndarray->array = buffer + offset;
return MP_OBJ_FROM_PTR(ndarray);
}
return mp_const_none;
}
MP_DEFINE_CONST_FUN_OBJ_KW(create_frombuffer_obj, 1, create_frombuffer);
#endif

84
python/port/mod/ulab/numpy/create.h Normal file
View File

@@ -0,0 +1,84 @@
/*
* 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
* 2019-2021 Zoltán Vörös
*/
#ifndef _CREATE_
#define _CREATE_
#include "../ulab.h"
#include "../ndarray.h"
#if ULAB_NUMPY_HAS_ARANGE
mp_obj_t create_arange(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_arange_obj);
#endif
#if ULAB_NUMPY_HAS_ASARRAY
mp_obj_t create_arange(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_asarray_obj);
#endif
#if ULAB_NUMPY_HAS_CONCATENATE
mp_obj_t create_concatenate(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_concatenate_obj);
#endif
#if ULAB_NUMPY_HAS_DIAG
mp_obj_t create_diag(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_diag_obj);
#endif
#if ULAB_MAX_DIMS > 1
#if ULAB_NUMPY_HAS_EYE
mp_obj_t create_eye(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_eye_obj);
#endif
#endif
#if ULAB_NUMPY_HAS_FULL
mp_obj_t create_full(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_full_obj);
#endif
#if ULAB_NUMPY_HAS_LINSPACE
mp_obj_t create_linspace(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_linspace_obj);
#endif
#if ULAB_NUMPY_HAS_LOGSPACE
mp_obj_t create_logspace(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_logspace_obj);
#endif
#if ULAB_NUMPY_HAS_ONES
mp_obj_t create_ones(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_ones_obj);
#endif
#if ULAB_NUMPY_HAS_ZEROS
mp_obj_t create_zeros(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_zeros_obj);
#endif
#if ULAB_NUMPY_HAS_FROMBUFFER
mp_obj_t create_frombuffer(size_t , const mp_obj_t *, mp_map_t *);
MP_DECLARE_CONST_FUN_OBJ_KW(create_frombuffer_obj);
#endif
#define ARANGE_LOOP(type_, ndarray, len, step, stop) \
({\
type_ *array = (type_ *)(ndarray)->array;\
for (size_t i = 0; i < (len) - 1; i++, (value) += (step)) {\
*array++ = (type_)(value);\
}\
*array = (type_)(stop);\
})
#endif

23
python/port/mod/ulab/numpy/fft/fft.c
View File

@@ -20,11 +20,13 @@
#include "py/obj.h"
#include "py/objarray.h"
#include "../carray/carray_tools.h"
#include "fft.h"
//| """Frequency-domain functions"""
//|
//| import ulab.numpy
//| import ulab.utils
//| def fft(r: ulab.numpy.ndarray, c: Optional[ulab.numpy.ndarray] = None) -> Tuple[ulab.numpy.ndarray, ulab.numpy.ndarray]:
@@ -35,10 +37,17 @@
//|
//| Perform a Fast Fourier Transform from the time domain into the frequency domain
//|
//| See also ~ulab.extras.spectrum, which computes the magnitude of the fft,
//| See also `ulab.utils.spectrogram`, which computes the magnitude of the fft,
//| rather than separately returning its real and imaginary parts."""
//| ...
//|
#if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE
static mp_obj_t fft_fft(mp_obj_t arg) {
return fft_fft_ifft_spectrogram(arg, FFT_FFT);
}
MP_DEFINE_CONST_FUN_OBJ_1(fft_fft_obj, fft_fft);
#else
static mp_obj_t fft_fft(size_t n_args, const mp_obj_t *args) {
if(n_args == 2) {
return fft_fft_ifft_spectrogram(n_args, args[0], args[1], FFT_FFT);
@@ -48,6 +57,7 @@ static mp_obj_t fft_fft(size_t n_args, const mp_obj_t *args) {
}
MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj, 1, 2, fft_fft);
#endif
//| def ifft(r: ulab.numpy.ndarray, c: Optional[ulab.numpy.ndarray] = None) -> Tuple[ulab.numpy.ndarray, ulab.numpy.ndarray]:
//| """
@@ -55,11 +65,19 @@ MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj, 1, 2, fft_fft);
//| :param ulab.numpy.ndarray c: An optional 1-dimension array of values whose size is a power of 2, giving the complex part of the value
//| :return tuple (r, c): The real and complex parts of the inverse FFT
//|
//| Perform an Inverse Fast Fourier Transform from the frequency domain into the time domain"""
//| Perform an Inverse Fast Fourier Transform from the frequeny domain into the time domain"""
//| ...
//|
#if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE
static mp_obj_t fft_ifft(mp_obj_t arg) {
return fft_fft_ifft_spectrogram(arg, FFT_IFFT);
}
MP_DEFINE_CONST_FUN_OBJ_1(fft_ifft_obj, fft_ifft);
#else
static mp_obj_t fft_ifft(size_t n_args, const mp_obj_t *args) {
NOT_IMPLEMENTED_FOR_COMPLEX()
if(n_args == 2) {
return fft_fft_ifft_spectrogram(n_args, args[0], args[1], FFT_IFFT);
} else {
@@ -68,6 +86,7 @@ static mp_obj_t fft_ifft(size_t n_args, const mp_obj_t *args) {
}
MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_ifft_obj, 1, 2, fft_ifft);
#endif
STATIC const mp_rom_map_elem_t ulab_fft_globals_table[] = {
{ MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_fft) },

6
python/port/mod/ulab/numpy/fft/fft.h
View File

@@ -19,6 +19,12 @@
extern mp_obj_module_t ulab_fft_module;
#if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE
MP_DECLARE_CONST_FUN_OBJ_3(fft_fft_obj);
MP_DECLARE_CONST_FUN_OBJ_3(fft_ifft_obj);
#else
MP_DECLARE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj);
MP_DECLARE_CONST_FUN_OBJ_VAR_BETWEEN(fft_ifft_obj);
#endif
#endif

136
python/port/mod/ulab/numpy/fft/fft_tools.c
View File

@@ -9,10 +9,12 @@
*/
#include <math.h>
#include <string.h>
#include "py/runtime.h"
#include "../../ndarray.h"
#include "../../ulab_tools.h"
#include "../carray/carray_tools.h"
#include "fft_tools.h"
#ifndef MP_PI
@@ -22,7 +24,8 @@
#define MP_E MICROPY_FLOAT_CONST(2.71828182845904523536)
#endif
/*
/* Kernel implementation for the case, when ulab has no complex support
* The following function takes two arrays, namely, the real and imaginary
* parts of a complex array, and calculates the Fourier transform in place.
*
@@ -31,6 +34,128 @@
* and can be used independent of ulab.
*/
#if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE
/* Kernel implementation for the complex case. Data are contained in data as
data[0], data[1], data[2], data[3], .... , data[2n - 2], data[2n-1]
real[0], imag[0], real[1], imag[1], .... , real[n-1], imag[n-1]
In general
real[i] = data[2i]
imag[i] = data[2i+1]
*/
void fft_kernel_complex(mp_float_t *data, size_t n, int isign) {
size_t j, m, mmax, istep;
mp_float_t tempr, tempi;
mp_float_t wtemp, wr, wpr, wpi, wi, theta;
j = 0;
for(size_t i = 0; i < n; i++) {
if (j > i) {
SWAP(mp_float_t, data[2*i], data[2*j]);
SWAP(mp_float_t, data[2*i+1], data[2*j+1]);
}
m = n >> 1;
while (j >= m && m > 0) {
j -= m;
m >>= 1;
}
j += m;
}
mmax = 1;
while (n > mmax) {
istep = mmax << 1;
theta = MICROPY_FLOAT_CONST(-2.0)*isign*MP_PI/istep;
wtemp = MICROPY_FLOAT_C_FUN(sin)(MICROPY_FLOAT_CONST(0.5) * theta);
wpr = MICROPY_FLOAT_CONST(-2.0) * wtemp * wtemp;
wpi = MICROPY_FLOAT_C_FUN(sin)(theta);
wr = MICROPY_FLOAT_CONST(1.0);
wi = MICROPY_FLOAT_CONST(0.0);
for(m = 0; m < mmax; m++) {
for(size_t i = m; i < n; i += istep) {
j = i + mmax;
tempr = wr * data[2*j] - wi * data[2*j+1];
tempi = wr * data[2*j+1] + wi * data[2*j];
data[2*j] = data[2*i] - tempr;
data[2*j+1] = data[2*i+1] - tempi;
data[2*i] += tempr;
data[2*i+1] += tempi;
}
wtemp = wr;
wr = wr*wpr - wi*wpi + wr;
wi = wi*wpr + wtemp*wpi + wi;
}
mmax = istep;
}
}
/*
* The following function is a helper interface to the python side.
* It has been factored out from fft.c, so that the same argument parsing
* routine can be called from scipy.signal.spectrogram.
*/
mp_obj_t fft_fft_ifft_spectrogram(mp_obj_t data_in, uint8_t type) {
if(!mp_obj_is_type(data_in, &ulab_ndarray_type)) {
mp_raise_NotImplementedError(translate("FFT is defined for ndarrays only"));
}
ndarray_obj_t *in = MP_OBJ_TO_PTR(data_in);
#if ULAB_MAX_DIMS > 1
if(in->ndim != 1) {
mp_raise_TypeError(translate("FFT is implemented for linear arrays only"));
}
#endif
size_t len = in->len;
// Check if input is of length of power of 2
if((len & (len-1)) != 0) {
mp_raise_ValueError(translate("input array length must be power of 2"));
}
ndarray_obj_t *out = ndarray_new_linear_array(len, NDARRAY_COMPLEX);
mp_float_t *data = (mp_float_t *)out->array;
uint8_t *array = (uint8_t *)in->array;
if(in->dtype == NDARRAY_COMPLEX) {
uint8_t sz = 2 * sizeof(mp_float_t);
uint8_t *data_ = (uint8_t *)out->array;
for(size_t i = 0; i < len; i++) {
memcpy(data_, array, sz);
array += in->strides[ULAB_MAX_DIMS - 1];
}
} else {
mp_float_t (*func)(void *) = ndarray_get_float_function(in->dtype);
for(size_t i = 0; i < len; i++) {
// real part; the imaginary part is 0, no need to assign
*data = func(array);
data += 2;
array += in->strides[ULAB_MAX_DIMS - 1];
}
}
data -= 2 * len;
if((type == FFT_FFT) || (type == FFT_SPECTROGRAM)) {
fft_kernel_complex(data, len, 1);
if(type == FFT_SPECTROGRAM) {
ndarray_obj_t *spectrum = ndarray_new_linear_array(len, NDARRAY_FLOAT);
mp_float_t *sarray = (mp_float_t *)spectrum->array;
for(size_t i = 0; i < len; i++) {
*sarray++ = MICROPY_FLOAT_C_FUN(sqrt)(data[0] * data[0] + data[1] * data[1]);
data += 2;
}
m_del(mp_float_t, data, 2 * len);
return MP_OBJ_FROM_PTR(spectrum);
}
} else { // inverse transform
fft_kernel_complex(data, len, -1);
// TODO: numpy accepts the norm keyword argument
for(size_t i = 0; i < len; i++) {
*data++ /= len;
}
}
return MP_OBJ_FROM_PTR(out);
}
#else /* ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE */
void fft_kernel(mp_float_t *real, mp_float_t *imag, size_t n, int isign) {
size_t j, m, mmax, istep;
mp_float_t tempr, tempi;
@@ -77,12 +202,6 @@ void fft_kernel(mp_float_t *real, mp_float_t *imag, size_t n, int isign) {
}
}
/*
* The following function is a helper interface to the python side.
* It has been factored out from fft.c, so that the same argument parsing
* routine can be called from scipy.signal.spectrogram.
*/
mp_obj_t fft_fft_ifft_spectrogram(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im, uint8_t type) {
if(!mp_obj_is_type(arg_re, &ulab_ndarray_type)) {
mp_raise_NotImplementedError(translate("FFT is defined for ndarrays only"));
@@ -95,6 +214,7 @@ mp_obj_t fft_fft_ifft_spectrogram(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_i
ndarray_obj_t *re = MP_OBJ_TO_PTR(arg_re);
#if ULAB_MAX_DIMS > 1
if(re->ndim != 1) {
COMPLEX_DTYPE_NOT_IMPLEMENTED(re->dtype)
mp_raise_TypeError(translate("FFT is implemented for linear arrays only"));
}
#endif
@@ -122,6 +242,7 @@ mp_obj_t fft_fft_ifft_spectrogram(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_i
ndarray_obj_t *im = MP_OBJ_TO_PTR(arg_im);
#if ULAB_MAX_DIMS > 1
if(im->ndim != 1) {
COMPLEX_DTYPE_NOT_IMPLEMENTED(im->dtype)
mp_raise_TypeError(translate("FFT is implemented for linear arrays only"));
}
#endif
@@ -163,3 +284,4 @@ mp_obj_t fft_fft_ifft_spectrogram(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_i
return mp_obj_new_tuple(2, tuple);
}
}
#endif /* ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE */

5
python/port/mod/ulab/numpy/fft/fft_tools.h
View File

@@ -17,7 +17,12 @@ enum FFT_TYPE {
FFT_SPECTROGRAM,
};
#if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE
void fft_kernel(mp_float_t *, size_t , int );
mp_obj_t fft_fft_ifft_spectrogram(mp_obj_t , uint8_t );
#else
void fft_kernel(mp_float_t *, mp_float_t *, size_t , int );
mp_obj_t fft_fft_ifft_spectrogram(size_t , mp_obj_t , mp_obj_t , uint8_t );
#endif /* ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE */
#endif /* _FFT_TOOLS_ */

66
python/port/mod/ulab/numpy/filter.c
View File

@@ -21,6 +21,7 @@
#include "../ulab.h"
#include "../scipy/signal/signal.h"
#include "carray/carray_tools.h"
#include "filter.h"
#if ULAB_NUMPY_HAS_CONVOLVE
@@ -53,30 +54,77 @@ mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_a
}
int len = len_a + len_c - 1; // convolve mode "full"
ndarray_obj_t *out = ndarray_new_linear_array(len, NDARRAY_FLOAT);
mp_float_t *outptr = (mp_float_t *)out->array;
int32_t off = len_c - 1;
uint8_t dtype = NDARRAY_FLOAT;
#if ULAB_SUPPORTS_COMPLEX
if((a->dtype == NDARRAY_COMPLEX) || (c->dtype == NDARRAY_COMPLEX)) {
dtype = NDARRAY_COMPLEX;
}
#endif
ndarray_obj_t *ndarray = ndarray_new_linear_array(len, dtype);
mp_float_t *array = (mp_float_t *)ndarray->array;
uint8_t *aarray = (uint8_t *)a->array;
uint8_t *carray = (uint8_t *)c->array;
int32_t off = len_c - 1;
int32_t as = a->strides[ULAB_MAX_DIMS - 1] / a->itemsize;
int32_t cs = c->strides[ULAB_MAX_DIMS - 1] / c->itemsize;
for(int32_t k=-off; k < len-off; k++) {
mp_float_t accum = (mp_float_t)0.0;
#if ULAB_SUPPORTS_COMPLEX
if(dtype == NDARRAY_COMPLEX) {
mp_float_t a_real, a_imag;
mp_float_t c_real, c_imag = MICROPY_FLOAT_CONST(0.0);
for(int32_t k = -off; k < len-off; k++) {
mp_float_t accum_real = MICROPY_FLOAT_CONST(0.0);
mp_float_t accum_imag = MICROPY_FLOAT_CONST(0.0);
int32_t top_n = MIN(len_c, len_a - k);
int32_t bot_n = MAX(-k, 0);
for(int32_t n = bot_n; n < top_n; n++) {
int32_t idx_c = (len_c - n - 1) * cs;
int32_t idx_a = (n + k) * as;
if(a->dtype != NDARRAY_COMPLEX) {
a_real = ndarray_get_float_index(aarray, a->dtype, idx_a);
a_imag = MICROPY_FLOAT_CONST(0.0);
} else {
a_real = ndarray_get_float_index(aarray, NDARRAY_FLOAT, 2 * idx_a);
a_imag = ndarray_get_float_index(aarray, NDARRAY_FLOAT, 2 * idx_a + 1);
}
if(c->dtype != NDARRAY_COMPLEX) {
c_real = ndarray_get_float_index(carray, c->dtype, idx_c);
c_imag = MICROPY_FLOAT_CONST(0.0);
} else {
c_real = ndarray_get_float_index(carray, NDARRAY_FLOAT, 2 * idx_c);
c_imag = ndarray_get_float_index(carray, NDARRAY_FLOAT, 2 * idx_c + 1);
}
accum_real += a_real * c_real - a_imag * c_imag;
accum_imag += a_real * c_imag + a_imag * c_real;
}
*array++ = accum_real;
*array++ = accum_imag;
}
return MP_OBJ_FROM_PTR(ndarray);
}
#endif
for(int32_t k = -off; k < len-off; k++) {
mp_float_t accum = MICROPY_FLOAT_CONST(0.0);
int32_t top_n = MIN(len_c, len_a - k);
int32_t bot_n = MAX(-k, 0);
for(int32_t n=bot_n; n < top_n; n++) {
for(int32_t n = bot_n; n < top_n; n++) {
int32_t idx_c = (len_c - n - 1) * cs;
int32_t idx_a = (n + k) * as;
mp_float_t ai = ndarray_get_float_index(aarray, a->dtype, idx_a);
mp_float_t ci = ndarray_get_float_index(carray, c->dtype, idx_c);
accum += ai * ci;
}
*outptr++ = accum;
*array++ = accum;
}
return out;
return MP_OBJ_FROM_PTR(ndarray);
}
MP_DEFINE_CONST_FUN_OBJ_KW(filter_convolve_obj, 2, filter_convolve);

817
python/port/mod/ulab/numpy/io/io.c Normal file
View File

@@ -0,0 +1,817 @@
/*
* This file is part of the micropython-ulab project,
*
* https://github.com/v923z/micropython-ulab
*
* The MIT License (MIT)
*
* Copyright (c) 2022 Zoltán Vörös
*/
#include <math.h>
#include <string.h>
#include "py/builtin.h"
#include "py/formatfloat.h"
#include "py/obj.h"
#include "py/parsenum.h"
#include "py/runtime.h"
#include "py/stream.h"
#include "extmod/vfs.h"
#include "../../ndarray.h"
#include "../../ulab_tools.h"
#include "io.h"
#define ULAB_IO_BUFFER_SIZE 128
#define ULAB_IO_CLIPBOARD_SIZE 32
#define ULAB_IO_MAX_ROWS 65535
#define ULAB_IO_NULL_ENDIAN 0
#define ULAB_IO_LITTLE_ENDIAN 1
#define ULAB_IO_BIG_ENDIAN 2
#if ULAB_NUMPY_HAS_LOAD
static void io_read_(mp_obj_t stream, const mp_stream_p_t *stream_p, char *buffer, char *string, uint16_t len, int *error) {
size_t read = stream_p->read(stream, buffer, len, error);
bool fail = false;
if(read == len) {
if(string != NULL) {
if(memcmp(buffer, string, len) != 0) {
fail = true;
}
}
} else {
fail = true;
}
if(fail) {
stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, error);
mp_raise_msg(&mp_type_RuntimeError, translate("corrupted file"));
}
}
static mp_obj_t io_load(mp_obj_t file) {
if(!mp_obj_is_str(file)) {
mp_raise_TypeError(translate("wrong input type"));
}
int error;
char *buffer = m_new(char, ULAB_IO_BUFFER_SIZE);
// test for endianness
uint16_t x = 1;
int8_t native_endianness = (x >> 8) == 1 ? ULAB_IO_BIG_ENDIAN : ULAB_IO_LITTLE_ENDIAN;
mp_obj_t open_args[2] = {
file,
MP_OBJ_NEW_QSTR(MP_QSTR_rb)
};
mp_obj_t stream = mp_builtin_open_obj.fun.kw(2, open_args, (mp_map_t *)&mp_const_empty_map);
const mp_stream_p_t *stream_p = mp_get_stream(stream);
// read header
// magic string
io_read_(stream, stream_p, buffer, "\x93NUMPY", 6, &error);
// simply discard the version number
io_read_(stream, stream_p, buffer, NULL, 2, &error);
// header length, represented as a little endian uint16 (0x76, 0x00)
io_read_(stream, stream_p, buffer, NULL, 2, &error);
uint16_t header_length = buffer[1];
header_length <<= 8;
header_length += buffer[0];
// beginning of the dictionary describing the array
io_read_(stream, stream_p, buffer, "{'descr': '", 11, &error);
uint8_t dtype;
io_read_(stream, stream_p, buffer, NULL, 1, &error);
uint8_t endianness = ULAB_IO_NULL_ENDIAN;
if(*buffer == '<') {
endianness = ULAB_IO_LITTLE_ENDIAN;
} else if(*buffer == '>') {
endianness = ULAB_IO_BIG_ENDIAN;
}
io_read_(stream, stream_p, buffer, NULL, 2, &error);
if(memcmp(buffer, "u1", 2) == 0) {
dtype = NDARRAY_UINT8;
} else if(memcmp(buffer, "i1", 2) == 0) {
dtype = NDARRAY_INT8;
} else if(memcmp(buffer, "u2", 2) == 0) {
dtype = NDARRAY_UINT16;
} else if(memcmp(buffer, "i2", 2) == 0) {
dtype = NDARRAY_INT16;
}
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
else if(memcmp(buffer, "f4", 2) == 0) {
dtype = NDARRAY_FLOAT;
}
#else
else if(memcmp(buffer, "f8", 2) == 0) {
dtype = NDARRAY_FLOAT;
}
#endif
#if ULAB_SUPPORTS_COMPLEX
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
else if(memcmp(buffer, "c8", 2) == 0) {
dtype = NDARRAY_COMPLEX;
}
#else
else if(memcmp(buffer, "c16", 3) == 0) {
dtype = NDARRAY_COMPLEX;
}
#endif
#endif /* ULAB_SUPPORT_COPMLEX */
else {
stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
mp_raise_TypeError(translate("wrong dtype"));
}
io_read_(stream, stream_p, buffer, "', 'fortran_order': False, 'shape': (", 37, &error);
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
uint16_t bytes_to_read = MIN(ULAB_IO_BUFFER_SIZE, header_length - 51);
// bytes_to_read is 128 at most. This should be enough to contain a
// maximum of 4 size_t numbers plus the delimiters
io_read_(stream, stream_p, buffer, NULL, bytes_to_read, &error);
char *needle = buffer;
uint8_t ndim = 0;
// find out the number of dimensions by counting the commas in the string
while(1) {
if(*needle == ',') {
ndim++;
if(needle[1] == ')') {
break;
}
} else if((*needle == ')') && (ndim > 0)) {
ndim++;
break;
}
needle++;
}
needle = buffer;
for(uint8_t i = 0; i < ndim; i++) {
size_t number = 0;
// trivial number parsing here
while(1) {
if((*needle == ' ') || (*needle == '\t')) {
needle++;
}
if((*needle > 47) && (*needle < 58)) {
number = number * 10 + (*needle - 48);
} else if((*needle == ',') || (*needle == ')')) {
break;
}
else {
stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
mp_raise_msg(&mp_type_RuntimeError, translate("corrupted file"));
}
needle++;
}
needle++;
shape[ULAB_MAX_DIMS - ndim + i] = number;
}
// strip the rest of the header
if((bytes_to_read + 51) < header_length) {
io_read_(stream, stream_p, buffer, NULL, header_length - (bytes_to_read + 51), &error);
}
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(ndim, shape, dtype);
char *array = (char *)ndarray->array;
size_t read = stream_p->read(stream, array, ndarray->len * ndarray->itemsize, &error);
if(read != ndarray->len * ndarray->itemsize) {
stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
mp_raise_msg(&mp_type_RuntimeError, translate("corrupted file"));
}
stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
m_del(char, buffer, ULAB_IO_BUFFER_SIZE);
// swap the bytes, if necessary
if((native_endianness != endianness) && (dtype != NDARRAY_UINT8) && (dtype != NDARRAY_INT8)) {
uint8_t sz = ndarray->itemsize;
char *tmpbuff = NULL;
#if ULAB_SUPPORTS_COMPLEX
if(dtype == NDARRAY_COMPLEX) {
// work with the floating point real and imaginary parts
sz /= 2;
tmpbuff = m_new(char, sz);
for(size_t i = 0; i < ndarray->len; i++) {
for(uint8_t k = 0; k < 2; k++) {
tmpbuff += sz;
for(uint8_t j = 0; j < sz; j++) {
memcpy(--tmpbuff, array++, 1);
}
memcpy(array-sz, tmpbuff, sz);
}
}
} else {
#endif
tmpbuff = m_new(char, sz);
for(size_t i = 0; i < ndarray->len; i++) {
tmpbuff += sz;
for(uint8_t j = 0; j < sz; j++) {
memcpy(--tmpbuff, array++, 1);
}
memcpy(array-sz, tmpbuff, sz);
}
#if ULAB_SUPPORTS_COMPLEX
}
#endif
m_del(char, tmpbuff, sz);
}
m_del(size_t, shape, ULAB_MAX_DIMS);
return MP_OBJ_FROM_PTR(ndarray);
}
MP_DEFINE_CONST_FUN_OBJ_1(io_load_obj, io_load);
#endif /* ULAB_NUMPY_HAS_LOAD */
#if ULAB_NUMPY_HAS_LOADTXT
static void io_assign_value(const char *clipboard, uint8_t len, ndarray_obj_t *ndarray, size_t *idx, uint8_t dtype) {
mp_obj_t value = mp_parse_num_decimal(clipboard, len, false, false, NULL);
if(dtype != NDARRAY_FLOAT) {
mp_float_t _value = mp_obj_get_float(value);
value = mp_obj_new_int((int32_t)MICROPY_FLOAT_C_FUN(round)(_value));
}
ndarray_set_value(dtype, ndarray->array, (*idx)++, value);
}
static mp_obj_t io_loadtxt(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_const_none } },
{ MP_QSTR_delimiter, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_comments, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_max_rows, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = -1 } },
{ MP_QSTR_usecols, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = NDARRAY_FLOAT } },
{ MP_QSTR_skiprows, MP_ARG_KW_ONLY | MP_ARG_INT, { .u_int = 0 } },
};
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 open_args[2] = {
args[0].u_obj,
MP_OBJ_NEW_QSTR(MP_QSTR_r)
};
mp_obj_t stream = mp_builtin_open_obj.fun.kw(2, open_args, (mp_map_t *)&mp_const_empty_map);
const mp_stream_p_t *stream_p = mp_get_stream(stream);
char *buffer = m_new(char, ULAB_IO_BUFFER_SIZE);
int error;
char delimiter = ' ';
if(args[1].u_obj != mp_const_none) {
size_t _len;
char *_delimiter = m_new(char, 8);
_delimiter = (char *)mp_obj_str_get_data(args[1].u_obj, &_len);
delimiter = _delimiter[0];
}
char comment_char = '#';
if(args[2].u_obj != mp_const_none) {
size_t _len;
char *_comment_char = m_new(char, 8);
_comment_char = (char *)mp_obj_str_get_data(args[2].u_obj, &_len);
comment_char = _comment_char[0];
}
uint16_t skiprows = args[6].u_int;
uint16_t max_rows = ULAB_IO_MAX_ROWS;
if((args[3].u_int > 0) && (args[3].u_int < ULAB_IO_MAX_ROWS)) {
max_rows = args[3].u_int + skiprows;
}
uint16_t *cols = NULL;
uint8_t used_columns = 0;
if(args[4].u_obj != mp_const_none) {
if(mp_obj_is_int(args[4].u_obj)) {
used_columns = 1;
cols = m_new(uint16_t, used_columns);
cols[0] = (uint16_t)mp_obj_get_int(args[4].u_obj);
} else {
#if ULAB_MAX_DIMS == 1
mp_raise_ValueError(translate("usecols keyword must be specified"));
#else
// assume that the argument is an iterable
used_columns = (uint16_t)mp_obj_get_int(mp_obj_len(args[4].u_obj));
cols = m_new(uint16_t, used_columns);
mp_obj_iter_buf_t iter_buf;
mp_obj_t item, iterable = mp_getiter(args[4].u_obj, &iter_buf);
while((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
*cols++ = (uint16_t)mp_obj_get_int(item);
}
cols -= used_columns;
#endif
}
}
uint8_t dtype = args[5].u_int;
// count the columns and rows
// we actually count only the rows and the items, and assume that
// the number of columns can be gotten by means of a simple division,
// i.e., that each row has the same number of columns
char *offset;
uint16_t rows = 0, items = 0, all_rows = 0;
uint8_t read;
uint8_t len = 0;
do {
read = (uint8_t)stream_p->read(stream, buffer, ULAB_IO_BUFFER_SIZE - 1, &error);
buffer[read] = '\0';
offset = buffer;
while(*offset != '\0') {
if(*offset == comment_char) {
// clear the line till the end, or the buffer's end
while((*offset != '\0')) {
offset++;
if(*offset == '\n') {
offset++;
all_rows++;
break;
}
}
}
// catch whitespaces here: if these are not on a comment line, then they delimit a number
if(*offset == '\n') {
all_rows++;
if(all_rows > skiprows) {
rows++;
items++;
len = 0;
}
if(all_rows == max_rows) {
break;
}
}
if((*offset == ' ') || (*offset == '\t') || (*offset == '\v') ||
(*offset == '\f') || (*offset == '\r') || (*offset == delimiter)) {
offset++;
while((*offset == ' ') || (*offset == '\t') || (*offset == '\v') || (*offset == '\f') || (*offset == '\r')) {
offset++;
}
if(len > 0) {
if(all_rows >= skiprows) {
items++;
}
len = 0;
}
} else {
offset++;
len++;
}
}
} while((read > 0) && (all_rows < max_rows));
if(rows == 0) {
mp_raise_ValueError(translate("empty file"));
}
uint16_t columns = items / rows;
if(columns < used_columns) {
mp_raise_ValueError(translate("usecols is too high"));
}
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
#if ULAB_MAX_DIMS == 1
shape[0] = rows;
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(1, shape, dtype);
#else
if(args[4].u_obj == mp_const_none) {
shape[ULAB_MAX_DIMS - 1] = columns;
} else {
shape[ULAB_MAX_DIMS - 1] = used_columns;
}
shape[ULAB_MAX_DIMS - 2] = rows;
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(2, shape, dtype);
#endif
struct mp_stream_seek_t seek_s;
seek_s.offset = 0;
seek_s.whence = MP_SEEK_SET;
stream_p->ioctl(stream, MP_STREAM_SEEK, (mp_uint_t)(uintptr_t)&seek_s, &error);
char *clipboard = m_new(char, ULAB_IO_CLIPBOARD_SIZE);
char *clipboard_origin = clipboard;
rows = 0;
columns = 0;
len = 0;
size_t idx = 0;
do {
read = stream_p->read(stream, buffer, ULAB_IO_BUFFER_SIZE - 1, &error);
buffer[read] = '\0';
offset = buffer;
while(*offset != '\0') {
if(*offset == comment_char) {
// clear the line till the end, or the buffer's end
while((*offset != '\0')) {
offset++;
if(*offset == '\n') {
rows++;
offset++;
break;
}
}
}
if(rows == max_rows) {
break;
}
if((*offset == ' ') || (*offset == '\t') || (*offset == '\v') ||
(*offset == '\f') || (*offset == '\r') || (*offset == '\n') || (*offset == delimiter)) {
offset++;
while((*offset == ' ') || (*offset == '\t') || (*offset == '\v') ||
(*offset == '\f') || (*offset == '\r') || (*offset == '\n')) {
offset++;
}
if(len > 0) {
clipboard = clipboard_origin;
if(rows >= skiprows) {
#if ULAB_MAX_DIMS == 1
if(columns == cols[0]) {
io_assign_value(clipboard, len, ndarray, &idx, dtype);
}
#else
if(args[4].u_obj == mp_const_none) {
io_assign_value(clipboard, len, ndarray, &idx, dtype);
} else {
for(uint8_t c = 0; c < used_columns; c++) {
if(columns == cols[c]) {
io_assign_value(clipboard, len, ndarray, &idx, dtype);
break;
}
}
}
#endif
}
columns++;
len = 0;
if(offset[-1] == '\n') {
columns = 0;
rows++;
}
}
} else {
*clipboard++ = *offset++;
len++;
}
}
} while((read > 0) && (rows < max_rows));
stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
m_del(size_t, shape, ULAB_MAX_DIMS);
m_del(char, buffer, ULAB_IO_BUFFER_SIZE);
m_del(char, clipboard, ULAB_IO_CLIPBOARD_SIZE);
m_del(uint16_t, cols, used_columns);
return MP_OBJ_FROM_PTR(ndarray);
}
MP_DEFINE_CONST_FUN_OBJ_KW(io_loadtxt_obj, 1, io_loadtxt);
#endif /* ULAB_NUMPY_HAS_LOADTXT */
#if ULAB_NUMPY_HAS_SAVE
static uint8_t io_sprintf(char *buffer, const char *comma, size_t x) {
uint8_t offset = 1;
char *buf = buffer;
// our own minimal implementation of sprintf for size_t types
// this is required on systems, where sprintf is not available
// find out, how many characters are required
// we could call log10 here...
for(size_t i = 10; i < 100000000; i *= 10) {
if(x < i) {
break;
}
buf++;
}
while(x > 0) {
uint8_t rem = x % 10;
*buf-- = '0' + rem;
x /= 10;
offset++;
}
buf += offset;
while(*comma != '\0') {
*buf++ = *comma++;
offset++;
}
return offset - 1;
}
static mp_obj_t io_save(mp_obj_t file, mp_obj_t ndarray_) {
if(!mp_obj_is_str(file) || !mp_obj_is_type(ndarray_, &ulab_ndarray_type)) {
mp_raise_TypeError(translate("wrong input type"));
}
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(ndarray_);
int error;
char *buffer = m_new(char, ULAB_IO_BUFFER_SIZE);
uint8_t offset = 0;
// test for endianness
uint16_t x = 1;
int8_t native_endianness = (x >> 8) == 1 ? '>' : '<';
mp_obj_t open_args[2] = {
file,
MP_OBJ_NEW_QSTR(MP_QSTR_wb)
};
mp_obj_t stream = mp_builtin_open_obj.fun.kw(2, open_args, (mp_map_t *)&mp_const_empty_map);
const mp_stream_p_t *stream_p = mp_get_stream(stream);
// write header;
// magic string + header length, which is always 128 - 10 = 118, represented as a little endian uint16 (0x76, 0x00)
// + beginning of the dictionary describing the array
memcpy(buffer, "\x93NUMPY\x01\x00\x76\x00{'descr': '", 21);
offset += 21;
buffer[offset] = native_endianness;
if((ndarray->dtype == NDARRAY_UINT8) || (ndarray->dtype == NDARRAY_INT8)) {
// for single-byte data, the endianness doesn't matter
buffer[offset] = '|';
}
offset++;
switch(ndarray->dtype) {
case NDARRAY_UINT8:
memcpy(buffer+offset, "u1", 2);
break;
case NDARRAY_INT8:
memcpy(buffer+offset, "i1", 2);
break;
case NDARRAY_UINT16:
memcpy(buffer+offset, "u2", 2);
break;
case NDARRAY_INT16:
memcpy(buffer+offset, "i2", 2);
break;
case NDARRAY_FLOAT:
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
memcpy(buffer+offset, "f4", 2);
#else
memcpy(buffer+offset, "f8", 2);
#endif
break;
#if ULAB_SUPPORTS_COMPLEX
case NDARRAY_COMPLEX:
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
memcpy(buffer+offset, "c8", 2);
#else
memcpy(buffer+offset, "c16", 3);
offset++;
#endif
break;
#endif
}
offset += 2;
memcpy(buffer+offset, "', 'fortran_order': False, 'shape': (", 37);
offset += 37;
if(ndarray->ndim == 1) {
offset += io_sprintf(buffer+offset, ",\0", ndarray->shape[ULAB_MAX_DIMS - 1]);
} else {
for(uint8_t i = ndarray->ndim; i > 1; i--) {
offset += io_sprintf(buffer+offset, ", \0", ndarray->shape[ULAB_MAX_DIMS - i]);
}
offset += io_sprintf(buffer+offset, "\0", ndarray->shape[ULAB_MAX_DIMS - 1]);
}
memcpy(buffer+offset, "), }", 4);
offset += 4;
// pad with space till the very end
memset(buffer+offset, 32, ULAB_IO_BUFFER_SIZE - offset - 1);
buffer[ULAB_IO_BUFFER_SIZE - 1] = '\n';
stream_p->write(stream, buffer, ULAB_IO_BUFFER_SIZE, &error);
// write the array data
uint8_t sz = ndarray->itemsize;
offset = 0;
uint8_t *array = (uint8_t *)ndarray->array;
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
memcpy(buffer+offset, array, sz);
offset += sz;
if(offset == ULAB_IO_BUFFER_SIZE) {
stream_p->write(stream, buffer, offset, &error);
offset = 0;
}
array += ndarray->strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < ndarray->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
array -= ndarray->strides[ULAB_MAX_DIMS - 1] * ndarray->shape[ULAB_MAX_DIMS-1];
array += ndarray->strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < ndarray->shape[ULAB_MAX_DIMS - 2]);
#endif
#if ULAB_MAX_DIMS > 2
array -= ndarray->strides[ULAB_MAX_DIMS - 2] * ndarray->shape[ULAB_MAX_DIMS-2];
array += ndarray->strides[ULAB_MAX_DIMS - 3];
j++;
} while(j < ndarray->shape[ULAB_MAX_DIMS - 3]);
#endif
#if ULAB_MAX_DIMS > 3
array -= ndarray->strides[ULAB_MAX_DIMS - 3] * ndarray->shape[ULAB_MAX_DIMS-3];
array += ndarray->strides[ULAB_MAX_DIMS - 4];
i++;
} while(i < ndarray->shape[ULAB_MAX_DIMS - 4]);
#endif
stream_p->write(stream, buffer, offset, &error);
stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
m_del(char, buffer, ULAB_IO_BUFFER_SIZE);
return mp_const_none;
}
MP_DEFINE_CONST_FUN_OBJ_2(io_save_obj, io_save);
#endif /* ULAB_NUMPY_HAS_SAVE */
#if ULAB_NUMPY_HAS_SAVETXT
static int8_t io_format_float(ndarray_obj_t *ndarray, mp_float_t (*func)(void *), uint8_t *array, char *buffer, char *delimiter) {
// own implementation of float formatting for platforms that don't have sprintf
int8_t offset = 0;
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
#if MICROPY_OBJ_REPR == MICROPY_OBJ_REPR_C
const int precision = 6;
#else
const int precision = 7;
#endif
#else
const int precision = 16;
#endif
#if ULAB_SUPPORTS_COMPLEX
if(ndarray->dtype == NDARRAY_COMPLEX) {
mp_float_t real = func(array);
mp_float_t imag = func(array + ndarray->itemsize / 2);
offset = mp_format_float(real, buffer, ULAB_IO_BUFFER_SIZE, 'f', precision, 'j');
if(imag >= MICROPY_FLOAT_CONST(0.0)) {
buffer[offset++] = '+';
} else {
buffer[offset++] = '-';
}
offset += mp_format_float(-imag, &buffer[offset], ULAB_IO_BUFFER_SIZE, 'f', precision, 'j');
}
#endif
offset = (uint8_t)mp_format_float(func(array), buffer, ULAB_IO_BUFFER_SIZE, 'f', precision, '\0');
#if ULAB_SUPPORTS_COMPLEX
if(ndarray->dtype != NDARRAY_COMPLEX) {
// complexes end with a 'j', floats with a '\0', so we have to wind back by one character
offset--;
}
#endif
while(*delimiter != '\0') {
buffer[offset++] = *delimiter++;
}
return offset;
}
static mp_obj_t io_savetxt(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_const_none } },
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_delimiter, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_header, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_footer, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_comments, 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(!mp_obj_is_str(args[0].u_obj) || !mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) {
mp_raise_TypeError(translate("wrong input type"));
}
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[1].u_obj);
#if ULAB_MAX_DIMS > 2
if(ndarray->ndim > 2) {
mp_raise_ValueError(translate("array has too many dimensions"));
}
#endif
mp_obj_t open_args[2] = {
args[0].u_obj,
MP_OBJ_NEW_QSTR(MP_QSTR_w)
};
mp_obj_t stream = mp_builtin_open_obj.fun.kw(2, open_args, (mp_map_t *)&mp_const_empty_map);
const mp_stream_p_t *stream_p = mp_get_stream(stream);
char *buffer = m_new(char, ULAB_IO_BUFFER_SIZE);
int error;
if(mp_obj_is_str(args[3].u_obj)) {
size_t _len;
if(mp_obj_is_str(args[5].u_obj)) {
const char *comments = mp_obj_str_get_data(args[5].u_obj, &_len);
stream_p->write(stream, comments, _len, &error);
} else {
stream_p->write(stream, "# ", 2, &error);
}
const char *header = mp_obj_str_get_data(args[3].u_obj, &_len);
stream_p->write(stream, header, _len, &error);
stream_p->write(stream, "\n", 1, &error);
}
uint8_t *array = (uint8_t *)ndarray->array;
mp_float_t (*func)(void *) = ndarray_get_float_function(ndarray->dtype);
char *delimiter = m_new(char, 8);
if(ndarray->ndim == 1) {
delimiter[0] = '\n';
delimiter[1] = '\0';
} else if(args[2].u_obj == mp_const_none) {
delimiter[0] = ' ';
delimiter[1] = '\0';
} else {
size_t delimiter_len;
delimiter = (char *)mp_obj_str_get_data(args[2].u_obj, &delimiter_len);
}
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
int8_t chars = io_format_float(ndarray, func, array, buffer, l == ndarray->shape[ULAB_MAX_DIMS - 1] - 1 ? "\n" : delimiter);
if(chars > 0) {
stream_p->write(stream, buffer, chars, &error);
}
array += ndarray->strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < ndarray->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
array -= ndarray->strides[ULAB_MAX_DIMS - 1] * ndarray->shape[ULAB_MAX_DIMS-1];
array += ndarray->strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < ndarray->shape[ULAB_MAX_DIMS - 2]);
#endif
if(mp_obj_is_str(args[4].u_obj)) {
size_t _len;
if(mp_obj_is_str(args[5].u_obj)) {
const char *comments = mp_obj_str_get_data(args[5].u_obj, &_len);
stream_p->write(stream, comments, _len, &error);
} else {
stream_p->write(stream, "# ", 2, &error);
}
const char *footer = mp_obj_str_get_data(args[4].u_obj, &_len);
stream_p->write(stream, footer, _len, &error);
stream_p->write(stream, "\n", 1, &error);
}
stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
return mp_const_none;
}
MP_DEFINE_CONST_FUN_OBJ_KW(io_savetxt_obj, 2, io_savetxt);
#endif /* ULAB_NUMPY_HAS_SAVETXT */

19
python/port/mod/ulab/numpy/io/io.h Normal file
View File

@@ -0,0 +1,19 @@
/*
* This file is part of the micropython-ulab project,
*
* https://github.com/v923z/micropython-ulab
*
* The MIT License (MIT)
*
* Copyright (c) 2022 Zoltán Vörös
*/
#ifndef _ULAB_IO_
#define _ULAB_IO_
MP_DECLARE_CONST_FUN_OBJ_1(io_load_obj);
MP_DECLARE_CONST_FUN_OBJ_KW(io_loadtxt_obj);
MP_DECLARE_CONST_FUN_OBJ_2(io_save_obj);
MP_DECLARE_CONST_FUN_OBJ_KW(io_savetxt_obj);
#endif

14
python/port/mod/ulab/numpy/linalg/linalg.c
View File

@@ -22,6 +22,7 @@
#include "../../ulab.h"
#include "../../ulab_tools.h"
#include "../carray/carray_tools.h"
#include "linalg.h"
#if ULAB_NUMPY_HAS_LINALG_MODULE
@@ -44,6 +45,7 @@
static mp_obj_t linalg_cholesky(mp_obj_t oin) {
ndarray_obj_t *ndarray = tools_object_is_square(oin);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
ndarray_obj_t *L = ndarray_new_dense_ndarray(2, ndarray_shape_vector(0, 0, ndarray->shape[ULAB_MAX_DIMS - 1], ndarray->shape[ULAB_MAX_DIMS - 1]), NDARRAY_FLOAT);
mp_float_t *Larray = (mp_float_t *)L->array;
@@ -110,6 +112,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(linalg_cholesky_obj, linalg_cholesky);
static mp_obj_t linalg_det(mp_obj_t oin) {
ndarray_obj_t *ndarray = tools_object_is_square(oin);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
uint8_t *array = (uint8_t *)ndarray->array;
size_t N = ndarray->shape[ULAB_MAX_DIMS - 1];
mp_float_t *tmp = m_new(mp_float_t, N * N);
@@ -182,6 +185,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(linalg_det_obj, linalg_det);
static mp_obj_t linalg_eig(mp_obj_t oin) {
ndarray_obj_t *in = tools_object_is_square(oin);
COMPLEX_DTYPE_NOT_IMPLEMENTED(in->dtype)
uint8_t *iarray = (uint8_t *)in->array;
size_t S = in->shape[ULAB_MAX_DIMS - 1];
mp_float_t *array = m_new(mp_float_t, S*S);
@@ -243,6 +247,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(linalg_eig_obj, linalg_eig);
//|
static mp_obj_t linalg_inv(mp_obj_t o_in) {
ndarray_obj_t *ndarray = tools_object_is_square(o_in);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
uint8_t *array = (uint8_t *)ndarray->array;
size_t N = ndarray->shape[ULAB_MAX_DIMS - 1];
ndarray_obj_t *inverted = ndarray_new_dense_ndarray(2, ndarray_shape_vector(0, 0, N, N), NDARRAY_FLOAT);
@@ -305,6 +310,7 @@ static mp_obj_t linalg_norm(size_t n_args, const mp_obj_t *pos_args, mp_map_t *k
return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(dot * (count - 1)));
} else if(mp_obj_is_type(x, &ulab_ndarray_type)) {
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(x);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
uint8_t *array = (uint8_t *)ndarray->array;
// always get a float, so that we don't have to resolve the dtype later
mp_float_t (*func)(void *) = ndarray_get_float_function(ndarray->dtype);
@@ -429,22 +435,22 @@ static mp_obj_t linalg_qr(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_
// [[c s],
// [s -c]]
if(MICROPY_FLOAT_C_FUN(fabs)(rarray[i * n + j]) < LINALG_EPSILON) { // r[i, j]
c = (rarray[(i - 1) * n + j] >= 0.0) ? 1.0 : -1.0; // r[i-1, j]
c = (rarray[(i - 1) * n + j] >= MICROPY_FLOAT_CONST(0.0)) ? MICROPY_FLOAT_CONST(1.0) : MICROPY_FLOAT_CONST(-1.0); // r[i-1, j]
s = 0.0;
} else if(MICROPY_FLOAT_C_FUN(fabs)(rarray[(i - 1) * n + j]) < LINALG_EPSILON) { // r[i-1, j]
c = 0.0;
s = (rarray[i * n + j] >= 0.0) ? -1.0 : 1.0; // r[i, j]
s = (rarray[i * n + j] >= MICROPY_FLOAT_CONST(0.0)) ? MICROPY_FLOAT_CONST(-1.0) : MICROPY_FLOAT_CONST(1.0); // r[i, j]
} else {
mp_float_t t, u;
if(MICROPY_FLOAT_C_FUN(fabs)(rarray[(i - 1) * n + j]) > MICROPY_FLOAT_C_FUN(fabs)(rarray[i * n + j])) { // r[i-1, j], r[i, j]
t = rarray[i * n + j] / rarray[(i - 1) * n + j]; // r[i, j]/r[i-1, j]
u = MICROPY_FLOAT_C_FUN(sqrt)(1 + t * t);
c = -1.0 / u;
c = MICROPY_FLOAT_CONST(-1.0) / u;
s = c * t;
} else {
t = rarray[(i - 1) * n + j] / rarray[i * n + j]; // r[i-1, j]/r[i, j]
u = MICROPY_FLOAT_C_FUN(sqrt)(1 + t * t);
s = -1.0 / u;
s = MICROPY_FLOAT_CONST(-1.0) / u;
c = s * t;
}
}

17
python/port/mod/ulab/numpy/linalg/linalg_tools.c
View File

@@ -14,8 +14,8 @@
#include "linalg_tools.h"
/*
* The following function inverts a matrix, whose entries are given in the input array
/*
* The following function inverts a matrix, whose entries are given in the input array
* The function has no dependencies beyond micropython itself (for the definition of mp_float_t),
* and can be used independent of ulab.
*/
@@ -26,10 +26,9 @@ bool linalg_invert_matrix(mp_float_t *data, size_t N) {
// initially, this is the unit matrix: the contents of this matrix is what
// will be returned after all the transformations
mp_float_t *unit = m_new(mp_float_t, N*N);
mp_float_t *unit = m_new0(mp_float_t, N*N);
mp_float_t elem = 1.0;
// initialise the unit matrix
memset(unit, 0, sizeof(mp_float_t)*N*N);
for(size_t m=0; m < N; m++) {
memcpy(&unit[m * (N+1)], &elem, sizeof(mp_float_t));
}
@@ -78,9 +77,9 @@ bool linalg_invert_matrix(mp_float_t *data, size_t N) {
return true;
}
/*
* The following function calculates the eigenvalues and eigenvectors of a symmetric
* real matrix, whose entries are given in the input array.
/*
* The following function calculates the eigenvalues and eigenvectors of a symmetric
* real matrix, whose entries are given in the input array.
* The function has no dependencies beyond micropython itself (for the definition of mp_float_t),
* and can be used independent of ulab.
*/
@@ -166,6 +165,6 @@ size_t linalg_jacobi_rotations(mp_float_t *array, mp_float_t *eigvectors, size_t
eigvectors[m * S + N] = s * vm + c * vn;
}
} while(iterations > 0);
return iterations;
}

231
python/port/mod/ulab/numpy/numerical.c
View File

@@ -22,6 +22,7 @@
#include "../ulab.h"
#include "../ulab_tools.h"
#include "./carray/carray_tools.h"
#include "numerical.h"
enum NUMERICAL_FUNCTION_TYPE {
@@ -130,33 +131,71 @@ static mp_obj_t numerical_all_any(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {
size_t l = 0;
if(axis == mp_const_none) {
do {
mp_float_t value = func(array);
if((value != MICROPY_FLOAT_CONST(0.0)) & !anytype) {
// optype = NUMERICAL_ANY
return mp_const_true;
} else if((value == MICROPY_FLOAT_CONST(0.0)) & anytype) {
// optype == NUMERICAL_ALL
return mp_const_false;
#if ULAB_SUPPORTS_COMPLEX
if(ndarray->dtype == NDARRAY_COMPLEX) {
mp_float_t real = *((mp_float_t *)array);
mp_float_t imag = *((mp_float_t *)(array + sizeof(mp_float_t)));
if(((real != MICROPY_FLOAT_CONST(0.0)) | (imag != MICROPY_FLOAT_CONST(0.0))) & !anytype) {
// optype = NUMERICAL_ANY
return mp_const_true;
} else if(((real == MICROPY_FLOAT_CONST(0.0)) & (imag == MICROPY_FLOAT_CONST(0.0))) & anytype) {
// optype == NUMERICAL_ALL
return mp_const_false;
}
} else {
#endif
mp_float_t value = func(array);
if((value != MICROPY_FLOAT_CONST(0.0)) & !anytype) {
// optype = NUMERICAL_ANY
return mp_const_true;
} else if((value == MICROPY_FLOAT_CONST(0.0)) & anytype) {
// optype == NUMERICAL_ALL
return mp_const_false;
}
#if ULAB_SUPPORTS_COMPLEX
}
#endif
array += _shape_strides.strides[0];
l++;
} while(l < _shape_strides.shape[0]);
} else { // a scalar axis keyword was supplied
do {
mp_float_t value = func(array);
if((value != MICROPY_FLOAT_CONST(0.0)) & !anytype) {
// optype == NUMERICAL_ANY
*rarray = 1;
// since we are breaking out of the loop, move the pointer forward
array += _shape_strides.strides[0] * (_shape_strides.shape[0] - l);
break;
} else if((value == MICROPY_FLOAT_CONST(0.0)) & anytype) {
// optype == NUMERICAL_ALL
*rarray = 0;
// since we are breaking out of the loop, move the pointer forward
array += _shape_strides.strides[0] * (_shape_strides.shape[0] - l);
break;
#if ULAB_SUPPORTS_COMPLEX
if(ndarray->dtype == NDARRAY_COMPLEX) {
mp_float_t real = *((mp_float_t *)array);
mp_float_t imag = *((mp_float_t *)(array + sizeof(mp_float_t)));
if(((real != MICROPY_FLOAT_CONST(0.0)) | (imag != MICROPY_FLOAT_CONST(0.0))) & !anytype) {
// optype = NUMERICAL_ANY
*rarray = 1;
// since we are breaking out of the loop, move the pointer forward
array += _shape_strides.strides[0] * (_shape_strides.shape[0] - l);
break;
} else if(((real == MICROPY_FLOAT_CONST(0.0)) & (imag == MICROPY_FLOAT_CONST(0.0))) & anytype) {
// optype == NUMERICAL_ALL
*rarray = 0;
// since we are breaking out of the loop, move the pointer forward
array += _shape_strides.strides[0] * (_shape_strides.shape[0] - l);
break;
}
} else {
#endif
mp_float_t value = func(array);
if((value != MICROPY_FLOAT_CONST(0.0)) & !anytype) {
// optype == NUMERICAL_ANY
*rarray = 1;
// since we are breaking out of the loop, move the pointer forward
array += _shape_strides.strides[0] * (_shape_strides.shape[0] - l);
break;
} else if((value == MICROPY_FLOAT_CONST(0.0)) & anytype) {
// optype == NUMERICAL_ALL
*rarray = 0;
// since we are breaking out of the loop, move the pointer forward
array += _shape_strides.strides[0] * (_shape_strides.shape[0] - l);
break;
}
#if ULAB_SUPPORTS_COMPLEX
}
#endif
array += _shape_strides.strides[0];
l++;
} while(l < _shape_strides.shape[0]);
@@ -234,6 +273,7 @@ static mp_obj_t numerical_sum_mean_std_iterable(mp_obj_t oin, uint8_t optype, si
}
static mp_obj_t numerical_sum_mean_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype, size_t ddof) {
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
uint8_t *array = (uint8_t *)ndarray->array;
shape_strides _shape_strides = tools_reduce_axes(ndarray, axis);
@@ -244,7 +284,7 @@ static mp_obj_t numerical_sum_mean_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t
return mp_obj_new_float(MICROPY_FLOAT_CONST(0.0));
}
mp_float_t (*func)(void *) = ndarray_get_float_function(ndarray->dtype);
mp_float_t M =MICROPY_FLOAT_CONST(0.0);
mp_float_t M = MICROPY_FLOAT_CONST(0.0);
mp_float_t m = MICROPY_FLOAT_CONST(0.0);
mp_float_t S = MICROPY_FLOAT_CONST(0.0);
mp_float_t s = MICROPY_FLOAT_CONST(0.0);
@@ -472,17 +512,12 @@ static mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t
}
}
} else {
int8_t ax = mp_obj_get_int(axis);
if(ax < 0) ax += ndarray->ndim;
if((ax < 0) || (ax > ndarray->ndim - 1)) {
mp_raise_ValueError(translate("axis is out of bounds"));
}
int8_t ax = tools_get_axis(axis, ndarray->ndim);
uint8_t *array = (uint8_t *)ndarray->array;
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
memset(shape, 0, sizeof(size_t)*ULAB_MAX_DIMS);
int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
memset(strides, 0, sizeof(uint32_t)*ULAB_MAX_DIMS);
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
int32_t *strides = m_new0(int32_t, ULAB_MAX_DIMS);
numerical_reduce_axes(ndarray, ax, shape, strides);
uint8_t index = ULAB_MAX_DIMS - ndarray->ndim + ax;
@@ -507,6 +542,9 @@ static mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t
} else {
RUN_ARGMIN(ndarray, mp_float_t, array, results, rarray, shape, strides, index, optype);
}
m_del(int32_t, strides, ULAB_MAX_DIMS);
if(results->len == 1) {
return mp_binary_get_val_array(results->dtype, results->array, 0);
}
@@ -555,9 +593,11 @@ static mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_m
case NUMERICAL_MAX:
case NUMERICAL_ARGMIN:
case NUMERICAL_ARGMAX:
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
return numerical_argmin_argmax_ndarray(ndarray, axis, optype);
case NUMERICAL_SUM:
case NUMERICAL_MEAN:
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
return numerical_sum_mean_std_ndarray(ndarray, axis, optype, 0);
default:
mp_raise_NotImplementedError(translate("operation is not implemented on ndarrays"));
@@ -580,6 +620,7 @@ static mp_obj_t numerical_sort_helper(mp_obj_t oin, mp_obj_t axis, uint8_t inpla
} else {
ndarray = ndarray_copy_view(MP_OBJ_TO_PTR(oin));
}
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
int8_t ax = 0;
if(axis == mp_const_none) {
@@ -594,30 +635,30 @@ static mp_obj_t numerical_sort_helper(mp_obj_t oin, mp_obj_t axis, uint8_t inpla
ndarray->ndim = 1;
#endif
} else {
ax = mp_obj_get_int(axis);
if(ax < 0) ax += ndarray->ndim;
if((ax < 0) || (ax > ndarray->ndim - 1)) {
mp_raise_ValueError(translate("index out of range"));
}
ax = tools_get_axis(axis, ndarray->ndim);
}
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
memset(shape, 0, sizeof(size_t)*ULAB_MAX_DIMS);
int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
memset(strides, 0, sizeof(uint32_t)*ULAB_MAX_DIMS);
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
int32_t *strides = m_new0(int32_t, ULAB_MAX_DIMS);
numerical_reduce_axes(ndarray, ax, shape, strides);
ax = ULAB_MAX_DIMS - ndarray->ndim + ax;
// we work with the typed array, so re-scale the stride
int32_t increment = ndarray->strides[ax] / ndarray->itemsize;
uint8_t *array = (uint8_t *)ndarray->array;
if((ndarray->dtype == NDARRAY_UINT8) || (ndarray->dtype == NDARRAY_INT8)) {
HEAPSORT(ndarray, uint8_t, array, shape, strides, ax, increment, ndarray->shape[ax]);
} else if((ndarray->dtype == NDARRAY_INT16) || (ndarray->dtype == NDARRAY_INT16)) {
HEAPSORT(ndarray, uint16_t, array, shape, strides, ax, increment, ndarray->shape[ax]);
} else {
HEAPSORT(ndarray, mp_float_t, array, shape, strides, ax, increment, ndarray->shape[ax]);
if(ndarray->shape[ax]) {
if((ndarray->dtype == NDARRAY_UINT8) || (ndarray->dtype == NDARRAY_INT8)) {
HEAPSORT(ndarray, uint8_t, array, shape, strides, ax, increment, ndarray->shape[ax]);
} else if((ndarray->dtype == NDARRAY_INT16) || (ndarray->dtype == NDARRAY_INT16)) {
HEAPSORT(ndarray, uint16_t, array, shape, strides, ax, increment, ndarray->shape[ax]);
} else {
HEAPSORT(ndarray, mp_float_t, array, shape, strides, ax, increment, ndarray->shape[ax]);
}
}
m_del(int32_t, strides, ULAB_MAX_DIMS);
if(inplace == 1) {
return mp_const_none;
} else {
@@ -682,6 +723,7 @@ mp_obj_t numerical_argsort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
}
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
if(args[1].u_obj == mp_const_none) {
// bail out, though dense arrays could still be sorted
mp_raise_NotImplementedError(translate("argsort is not implemented for flattened arrays"));
@@ -693,22 +735,17 @@ mp_obj_t numerical_argsort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
mp_raise_ValueError(translate("axis too long"));
}
}
int8_t ax = mp_obj_get_int(args[1].u_obj);
if(ax < 0) ax += ndarray->ndim;
if((ax < 0) || (ax > ndarray->ndim - 1)) {
mp_raise_ValueError(translate("index out of range"));
}
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
memset(shape, 0, sizeof(size_t)*ULAB_MAX_DIMS);
int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
memset(strides, 0, sizeof(uint32_t)*ULAB_MAX_DIMS);
int8_t ax = tools_get_axis(args[1].u_obj, ndarray->ndim);
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
int32_t *strides = m_new0(int32_t, ULAB_MAX_DIMS);
numerical_reduce_axes(ndarray, ax, shape, strides);
// We could return an NDARRAY_UINT8 array, if all lengths are shorter than 256
ndarray_obj_t *indices = ndarray_new_ndarray(ndarray->ndim, ndarray->shape, NULL, NDARRAY_UINT16);
int32_t *istrides = m_new(int32_t, ULAB_MAX_DIMS);
memset(istrides, 0, sizeof(uint32_t)*ULAB_MAX_DIMS);
int32_t *istrides = m_new0(int32_t, ULAB_MAX_DIMS);
numerical_reduce_axes(indices, ax, shape, istrides);
for(uint8_t i=0; i < ULAB_MAX_DIMS; i++) {
istrides[i] /= sizeof(uint16_t);
}
@@ -760,13 +797,20 @@ mp_obj_t numerical_argsort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
// reset the array
iarray = indices->array;
if((ndarray->dtype == NDARRAY_UINT8) || (ndarray->dtype == NDARRAY_INT8)) {
HEAP_ARGSORT(ndarray, uint8_t, array, shape, strides, ax, increment, ndarray->shape[ax], iarray, istrides, iincrement);
} else if((ndarray->dtype == NDARRAY_UINT16) || (ndarray->dtype == NDARRAY_INT16)) {
HEAP_ARGSORT(ndarray, uint16_t, array, shape, strides, ax, increment, ndarray->shape[ax], iarray, istrides, iincrement);
} else {
HEAP_ARGSORT(ndarray, mp_float_t, array, shape, strides, ax, increment, ndarray->shape[ax], iarray, istrides, iincrement);
if(ndarray->shape[ax]) {
if((ndarray->dtype == NDARRAY_UINT8) || (ndarray->dtype == NDARRAY_INT8)) {
HEAP_ARGSORT(ndarray, uint8_t, array, shape, strides, ax, increment, ndarray->shape[ax], iarray, istrides, iincrement);
} else if((ndarray->dtype == NDARRAY_UINT16) || (ndarray->dtype == NDARRAY_INT16)) {
HEAP_ARGSORT(ndarray, uint16_t, array, shape, strides, ax, increment, ndarray->shape[ax], iarray, istrides, iincrement);
} else {
HEAP_ARGSORT(ndarray, mp_float_t, array, shape, strides, ax, increment, ndarray->shape[ax], iarray, istrides, iincrement);
}
}
m_del(size_t, shape, ULAB_MAX_DIMS);
m_del(int32_t, strides, ULAB_MAX_DIMS);
m_del(int32_t, istrides, ULAB_MAX_DIMS);
return MP_OBJ_FROM_PTR(indices);
}
@@ -785,6 +829,8 @@ static mp_obj_t numerical_cross(mp_obj_t _a, mp_obj_t _b) {
}
ndarray_obj_t *a = MP_OBJ_TO_PTR(_a);
ndarray_obj_t *b = MP_OBJ_TO_PTR(_b);
COMPLEX_DTYPE_NOT_IMPLEMENTED(a->dtype)
COMPLEX_DTYPE_NOT_IMPLEMENTED(b->dtype)
if((a->ndim != 1) || (b->ndim != 1) || (a->len != b->len) || (a->len != 3)) {
mp_raise_ValueError(translate("cross is defined for 1D arrays of length 3"));
}
@@ -873,6 +919,7 @@ mp_obj_t numerical_diff(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
}
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
int8_t ax = args[2].u_int;
if(ax < 0) ax += ndarray->ndim;
@@ -891,13 +938,12 @@ mp_obj_t numerical_diff(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
int8_t *stencil = m_new(int8_t, N+1);
stencil[0] = 1;
for(uint8_t i=1; i < N+1; i++) {
for(uint8_t i = 1; i < N+1; i++) {
stencil[i] = -stencil[i-1]*(N-i+1)/i;
}
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
memset(shape, 0, sizeof(size_t)*ULAB_MAX_DIMS);
for(uint8_t i=0; i < ULAB_MAX_DIMS; i++) {
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
for(uint8_t i = 0; i < ULAB_MAX_DIMS; i++) {
shape[i] = ndarray->shape[i];
if(i == index) {
shape[i] -= N;
@@ -908,8 +954,7 @@ mp_obj_t numerical_diff(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
uint8_t *rarray = (uint8_t *)results->array;
memset(shape, 0, sizeof(size_t)*ULAB_MAX_DIMS);
int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
memset(strides, 0, sizeof(int32_t)*ULAB_MAX_DIMS);
int32_t *strides = m_new0(int32_t, ULAB_MAX_DIMS);
numerical_reduce_axes(ndarray, ax, shape, strides);
if(ndarray->dtype == NDARRAY_UINT8) {
@@ -956,17 +1001,14 @@ mp_obj_t numerical_flip(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
if(args[1].u_obj == mp_const_none) { // flip the flattened array
results = ndarray_new_linear_array(ndarray->len, ndarray->dtype);
ndarray_copy_array(ndarray, results);
ndarray_copy_array(ndarray, results, 0);
uint8_t *rarray = (uint8_t *)results->array;
rarray += (results->len - 1) * results->itemsize;
results->array = rarray;
results->strides[ULAB_MAX_DIMS - 1] = -results->strides[ULAB_MAX_DIMS - 1];
} else if(mp_obj_is_int(args[1].u_obj)){
int8_t ax = mp_obj_get_int(args[1].u_obj);
if(ax < 0) ax += ndarray->ndim;
if((ax < 0) || (ax > ndarray->ndim - 1)) {
mp_raise_ValueError(translate("index out of range"));
}
int8_t ax = tools_get_axis(args[1].u_obj, ndarray->ndim);
ax = ULAB_MAX_DIMS - ndarray->ndim + ax;
int32_t offset = (ndarray->shape[ax] - 1) * ndarray->strides[ax];
results = ndarray_new_view(ndarray, ndarray->ndim, ndarray->shape, ndarray->strides, offset);
@@ -1044,17 +1086,16 @@ mp_obj_t numerical_median(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_
}
return mp_obj_new_float(median);
} else {
int8_t ax = mp_obj_get_int(args[1].u_obj);
if(ax < 0) ax += ndarray->ndim;
// here we can save the exception, because if the axis is out of range,
// then numerical_sort_helper has already taken care of the issue
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
memset(shape, 0, sizeof(size_t)*ULAB_MAX_DIMS);
int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
memset(strides, 0, sizeof(uint32_t)*ULAB_MAX_DIMS);
int8_t ax = tools_get_axis(args[1].u_obj, ndarray->ndim);
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
int32_t *strides = m_new0(int32_t, ULAB_MAX_DIMS);
numerical_reduce_axes(ndarray, ax, shape, strides);
ax = ULAB_MAX_DIMS - ndarray->ndim + ax;
ndarray_obj_t *results = ndarray_new_dense_ndarray(ndarray->ndim-1, shape, NDARRAY_FLOAT);
m_del(size_t, shape, ULAB_MAX_DIMS);
mp_float_t *rarray = (mp_float_t *)results->array;
uint8_t *array = (uint8_t *)ndarray->array;
@@ -1200,21 +1241,14 @@ mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
} while(i < ndarray->shape[ULAB_MAX_DIMS - 4]);
#endif
} else if(mp_obj_is_int(args[2].u_obj)){
int8_t ax = mp_obj_get_int(args[2].u_obj);
if(ax < 0) ax += ndarray->ndim;
if((ax < 0) || (ax > ndarray->ndim - 1)) {
mp_raise_ValueError(translate("index out of range"));
}
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
memset(shape, 0, sizeof(size_t)*ULAB_MAX_DIMS);
int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
memset(strides, 0, sizeof(int32_t)*ULAB_MAX_DIMS);
int8_t ax = tools_get_axis(args[2].u_obj, ndarray->ndim);
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
int32_t *strides = m_new0(int32_t, ULAB_MAX_DIMS);
numerical_reduce_axes(ndarray, ax, shape, strides);
size_t *rshape = m_new(size_t, ULAB_MAX_DIMS);
memset(rshape, 0, sizeof(size_t)*ULAB_MAX_DIMS);
int32_t *rstrides = m_new(int32_t, ULAB_MAX_DIMS);
memset(rstrides, 0, sizeof(int32_t)*ULAB_MAX_DIMS);
size_t *rshape = m_new0(size_t, ULAB_MAX_DIMS);
int32_t *rstrides = m_new0(int32_t, ULAB_MAX_DIMS);
numerical_reduce_axes(results, ax, rshape, rstrides);
ax = ULAB_MAX_DIMS - ndarray->ndim + ax;
@@ -1275,9 +1309,16 @@ mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
i++;
} while(i < shape[ULAB_MAX_DIMS - 3]);
#endif
m_del(size_t, shape, ULAB_MAX_DIMS);
m_del(int32_t, strides, ULAB_MAX_DIMS);
m_del(size_t, rshape, ULAB_MAX_DIMS);
m_del(int32_t, rstrides, ULAB_MAX_DIMS);
} else {
mp_raise_TypeError(translate("wrong axis index"));
}
return results;
}

1
python/port/mod/ulab/numpy/numerical.h
View File

@@ -155,6 +155,7 @@
type *_array = (type *)array;\
type tmp;\
uint16_t itmp, c, q = (N), p, r = (N) >> 1;\
assert(N);\
for (;;) {\
if (r > 0) {\
r--;\

98
python/port/mod/ulab/numpy/numpy.c
View File

@@ -8,7 +8,7 @@
*
* Copyright (c) 2020 Jeff Epler for Adafruit Industries
* 2020 Scott Shawcroft for Adafruit Industries
* 2020-2021 Zoltán Vörös
* 2020-2022 Zoltán Vörös
* 2020 Taku Fukada
*/
@@ -17,11 +17,13 @@
#include "py/runtime.h"
#include "numpy.h"
#include "../ulab_create.h"
#include "approx.h"
#include "carray/carray.h"
#include "compare.h"
#include "create.h"
#include "fft/fft.h"
#include "filter.h"
#include "io/io.h"
#include "linalg/linalg.h"
#include "numerical.h"
#include "stats.h"
@@ -125,6 +127,9 @@ static const mp_rom_map_elem_t ulab_numpy_globals_table[] = {
{ MP_ROM_QSTR(MP_QSTR_uint16), MP_ROM_INT(NDARRAY_UINT16) },
{ MP_ROM_QSTR(MP_QSTR_int16), MP_ROM_INT(NDARRAY_INT16) },
{ MP_ROM_QSTR(MP_QSTR_float), MP_ROM_INT(NDARRAY_FLOAT) },
#if ULAB_SUPPORTS_COMPLEX
{ MP_ROM_QSTR(MP_QSTR_complex), MP_ROM_INT(NDARRAY_COMPLEX) },
#endif
// modules of numpy
#if ULAB_NUMPY_HAS_FFT_MODULE
{ MP_ROM_QSTR(MP_QSTR_fft), MP_ROM_PTR(&ulab_fft_module) },
@@ -142,9 +147,15 @@ static const mp_rom_map_elem_t ulab_numpy_globals_table[] = {
#if ULAB_NUMPY_HAS_ARANGE
{ MP_ROM_QSTR(MP_QSTR_arange), (mp_obj_t)&create_arange_obj },
#endif
#if ULAB_NUMPY_HAS_COMPRESS
{ MP_ROM_QSTR(MP_QSTR_compress), (mp_obj_t)&transform_compress_obj },
#endif
#if ULAB_NUMPY_HAS_CONCATENATE
{ MP_ROM_QSTR(MP_QSTR_concatenate), (mp_obj_t)&create_concatenate_obj },
#endif
#if ULAB_NUMPY_HAS_DELETE
{ MP_ROM_QSTR(MP_QSTR_delete), (mp_obj_t)&transform_delete_obj },
#endif
#if ULAB_NUMPY_HAS_DIAG
#if ULAB_MAX_DIMS > 1
{ MP_ROM_QSTR(MP_QSTR_diag), (mp_obj_t)&create_diag_obj },
@@ -224,6 +235,9 @@ static const mp_rom_map_elem_t ulab_numpy_globals_table[] = {
#if ULAB_NUMPY_HAS_ARGSORT
{ MP_OBJ_NEW_QSTR(MP_QSTR_argsort), (mp_obj_t)&numerical_argsort_obj },
#endif
#if ULAB_NUMPY_HAS_ASARRAY
{ MP_OBJ_NEW_QSTR(MP_QSTR_asarray), (mp_obj_t)&create_asarray_obj },
#endif
#if ULAB_NUMPY_HAS_CROSS
{ MP_OBJ_NEW_QSTR(MP_QSTR_cross), (mp_obj_t)&numerical_cross_obj },
#endif
@@ -243,6 +257,12 @@ static const mp_rom_map_elem_t ulab_numpy_globals_table[] = {
#if ULAB_NUMPY_HAS_FLIP
{ MP_OBJ_NEW_QSTR(MP_QSTR_flip), (mp_obj_t)&numerical_flip_obj },
#endif
#if ULAB_NUMPY_HAS_LOAD
{ MP_OBJ_NEW_QSTR(MP_QSTR_load), (mp_obj_t)&io_load_obj },
#endif
#if ULAB_NUMPY_HAS_LOADTXT
{ MP_OBJ_NEW_QSTR(MP_QSTR_loadtxt), (mp_obj_t)&io_loadtxt_obj },
#endif
#if ULAB_NUMPY_HAS_MINMAX
{ MP_OBJ_NEW_QSTR(MP_QSTR_max), (mp_obj_t)&numerical_max_obj },
#endif
@@ -258,6 +278,15 @@ static const mp_rom_map_elem_t ulab_numpy_globals_table[] = {
#if ULAB_NUMPY_HAS_ROLL
{ MP_OBJ_NEW_QSTR(MP_QSTR_roll), (mp_obj_t)&numerical_roll_obj },
#endif
#if ULAB_NUMPY_HAS_SAVE
{ MP_OBJ_NEW_QSTR(MP_QSTR_save), (mp_obj_t)&io_save_obj },
#endif
#if ULAB_NUMPY_HAS_SAVETXT
{ MP_OBJ_NEW_QSTR(MP_QSTR_savetxt), (mp_obj_t)&io_savetxt_obj },
#endif
#if ULAB_NUMPY_HAS_SIZE
{ MP_OBJ_NEW_QSTR(MP_QSTR_size), (mp_obj_t)&transform_size_obj },
#endif
#if ULAB_NUMPY_HAS_SORT
{ MP_OBJ_NEW_QSTR(MP_QSTR_sort), (mp_obj_t)&numerical_sort_obj },
#endif
@@ -276,81 +305,94 @@ static const mp_rom_map_elem_t ulab_numpy_globals_table[] = {
#endif
// functions of the vector sub-module
#if ULAB_NUMPY_HAS_ACOS
{ MP_OBJ_NEW_QSTR(MP_QSTR_acos), (mp_obj_t)&vectorise_acos_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_acos), (mp_obj_t)&vector_acos_obj },
#endif
#if ULAB_NUMPY_HAS_ACOSH
{ MP_OBJ_NEW_QSTR(MP_QSTR_acosh), (mp_obj_t)&vectorise_acosh_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_acosh), (mp_obj_t)&vector_acosh_obj },
#endif
#if ULAB_NUMPY_HAS_ARCTAN2
{ MP_OBJ_NEW_QSTR(MP_QSTR_arctan2), (mp_obj_t)&vectorise_arctan2_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_arctan2), (mp_obj_t)&vector_arctan2_obj },
#endif
#if ULAB_NUMPY_HAS_AROUND
{ MP_OBJ_NEW_QSTR(MP_QSTR_around), (mp_obj_t)&vectorise_around_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_around), (mp_obj_t)&vector_around_obj },
#endif
#if ULAB_NUMPY_HAS_ASIN
{ MP_OBJ_NEW_QSTR(MP_QSTR_asin), (mp_obj_t)&vectorise_asin_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_asin), (mp_obj_t)&vector_asin_obj },
#endif
#if ULAB_NUMPY_HAS_ASINH
{ MP_OBJ_NEW_QSTR(MP_QSTR_asinh), (mp_obj_t)&vectorise_asinh_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_asinh), (mp_obj_t)&vector_asinh_obj },
#endif
#if ULAB_NUMPY_HAS_ATAN
{ MP_OBJ_NEW_QSTR(MP_QSTR_atan), (mp_obj_t)&vectorise_atan_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_atan), (mp_obj_t)&vector_atan_obj },
#endif
#if ULAB_NUMPY_HAS_ATANH
{ MP_OBJ_NEW_QSTR(MP_QSTR_atanh), (mp_obj_t)&vectorise_atanh_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_atanh), (mp_obj_t)&vector_atanh_obj },
#endif
#if ULAB_NUMPY_HAS_CEIL
{ MP_OBJ_NEW_QSTR(MP_QSTR_ceil), (mp_obj_t)&vectorise_ceil_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_ceil), (mp_obj_t)&vector_ceil_obj },
#endif
#if ULAB_NUMPY_HAS_COS
{ MP_OBJ_NEW_QSTR(MP_QSTR_cos), (mp_obj_t)&vectorise_cos_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_cos), (mp_obj_t)&vector_cos_obj },
#endif
#if ULAB_NUMPY_HAS_COSH
{ MP_OBJ_NEW_QSTR(MP_QSTR_cosh), (mp_obj_t)&vectorise_cosh_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_cosh), (mp_obj_t)&vector_cosh_obj },
#endif
#if ULAB_NUMPY_HAS_DEGREES
{ MP_OBJ_NEW_QSTR(MP_QSTR_degrees), (mp_obj_t)&vectorise_degrees_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_degrees), (mp_obj_t)&vector_degrees_obj },
#endif
#if ULAB_NUMPY_HAS_EXP
{ MP_OBJ_NEW_QSTR(MP_QSTR_exp), (mp_obj_t)&vectorise_exp_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_exp), (mp_obj_t)&vector_exp_obj },
#endif
#if ULAB_NUMPY_HAS_EXPM1
{ MP_OBJ_NEW_QSTR(MP_QSTR_expm1), (mp_obj_t)&vectorise_expm1_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_expm1), (mp_obj_t)&vector_expm1_obj },
#endif
#if ULAB_NUMPY_HAS_FLOOR
{ MP_OBJ_NEW_QSTR(MP_QSTR_floor), (mp_obj_t)&vectorise_floor_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_floor), (mp_obj_t)&vector_floor_obj },
#endif
#if ULAB_NUMPY_HAS_LOG
{ MP_OBJ_NEW_QSTR(MP_QSTR_log), (mp_obj_t)&vectorise_log_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_log), (mp_obj_t)&vector_log_obj },
#endif
#if ULAB_NUMPY_HAS_LOG10
{ MP_OBJ_NEW_QSTR(MP_QSTR_log10), (mp_obj_t)&vectorise_log10_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_log10), (mp_obj_t)&vector_log10_obj },
#endif
#if ULAB_NUMPY_HAS_LOG2
{ MP_OBJ_NEW_QSTR(MP_QSTR_log2), (mp_obj_t)&vectorise_log2_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_log2), (mp_obj_t)&vector_log2_obj },
#endif
#if ULAB_NUMPY_HAS_RADIANS
{ MP_OBJ_NEW_QSTR(MP_QSTR_radians), (mp_obj_t)&vectorise_radians_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_radians), (mp_obj_t)&vector_radians_obj },
#endif
#if ULAB_NUMPY_HAS_SIN
{ MP_OBJ_NEW_QSTR(MP_QSTR_sin), (mp_obj_t)&vectorise_sin_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_sin), (mp_obj_t)&vector_sin_obj },
#endif
#if ULAB_NUMPY_HAS_SINH
{ MP_OBJ_NEW_QSTR(MP_QSTR_sinh), (mp_obj_t)&vectorise_sinh_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_sinh), (mp_obj_t)&vector_sinh_obj },
#endif
#if ULAB_NUMPY_HAS_SQRT
{ MP_OBJ_NEW_QSTR(MP_QSTR_sqrt), (mp_obj_t)&vectorise_sqrt_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_sqrt), (mp_obj_t)&vector_sqrt_obj },
#endif
#if ULAB_NUMPY_HAS_TAN
{ MP_OBJ_NEW_QSTR(MP_QSTR_tan), (mp_obj_t)&vectorise_tan_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_tan), (mp_obj_t)&vector_tan_obj },
#endif
#if ULAB_NUMPY_HAS_TANH
{ MP_OBJ_NEW_QSTR(MP_QSTR_tanh), (mp_obj_t)&vectorise_tanh_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_tanh), (mp_obj_t)&vector_tanh_obj },
#endif
#if ULAB_NUMPY_HAS_VECTORIZE
{ MP_OBJ_NEW_QSTR(MP_QSTR_vectorize), (mp_obj_t)&vectorise_vectorize_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_vectorize), (mp_obj_t)&vector_vectorize_obj },
#endif
#if ULAB_SUPPORTS_COMPLEX
#if ULAB_NUMPY_HAS_REAL
{ MP_OBJ_NEW_QSTR(MP_QSTR_real), (mp_obj_t)&carray_real_obj },
#endif
#if ULAB_NUMPY_HAS_IMAG
{ MP_OBJ_NEW_QSTR(MP_QSTR_imag), (mp_obj_t)&carray_imag_obj },
#endif
#if ULAB_NUMPY_HAS_CONJUGATE
{ MP_ROM_QSTR(MP_QSTR_conjugate), (mp_obj_t)&carray_conjugate_obj },
#endif
#if ULAB_NUMPY_HAS_SORT_COMPLEX
{ MP_ROM_QSTR(MP_QSTR_sort_complex), (mp_obj_t)&carray_sort_complex_obj },
#endif
#endif
};
static MP_DEFINE_CONST_DICT(mp_module_ulab_numpy_globals, ulab_numpy_globals_table);

20
python/port/mod/ulab/numpy/poly.c
View File

@@ -19,6 +19,7 @@
#include "../ulab.h"
#include "linalg/linalg_tools.h"
#include "../ulab_tools.h"
#include "carray/carray_tools.h"
#include "poly.h"
#if ULAB_NUMPY_HAS_POLYFIT
@@ -27,6 +28,12 @@ mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {
if(!ndarray_object_is_array_like(args[0])) {
mp_raise_ValueError(translate("input data must be an iterable"));
}
#if ULAB_SUPPORTS_COMPLEX
if(mp_obj_is_type(args[0], &ulab_ndarray_type)) {
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0]);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
}
#endif
size_t lenx = 0, leny = 0;
uint8_t deg = 0;
mp_float_t *x, *XT, *y, *prod;
@@ -142,6 +149,17 @@ mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {
if(!ndarray_object_is_array_like(o_p) || !ndarray_object_is_array_like(o_x)) {
mp_raise_TypeError(translate("inputs are not iterable"));
}
#if ULAB_SUPPORTS_COMPLEX
ndarray_obj_t *input;
if(mp_obj_is_type(o_p, &ulab_ndarray_type)) {
input = MP_OBJ_TO_PTR(o_p);
COMPLEX_DTYPE_NOT_IMPLEMENTED(input->dtype)
}
if(mp_obj_is_type(o_x, &ulab_ndarray_type)) {
input = MP_OBJ_TO_PTR(o_x);
COMPLEX_DTYPE_NOT_IMPLEMENTED(input->dtype)
}
#endif
// p had better be a one-dimensional standard iterable
uint8_t plen = mp_obj_get_int(mp_obj_len_maybe(o_p));
mp_float_t *p = m_new(mp_float_t, plen);
@@ -164,7 +182,7 @@ mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {
mp_float_t (*func)(void *) = ndarray_get_float_function(source->dtype);
// TODO: these loops are really nothing, but the re-implementation of
// TODO: these loops are really nothing, but the re-impplementation of
// ITERATE_VECTOR from vectorise.c. We could pass a function pointer here
#if ULAB_MAX_DIMS > 3
size_t i = 0;

2
python/port/mod/ulab/numpy/stats.c
View File

@@ -21,6 +21,7 @@
#include "../ulab.h"
#include "../ulab_tools.h"
#include "carray/carray_tools.h"
#include "stats.h"
#if ULAB_MAX_DIMS > 1
@@ -36,6 +37,7 @@
static mp_obj_t stats_trace(mp_obj_t oin) {
ndarray_obj_t *ndarray = tools_object_is_square(oin);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
mp_float_t trace = 0.0;
for(size_t i=0; i < ndarray->shape[ULAB_MAX_DIMS - 1]; i++) {
int32_t pos = i * (ndarray->strides[ULAB_MAX_DIMS - 1] + ndarray->strides[ULAB_MAX_DIMS - 2]);

364
python/port/mod/ulab/numpy/transform.c
View File

@@ -9,17 +9,337 @@
*
*/
#include <unistd.h>
#include <math.h>
#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 "carray/carray_tools.h"
#include "numerical.h"
#include "transform.h"
#if ULAB_NUMPY_HAS_COMPRESS
static mp_obj_t transform_compress(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_const_none } },
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_axis, 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 condition = args[0].u_obj;
if(!mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) {
mp_raise_TypeError(translate("wrong input type"));
}
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[1].u_obj);
uint8_t *array = (uint8_t *)ndarray->array;
mp_obj_t axis = args[2].u_obj;
size_t len = MP_OBJ_SMALL_INT_VALUE(mp_obj_len_maybe(condition));
int8_t ax, shift_ax = 0;
if(axis != mp_const_none) {
ax = tools_get_axis(axis, ndarray->ndim);
shift_ax = ULAB_MAX_DIMS - ndarray->ndim + ax;
}
if(((axis == mp_const_none) && (len != ndarray->len)) ||
((axis != mp_const_none) && (len != ndarray->shape[shift_ax]))) {
mp_raise_ValueError(translate("wrong length of condition array"));
}
size_t true_count = 0;
mp_obj_iter_buf_t iter_buf;
mp_obj_t item, iterable = mp_getiter(condition, &iter_buf);
while((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
if(mp_obj_is_true(item)) {
true_count++;
}
}
iterable = mp_getiter(condition, &iter_buf);
ndarray_obj_t *result = NULL;
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
memcpy(shape, ndarray->shape, ULAB_MAX_DIMS * sizeof(size_t));
size_t *rshape = m_new(size_t, ULAB_MAX_DIMS);
memcpy(rshape, ndarray->shape, ULAB_MAX_DIMS * sizeof(size_t));
int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
memcpy(strides, ndarray->strides, ULAB_MAX_DIMS * sizeof(int32_t));
int32_t *rstrides = m_new0(int32_t, ULAB_MAX_DIMS);
if(axis == mp_const_none) {
result = ndarray_new_linear_array(true_count, ndarray->dtype);
rstrides[ULAB_MAX_DIMS - 1] = ndarray->itemsize;
rshape[ULAB_MAX_DIMS - 1] = 0;
} else {
rshape[shift_ax] = true_count;
result = ndarray_new_dense_ndarray(ndarray->ndim, rshape, ndarray->dtype);
SWAP(size_t, shape[shift_ax], shape[ULAB_MAX_DIMS - 1]);
SWAP(size_t, rshape[shift_ax], rshape[ULAB_MAX_DIMS - 1]);
SWAP(int32_t, strides[shift_ax], strides[ULAB_MAX_DIMS - 1]);
memcpy(rstrides, result->strides, ULAB_MAX_DIMS * sizeof(int32_t));
SWAP(int32_t, rstrides[shift_ax], rstrides[ULAB_MAX_DIMS - 1]);
}
uint8_t *rarray = (uint8_t *)result->array;
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
if(axis != mp_const_none) {
iterable = mp_getiter(condition, &iter_buf);
}
do {
item = mp_iternext(iterable);
if(mp_obj_is_true(item)) {
memcpy(rarray, array, ndarray->itemsize);
rarray += rstrides[ULAB_MAX_DIMS - 1];
}
array += strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
array -= strides[ULAB_MAX_DIMS - 1] * shape[ULAB_MAX_DIMS - 1];
array += strides[ULAB_MAX_DIMS - 2];
rarray -= rstrides[ULAB_MAX_DIMS - 1] * rshape[ULAB_MAX_DIMS - 1];
rarray += rstrides[ULAB_MAX_DIMS - 2];
k++;
} while(k < shape[ULAB_MAX_DIMS - 2]);
#endif
#if ULAB_MAX_DIMS > 2
array -= strides[ULAB_MAX_DIMS - 2] * shape[ULAB_MAX_DIMS - 2];
array += strides[ULAB_MAX_DIMS - 3];
rarray -= rstrides[ULAB_MAX_DIMS - 2] * rshape[ULAB_MAX_DIMS - 2];
rarray += rstrides[ULAB_MAX_DIMS - 3];
j++;
} while(j < shape[ULAB_MAX_DIMS - 3]);
#endif
#if ULAB_MAX_DIMS > 3
array -= strides[ULAB_MAX_DIMS - 3] * shape[ULAB_MAX_DIMS - 3];
array += strides[ULAB_MAX_DIMS - 4];
rarray -= rstrides[ULAB_MAX_DIMS - 2] * rshape[ULAB_MAX_DIMS - 2];
rarray += rstrides[ULAB_MAX_DIMS - 3];
i++;
} while(i < shape[ULAB_MAX_DIMS - 4]);
#endif
m_del(size_t, shape, ULAB_MAX_DIMS);
m_del(size_t, rshape, ULAB_MAX_DIMS);
m_del(int32_t, strides, ULAB_MAX_DIMS);
m_del(int32_t, rstrides, ULAB_MAX_DIMS);
return result;
}
MP_DEFINE_CONST_FUN_OBJ_KW(transform_compress_obj, 2, transform_compress);
#endif /* ULAB_NUMPY_HAS_COMPRESS */
#if ULAB_NUMPY_HAS_DELETE
static mp_obj_t transform_delete(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_const_none } },
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = mp_const_none } },
{ MP_QSTR_axis, 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(!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type)) {
mp_raise_TypeError(translate("first argument must be an ndarray"));
}
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
uint8_t *array = (uint8_t *)ndarray->array;
mp_obj_t indices = args[1].u_obj;
mp_obj_t axis = args[2].u_obj;
int8_t shift_ax;
size_t axis_len;
if(axis != mp_const_none) {
int8_t ax = tools_get_axis(axis, ndarray->ndim);
shift_ax = ULAB_MAX_DIMS - ndarray->ndim + ax;
axis_len = ndarray->shape[shift_ax];
} else {
axis_len = ndarray->len;
}
size_t index_len;
if(mp_obj_is_int(indices)) {
index_len = 1;
} else {
if(mp_obj_len_maybe(indices) == MP_OBJ_NULL) {
mp_raise_TypeError(translate("wrong index type"));
}
index_len = MP_OBJ_SMALL_INT_VALUE(mp_obj_len_maybe(indices));
}
if(index_len > axis_len) {
mp_raise_ValueError(translate("wrong length of index array"));
}
size_t *index_array = m_new(size_t, index_len);
if(mp_obj_is_int(indices)) {
ssize_t value = (ssize_t)mp_obj_get_int(indices);
if(value < 0) {
value += axis_len;
}
if((value < 0) || (value > (ssize_t)axis_len)) {
mp_raise_ValueError(translate("index is out of bounds"));
} else {
*index_array++ = (size_t)value;
}
} else {
mp_obj_iter_buf_t iter_buf;
mp_obj_t item, iterable = mp_getiter(indices, &iter_buf);
while((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
ssize_t value = (ssize_t)mp_obj_get_int(item);
if(value < 0) {
value += axis_len;
}
if((value < 0) || (value > (ssize_t)axis_len)) {
mp_raise_ValueError(translate("index is out of bounds"));
} else {
*index_array++ = (size_t)value;
}
}
}
// sort the array, since it is not guaranteed that the input is sorted
HEAPSORT1(size_t, index_array, 1, index_len);
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
memcpy(shape, ndarray->shape, ULAB_MAX_DIMS * sizeof(size_t));
size_t *rshape = m_new(size_t, ULAB_MAX_DIMS);
memcpy(rshape, ndarray->shape, ULAB_MAX_DIMS * sizeof(size_t));
int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
memcpy(strides, ndarray->strides, ULAB_MAX_DIMS * sizeof(int32_t));
int32_t *rstrides = m_new0(int32_t, ULAB_MAX_DIMS);
ndarray_obj_t *result = NULL;
if(axis == mp_const_none) {
result = ndarray_new_linear_array(ndarray->len - index_len, ndarray->dtype);
rstrides[ULAB_MAX_DIMS - 1] = ndarray->itemsize;
memset(rshape, 0, sizeof(size_t) * ULAB_MAX_DIMS);
} else {
rshape[shift_ax] = shape[shift_ax] - index_len;
result = ndarray_new_dense_ndarray(ndarray->ndim, rshape, ndarray->dtype);
SWAP(size_t, shape[shift_ax], shape[ULAB_MAX_DIMS - 1]);
SWAP(size_t, rshape[shift_ax], rshape[ULAB_MAX_DIMS - 1]);
SWAP(int32_t, strides[shift_ax], strides[ULAB_MAX_DIMS - 1]);
memcpy(rstrides, result->strides, ULAB_MAX_DIMS * sizeof(int32_t));
SWAP(int32_t, rstrides[shift_ax], rstrides[ULAB_MAX_DIMS - 1]);
}
uint8_t *rarray = (uint8_t *)result->array;
index_array -= index_len;
size_t count = 0;
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
if(count == *index_array) {
index_array++;
} else {
memcpy(rarray, array, ndarray->itemsize);
rarray += rstrides[ULAB_MAX_DIMS - 1];
}
array += strides[ULAB_MAX_DIMS - 1];
l++;
count++;
} while(l < shape[ULAB_MAX_DIMS - 1]);
if(axis != mp_const_none) {
index_array -= index_len;
count = 0;
}
#if ULAB_MAX_DIMS > 1
array -= strides[ULAB_MAX_DIMS - 1] * shape[ULAB_MAX_DIMS - 1];
array += strides[ULAB_MAX_DIMS - 2];
rarray -= rstrides[ULAB_MAX_DIMS - 1] * rshape[ULAB_MAX_DIMS - 1];
rarray += rstrides[ULAB_MAX_DIMS - 2];
k++;
} while(k < shape[ULAB_MAX_DIMS - 2]);
#endif
#if ULAB_MAX_DIMS > 2
array -= strides[ULAB_MAX_DIMS - 2] * shape[ULAB_MAX_DIMS - 2];
array += strides[ULAB_MAX_DIMS - 3];
rarray -= rstrides[ULAB_MAX_DIMS - 2] * rshape[ULAB_MAX_DIMS - 2];
rarray += rstrides[ULAB_MAX_DIMS - 3];
j++;
} while(j < shape[ULAB_MAX_DIMS - 3]);
#endif
#if ULAB_MAX_DIMS > 3
array -= strides[ULAB_MAX_DIMS - 3] * shape[ULAB_MAX_DIMS - 3];
array += strides[ULAB_MAX_DIMS - 4];
rarray -= rstrides[ULAB_MAX_DIMS - 3] * rshape[ULAB_MAX_DIMS - 3];
rarray += rstrides[ULAB_MAX_DIMS - 4];
i++;
} while(i < shape[ULAB_MAX_DIMS - 4]);
#endif
// TODO: deleting shape generates a seg fault
// m_del(size_t, shape, ULAB_MAX_DIMS);
m_del(size_t, rshape, ULAB_MAX_DIMS);
m_del(int32_t, strides, ULAB_MAX_DIMS);
m_del(int32_t, rstrides, ULAB_MAX_DIMS);
return MP_OBJ_FROM_PTR(result);
}
MP_DEFINE_CONST_FUN_OBJ_KW(transform_delete_obj, 2, transform_delete);
#endif /* ULAB_NUMPY_HAS_DELETE */
#if ULAB_MAX_DIMS > 1
#if ULAB_NUMPY_HAS_DOT
//| def dot(m1: ulab.numpy.ndarray, m2: ulab.numpy.ndarray) -> Union[ulab.numpy.ndarray, _float]:
@@ -39,6 +359,9 @@ mp_obj_t transform_dot(mp_obj_t _m1, mp_obj_t _m2) {
}
ndarray_obj_t *m1 = MP_OBJ_TO_PTR(_m1);
ndarray_obj_t *m2 = MP_OBJ_TO_PTR(_m2);
COMPLEX_DTYPE_NOT_IMPLEMENTED(m1->dtype)
COMPLEX_DTYPE_NOT_IMPLEMENTED(m2->dtype)
uint8_t *array1 = (uint8_t *)m1->array;
uint8_t *array2 = (uint8_t *)m2->array;
@@ -86,5 +409,42 @@ mp_obj_t transform_dot(mp_obj_t _m1, mp_obj_t _m2) {
}
MP_DEFINE_CONST_FUN_OBJ_2(transform_dot_obj, transform_dot);
#endif /* ULAB_NUMPY_HAS_DOT */
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_NUMPY_HAS_SIZE
static mp_obj_t transform_size(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_const_none } },
{ MP_QSTR_axis, 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(ulab_tools_mp_obj_is_scalar(args[0].u_obj)) {
return mp_obj_new_int(1);
}
if(!ndarray_object_is_array_like(args[0].u_obj)) {
mp_raise_TypeError(translate("first argument must be an ndarray"));
}
if(!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type)) {
return mp_obj_len_maybe(args[0].u_obj);
}
// at this point, the args[0] is most certainly an ndarray
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
mp_obj_t axis = args[1].u_obj;
size_t len;
if(axis != mp_const_none) {
int8_t ax = tools_get_axis(axis, ndarray->ndim);
len = ndarray->shape[ULAB_MAX_DIMS - ndarray->ndim + ax];
} else {
len = ndarray->len;
}
return mp_obj_new_int(len);
}
MP_DEFINE_CONST_FUN_OBJ_KW(transform_size_obj, 1, transform_size);
#endif
#endif

4
python/port/mod/ulab/numpy/transform.h
View File

@@ -21,8 +21,10 @@
#include "../ulab.h"
#include "../ulab_tools.h"
#include "transform.h"
MP_DECLARE_CONST_FUN_OBJ_KW(transform_compress_obj);
MP_DECLARE_CONST_FUN_OBJ_KW(transform_delete_obj);
MP_DECLARE_CONST_FUN_OBJ_2(transform_dot_obj);
MP_DECLARE_CONST_FUN_OBJ_KW(transform_size_obj);
#endif

304
python/port/mod/ulab/numpy/vector.c
View File

@@ -22,6 +22,7 @@
#include "../ulab.h"
#include "../ulab_tools.h"
#include "carray/carray_tools.h"
#include "vector.h"
//| """Element-by-element functions
@@ -31,7 +32,7 @@
//| much more efficient than expressing the same operation as a Python loop."""
//|
static mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) {
static mp_obj_t vector_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) {
// Return a single value, if o_in is not iterable
if(mp_obj_is_float(o_in) || mp_obj_is_int(o_in)) {
return mp_obj_new_float(f(mp_obj_get_float(o_in)));
@@ -39,6 +40,7 @@ static mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float
ndarray_obj_t *ndarray = NULL;
if(mp_obj_is_type(o_in, &ulab_ndarray_type)) {
ndarray_obj_t *source = MP_OBJ_TO_PTR(o_in);
COMPLEX_DTYPE_NOT_IMPLEMENTED(source->dtype)
uint8_t *sarray = (uint8_t *)source->array;
ndarray = ndarray_new_dense_ndarray(source->ndim, source->shape, NDARRAY_FLOAT);
mp_float_t *array = (mp_float_t *)ndarray->array;
@@ -99,10 +101,10 @@ static mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float
#endif /* ULAB_VECTORISE_USES_FUN_POINTER */
} else {
ndarray = ndarray_from_mp_obj(o_in, 0);
mp_float_t *array = (mp_float_t *)ndarray->array;
mp_float_t *narray = (mp_float_t *)ndarray->array;
for(size_t i = 0; i < ndarray->len; i++) {
*array = f(*array);
array++;
*narray = f(*narray);
narray++;
}
}
return MP_OBJ_FROM_PTR(ndarray);
@@ -115,7 +117,7 @@ static mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float
//|
MATH_FUN_1(acos, acos);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acos_obj, vectorise_acos);
MP_DEFINE_CONST_FUN_OBJ_1(vector_acos_obj, vector_acos);
#endif
#if ULAB_NUMPY_HAS_ACOSH
@@ -125,7 +127,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acos_obj, vectorise_acos);
//|
MATH_FUN_1(acosh, acosh);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acosh_obj, vectorise_acosh);
MP_DEFINE_CONST_FUN_OBJ_1(vector_acosh_obj, vector_acosh);
#endif
#if ULAB_NUMPY_HAS_ASIN
@@ -135,7 +137,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acosh_obj, vectorise_acosh);
//|
MATH_FUN_1(asin, asin);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_asin_obj, vectorise_asin);
MP_DEFINE_CONST_FUN_OBJ_1(vector_asin_obj, vector_asin);
#endif
#if ULAB_NUMPY_HAS_ASINH
@@ -145,7 +147,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_asin_obj, vectorise_asin);
//|
MATH_FUN_1(asinh, asinh);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_asinh_obj, vectorise_asinh);
MP_DEFINE_CONST_FUN_OBJ_1(vector_asinh_obj, vector_asinh);
#endif
#if ULAB_NUMPY_HAS_AROUND
@@ -155,7 +157,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_asinh_obj, vectorise_asinh);
//| ...
//|
mp_obj_t vectorise_around(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
mp_obj_t vector_around(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_const_none} },
{ MP_QSTR_decimals, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0 } }
@@ -169,6 +171,7 @@ mp_obj_t vectorise_around(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_
int8_t n = args[1].u_int;
mp_float_t mul = MICROPY_FLOAT_C_FUN(pow)(10.0, n);
ndarray_obj_t *source = MP_OBJ_TO_PTR(args[0].u_obj);
COMPLEX_DTYPE_NOT_IMPLEMENTED(source->dtype)
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(source->ndim, source->shape, NDARRAY_FLOAT);
mp_float_t *narray = (mp_float_t *)ndarray->array;
uint8_t *sarray = (uint8_t *)source->array;
@@ -215,7 +218,7 @@ mp_obj_t vectorise_around(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_
return MP_OBJ_FROM_PTR(ndarray);
}
MP_DEFINE_CONST_FUN_OBJ_KW(vectorise_around_obj, 1, vectorise_around);
MP_DEFINE_CONST_FUN_OBJ_KW(vector_around_obj, 1, vector_around);
#endif
#if ULAB_NUMPY_HAS_ATAN
@@ -226,7 +229,7 @@ MP_DEFINE_CONST_FUN_OBJ_KW(vectorise_around_obj, 1, vectorise_around);
//|
MATH_FUN_1(atan, atan);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_atan_obj, vectorise_atan);
MP_DEFINE_CONST_FUN_OBJ_1(vector_atan_obj, vector_atan);
#endif
#if ULAB_NUMPY_HAS_ARCTAN2
@@ -236,9 +239,12 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_atan_obj, vectorise_atan);
//| ...
//|
mp_obj_t vectorise_arctan2(mp_obj_t y, mp_obj_t x) {
mp_obj_t vector_arctan2(mp_obj_t y, mp_obj_t x) {
ndarray_obj_t *ndarray_x = ndarray_from_mp_obj(x, 0);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray_x->dtype)
ndarray_obj_t *ndarray_y = ndarray_from_mp_obj(y, 0);
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray_y->dtype)
uint8_t ndim = 0;
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
@@ -309,7 +315,7 @@ mp_obj_t vectorise_arctan2(mp_obj_t y, mp_obj_t x) {
return MP_OBJ_FROM_PTR(results);
}
MP_DEFINE_CONST_FUN_OBJ_2(vectorise_arctan2_obj, vectorise_arctan2);
MP_DEFINE_CONST_FUN_OBJ_2(vector_arctan2_obj, vector_arctan2);
#endif /* ULAB_VECTORISE_HAS_ARCTAN2 */
#if ULAB_NUMPY_HAS_ATANH
@@ -319,7 +325,7 @@ MP_DEFINE_CONST_FUN_OBJ_2(vectorise_arctan2_obj, vectorise_arctan2);
//|
MATH_FUN_1(atanh, atanh);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_atanh_obj, vectorise_atanh);
MP_DEFINE_CONST_FUN_OBJ_1(vector_atanh_obj, vector_atanh);
#endif
#if ULAB_NUMPY_HAS_CEIL
@@ -329,7 +335,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_atanh_obj, vectorise_atanh);
//|
MATH_FUN_1(ceil, ceil);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_ceil_obj, vectorise_ceil);
MP_DEFINE_CONST_FUN_OBJ_1(vector_ceil_obj, vector_ceil);
#endif
#if ULAB_NUMPY_HAS_COS
@@ -339,7 +345,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_ceil_obj, vectorise_ceil);
//|
MATH_FUN_1(cos, cos);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_cos_obj, vectorise_cos);
MP_DEFINE_CONST_FUN_OBJ_1(vector_cos_obj, vector_cos);
#endif
#if ULAB_NUMPY_HAS_COSH
@@ -349,7 +355,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_cos_obj, vectorise_cos);
//|
MATH_FUN_1(cosh, cosh);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_cosh_obj, vectorise_cosh);
MP_DEFINE_CONST_FUN_OBJ_1(vector_cosh_obj, vector_cosh);
#endif
#if ULAB_NUMPY_HAS_DEGREES
@@ -358,15 +364,15 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_cosh_obj, vectorise_cosh);
//| ...
//|
static mp_float_t vectorise_degrees_(mp_float_t value) {
static mp_float_t vector_degrees_(mp_float_t value) {
return value * MICROPY_FLOAT_CONST(180.0) / MP_PI;
}
static mp_obj_t vectorise_degrees(mp_obj_t x_obj) {
return vectorise_generic_vector(x_obj, vectorise_degrees_);
static mp_obj_t vector_degrees(mp_obj_t x_obj) {
return vector_generic_vector(x_obj, vector_degrees_);
}
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_degrees_obj, vectorise_degrees);
MP_DEFINE_CONST_FUN_OBJ_1(vector_degrees_obj, vector_degrees);
#endif
#if ULAB_SCIPY_SPECIAL_HAS_ERF
@@ -376,7 +382,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_degrees_obj, vectorise_degrees);
//|
MATH_FUN_1(erf, erf);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_erf_obj, vectorise_erf);
MP_DEFINE_CONST_FUN_OBJ_1(vector_erf_obj, vector_erf);
#endif
#if ULAB_SCIPY_SPECIAL_HAS_ERFC
@@ -386,7 +392,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_erf_obj, vectorise_erf);
//|
MATH_FUN_1(erfc, erfc);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_erfc_obj, vectorise_erfc);
MP_DEFINE_CONST_FUN_OBJ_1(vector_erfc_obj, vector_erfc);
#endif
#if ULAB_NUMPY_HAS_EXP
@@ -395,8 +401,69 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_erfc_obj, vectorise_erfc);
//| ...
//|
MATH_FUN_1(exp, exp);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_exp_obj, vectorise_exp);
static mp_obj_t vector_exp(mp_obj_t o_in) {
#if ULAB_SUPPORTS_COMPLEX
if(mp_obj_is_type(o_in, &mp_type_complex)) {
mp_float_t real, imag;
mp_obj_get_complex(o_in, &real, &imag);
mp_float_t exp_real = MICROPY_FLOAT_C_FUN(exp)(real);
return mp_obj_new_complex(exp_real * MICROPY_FLOAT_C_FUN(cos)(imag), exp_real * MICROPY_FLOAT_C_FUN(sin)(imag));
} else if(mp_obj_is_type(o_in, &ulab_ndarray_type)) {
ndarray_obj_t *source = MP_OBJ_TO_PTR(o_in);
if(source->dtype == NDARRAY_COMPLEX) {
uint8_t *sarray = (uint8_t *)source->array;
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(source->ndim, source->shape, NDARRAY_COMPLEX);
mp_float_t *array = (mp_float_t *)ndarray->array;
uint8_t itemsize = sizeof(mp_float_t);
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
mp_float_t real = *(mp_float_t *)sarray;
mp_float_t imag = *(mp_float_t *)(sarray + itemsize);
mp_float_t exp_real = MICROPY_FLOAT_C_FUN(exp)(real);
*array++ = exp_real * MICROPY_FLOAT_C_FUN(cos)(imag);
*array++ = exp_real * MICROPY_FLOAT_C_FUN(sin)(imag);
sarray += source->strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < source->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
sarray -= source->strides[ULAB_MAX_DIMS - 1] * source->shape[ULAB_MAX_DIMS-1];
sarray += source->strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < source->shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_MAX_DIMS > 2
sarray -= source->strides[ULAB_MAX_DIMS - 2] * source->shape[ULAB_MAX_DIMS-2];
sarray += source->strides[ULAB_MAX_DIMS - 3];
j++;
} while(j < source->shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS > 2 */
#if ULAB_MAX_DIMS > 3
sarray -= source->strides[ULAB_MAX_DIMS - 3] * source->shape[ULAB_MAX_DIMS-3];
sarray += source->strides[ULAB_MAX_DIMS - 4];
i++;
} while(i < source->shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS > 3 */
return MP_OBJ_FROM_PTR(ndarray);
}
}
#endif
return vector_generic_vector(o_in, MICROPY_FLOAT_C_FUN(exp));
}
MP_DEFINE_CONST_FUN_OBJ_1(vector_exp_obj, vector_exp);
#endif
#if ULAB_NUMPY_HAS_EXPM1
@@ -406,7 +473,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_exp_obj, vectorise_exp);
//|
MATH_FUN_1(expm1, expm1);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_expm1_obj, vectorise_expm1);
MP_DEFINE_CONST_FUN_OBJ_1(vector_expm1_obj, vector_expm1);
#endif
#if ULAB_NUMPY_HAS_FLOOR
@@ -416,7 +483,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_expm1_obj, vectorise_expm1);
//|
MATH_FUN_1(floor, floor);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_floor_obj, vectorise_floor);
MP_DEFINE_CONST_FUN_OBJ_1(vector_floor_obj, vector_floor);
#endif
#if ULAB_SCIPY_SPECIAL_HAS_GAMMA
@@ -426,7 +493,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_floor_obj, vectorise_floor);
//|
MATH_FUN_1(gamma, tgamma);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_gamma_obj, vectorise_gamma);
MP_DEFINE_CONST_FUN_OBJ_1(vector_gamma_obj, vector_gamma);
#endif
#if ULAB_SCIPY_SPECIAL_HAS_GAMMALN
@@ -436,7 +503,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_gamma_obj, vectorise_gamma);
//|
MATH_FUN_1(lgamma, lgamma);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_lgamma_obj, vectorise_lgamma);
MP_DEFINE_CONST_FUN_OBJ_1(vector_lgamma_obj, vector_lgamma);
#endif
#if ULAB_NUMPY_HAS_LOG
@@ -446,7 +513,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_lgamma_obj, vectorise_lgamma);
//|
MATH_FUN_1(log, log);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log_obj, vectorise_log);
MP_DEFINE_CONST_FUN_OBJ_1(vector_log_obj, vector_log);
#endif
#if ULAB_NUMPY_HAS_LOG10
@@ -456,7 +523,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log_obj, vectorise_log);
//|
MATH_FUN_1(log10, log10);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log10_obj, vectorise_log10);
MP_DEFINE_CONST_FUN_OBJ_1(vector_log10_obj, vector_log10);
#endif
#if ULAB_NUMPY_HAS_LOG2
@@ -466,7 +533,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log10_obj, vectorise_log10);
//|
MATH_FUN_1(log2, log2);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log2_obj, vectorise_log2);
MP_DEFINE_CONST_FUN_OBJ_1(vector_log2_obj, vector_log2);
#endif
#if ULAB_NUMPY_HAS_RADIANS
@@ -475,15 +542,15 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log2_obj, vectorise_log2);
//| ...
//|
static mp_float_t vectorise_radians_(mp_float_t value) {
static mp_float_t vector_radians_(mp_float_t value) {
return value * MP_PI / MICROPY_FLOAT_CONST(180.0);
}
static mp_obj_t vectorise_radians(mp_obj_t x_obj) {
return vectorise_generic_vector(x_obj, vectorise_radians_);
static mp_obj_t vector_radians(mp_obj_t x_obj) {
return vector_generic_vector(x_obj, vector_radians_);
}
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_radians_obj, vectorise_radians);
MP_DEFINE_CONST_FUN_OBJ_1(vector_radians_obj, vector_radians);
#endif
#if ULAB_NUMPY_HAS_SIN
@@ -493,7 +560,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_radians_obj, vectorise_radians);
//|
MATH_FUN_1(sin, sin);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sin_obj, vectorise_sin);
MP_DEFINE_CONST_FUN_OBJ_1(vector_sin_obj, vector_sin);
#endif
#if ULAB_NUMPY_HAS_SINH
@@ -503,18 +570,158 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sin_obj, vectorise_sin);
//|
MATH_FUN_1(sinh, sinh);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sinh_obj, vectorise_sinh);
MP_DEFINE_CONST_FUN_OBJ_1(vector_sinh_obj, vector_sinh);
#endif
#if ULAB_NUMPY_HAS_SQRT
//| def sqrt(a: _ArrayLike) -> ulab.numpy.ndarray:
//| """Computes the square root"""
//| ...
//|
#if ULAB_SUPPORTS_COMPLEX
mp_obj_t vector_sqrt(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_const_none } },
{ MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_INT(NDARRAY_FLOAT) } },
};
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 o_in = args[0].u_obj;
uint8_t dtype = mp_obj_get_int(args[1].u_obj);
if((dtype != NDARRAY_FLOAT) && (dtype != NDARRAY_COMPLEX)) {
mp_raise_TypeError(translate("dtype must be float, or complex"));
}
if(mp_obj_is_type(o_in, &mp_type_complex)) {
mp_float_t real, imag;
mp_obj_get_complex(o_in, &real, &imag);
mp_float_t sqrt_abs = MICROPY_FLOAT_C_FUN(sqrt)(real * real + imag * imag);
sqrt_abs = MICROPY_FLOAT_C_FUN(sqrt)(sqrt_abs);
mp_float_t theta = MICROPY_FLOAT_CONST(0.5) * MICROPY_FLOAT_C_FUN(atan2)(imag, real);
return mp_obj_new_complex(sqrt_abs * MICROPY_FLOAT_C_FUN(cos)(theta), sqrt_abs * MICROPY_FLOAT_C_FUN(sin)(theta));
} else if(mp_obj_is_type(o_in, &ulab_ndarray_type)) {
ndarray_obj_t *source = MP_OBJ_TO_PTR(o_in);
if((source->dtype == NDARRAY_COMPLEX) && (dtype == NDARRAY_FLOAT)) {
mp_raise_TypeError(translate("can't convert complex to float"));
}
if(dtype == NDARRAY_COMPLEX) {
if(source->dtype == NDARRAY_COMPLEX) {
uint8_t *sarray = (uint8_t *)source->array;
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(source->ndim, source->shape, NDARRAY_COMPLEX);
mp_float_t *array = (mp_float_t *)ndarray->array;
uint8_t itemsize = sizeof(mp_float_t);
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
mp_float_t real = *(mp_float_t *)sarray;
mp_float_t imag = *(mp_float_t *)(sarray + itemsize);
mp_float_t sqrt_abs = MICROPY_FLOAT_C_FUN(sqrt)(real * real + imag * imag);
sqrt_abs = MICROPY_FLOAT_C_FUN(sqrt)(sqrt_abs);
mp_float_t theta = MICROPY_FLOAT_CONST(0.5) * MICROPY_FLOAT_C_FUN(atan2)(imag, real);
*array++ = sqrt_abs * MICROPY_FLOAT_C_FUN(cos)(theta);
*array++ = sqrt_abs * MICROPY_FLOAT_C_FUN(sin)(theta);
sarray += source->strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < source->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
sarray -= source->strides[ULAB_MAX_DIMS - 1] * source->shape[ULAB_MAX_DIMS-1];
sarray += source->strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < source->shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_MAX_DIMS > 2
sarray -= source->strides[ULAB_MAX_DIMS - 2] * source->shape[ULAB_MAX_DIMS-2];
sarray += source->strides[ULAB_MAX_DIMS - 3];
j++;
} while(j < source->shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS > 2 */
#if ULAB_MAX_DIMS > 3
sarray -= source->strides[ULAB_MAX_DIMS - 3] * source->shape[ULAB_MAX_DIMS-3];
sarray += source->strides[ULAB_MAX_DIMS - 4];
i++;
} while(i < source->shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS > 3 */
return MP_OBJ_FROM_PTR(ndarray);
} else if(source->dtype == NDARRAY_FLOAT) {
uint8_t *sarray = (uint8_t *)source->array;
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(source->ndim, source->shape, NDARRAY_COMPLEX);
mp_float_t *array = (mp_float_t *)ndarray->array;
#if ULAB_MAX_DIMS > 3
size_t i = 0;
do {
#endif
#if ULAB_MAX_DIMS > 2
size_t j = 0;
do {
#endif
#if ULAB_MAX_DIMS > 1
size_t k = 0;
do {
#endif
size_t l = 0;
do {
mp_float_t value = *(mp_float_t *)sarray;
if(value >= MICROPY_FLOAT_CONST(0.0)) {
*array++ = MICROPY_FLOAT_C_FUN(sqrt)(value);
array++;
} else {
array++;
*array++ = MICROPY_FLOAT_C_FUN(sqrt)(-value);
}
sarray += source->strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < source->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
sarray -= source->strides[ULAB_MAX_DIMS - 1] * source->shape[ULAB_MAX_DIMS-1];
sarray += source->strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < source->shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS > 1 */
#if ULAB_MAX_DIMS > 2
sarray -= source->strides[ULAB_MAX_DIMS - 2] * source->shape[ULAB_MAX_DIMS-2];
sarray += source->strides[ULAB_MAX_DIMS - 3];
j++;
} while(j < source->shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS > 2 */
#if ULAB_MAX_DIMS > 3
sarray -= source->strides[ULAB_MAX_DIMS - 3] * source->shape[ULAB_MAX_DIMS-3];
sarray += source->strides[ULAB_MAX_DIMS - 4];
i++;
} while(i < source->shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS > 3 */
return MP_OBJ_FROM_PTR(ndarray);
} else {
mp_raise_TypeError(translate("input dtype must be float or complex"));
}
}
}
return vector_generic_vector(o_in, MICROPY_FLOAT_C_FUN(sqrt));
}
MP_DEFINE_CONST_FUN_OBJ_KW(vector_sqrt_obj, 1, vector_sqrt);
#else
MATH_FUN_1(sqrt, sqrt);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sqrt_obj, vectorise_sqrt);
#endif
MP_DEFINE_CONST_FUN_OBJ_1(vector_sqrt_obj, vector_sqrt);
#endif /* ULAB_SUPPORTS_COMPLEX */
#endif /* ULAB_NUMPY_HAS_SQRT */
#if ULAB_NUMPY_HAS_TAN
//| def tan(a: _ArrayLike) -> ulab.numpy.ndarray:
@@ -523,7 +730,7 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sqrt_obj, vectorise_sqrt);
//|
MATH_FUN_1(tan, tan);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tan_obj, vectorise_tan);
MP_DEFINE_CONST_FUN_OBJ_1(vector_tan_obj, vector_tan);
#endif
#if ULAB_NUMPY_HAS_TANH
@@ -532,11 +739,11 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tan_obj, vectorise_tan);
//| ...
MATH_FUN_1(tanh, tanh);
MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tanh_obj, vectorise_tanh);
MP_DEFINE_CONST_FUN_OBJ_1(vector_tanh_obj, vector_tanh);
#endif
#if ULAB_NUMPY_HAS_VECTORIZE
static mp_obj_t vectorise_vectorized_function_call(mp_obj_t self_in, size_t n_args, size_t n_kw, const mp_obj_t *args) {
static mp_obj_t vector_vectorized_function_call(mp_obj_t self_in, size_t n_args, size_t n_kw, const mp_obj_t *args) {
(void) n_args;
(void) n_kw;
vectorized_function_obj_t *self = MP_OBJ_TO_PTR(self_in);
@@ -544,6 +751,7 @@ static mp_obj_t vectorise_vectorized_function_call(mp_obj_t self_in, size_t n_ar
mp_obj_t fvalue;
if(mp_obj_is_type(args[0], &ulab_ndarray_type)) {
ndarray_obj_t *source = MP_OBJ_TO_PTR(args[0]);
COMPLEX_DTYPE_NOT_IMPLEMENTED(source->dtype)
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(source->ndim, source->shape, self->otypes);
for(size_t i=0; i < source->len; i++) {
avalue[0] = mp_binary_get_val_array(source->dtype, source->array, i);
@@ -575,12 +783,12 @@ static mp_obj_t vectorise_vectorized_function_call(mp_obj_t self_in, size_t n_ar
return mp_const_none;
}
const mp_obj_type_t vectorise_function_type = {
const mp_obj_type_t vector_function_type = {
{ &mp_type_type },
.flags = MP_TYPE_FLAG_EXTENDED,
.name = MP_QSTR_,
MP_TYPE_EXTENDED_FIELDS(
.call = vectorise_vectorized_function_call,
.call = vector_vectorized_function_call,
)
};
@@ -598,7 +806,7 @@ const mp_obj_type_t vectorise_function_type = {
//| ...
//|
static mp_obj_t vectorise_vectorize(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
static mp_obj_t vector_vectorize(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_const_none} },
{ MP_QSTR_otypes, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} }
@@ -625,12 +833,12 @@ static mp_obj_t vectorise_vectorize(size_t n_args, const mp_obj_t *pos_args, mp_
mp_raise_ValueError(translate("wrong output type"));
}
vectorized_function_obj_t *function = m_new_obj(vectorized_function_obj_t);
function->base.type = &vectorise_function_type;
function->base.type = &vector_function_type;
function->otypes = otypes;
function->fun = args[0].u_obj;
function->type = type;
return MP_OBJ_FROM_PTR(function);
}
MP_DEFINE_CONST_FUN_OBJ_KW(vectorise_vectorize_obj, 1, vectorise_vectorize);
MP_DEFINE_CONST_FUN_OBJ_KW(vector_vectorize_obj, 1, vector_vectorize);
#endif

71
python/port/mod/ulab/numpy/vector.h
View File

@@ -15,35 +15,39 @@
#include "../ulab.h"
#include "../ndarray.h"
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_acos_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_acosh_obj);
MP_DECLARE_CONST_FUN_OBJ_2(vectorise_arctan2_obj);
MP_DECLARE_CONST_FUN_OBJ_KW(vectorise_around_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_asin_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_asinh_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_atan_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_atanh_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_ceil_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_cos_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_cosh_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_degrees_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_erf_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_erfc_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_exp_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_expm1_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_floor_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_gamma_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_lgamma_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_log_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_log10_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_log2_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_radians_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_sin_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_sinh_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_sqrt_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_tan_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vectorise_tanh_obj);
MP_DECLARE_CONST_FUN_OBJ_KW(vectorise_vectorize_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_acos_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_acosh_obj);
MP_DECLARE_CONST_FUN_OBJ_2(vector_arctan2_obj);
MP_DECLARE_CONST_FUN_OBJ_KW(vector_around_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_asin_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_asinh_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_atan_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_atanh_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_ceil_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_cos_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_cosh_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_degrees_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_erf_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_erfc_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_exp_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_expm1_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_floor_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_gamma_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_lgamma_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_log_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_log10_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_log2_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_radians_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_sin_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_sinh_obj);
#if ULAB_SUPPORTS_COMPLEX
MP_DECLARE_CONST_FUN_OBJ_KW(vector_sqrt_obj);
#else
MP_DECLARE_CONST_FUN_OBJ_1(vector_sqrt_obj);
#endif
MP_DECLARE_CONST_FUN_OBJ_1(vector_tan_obj);
MP_DECLARE_CONST_FUN_OBJ_1(vector_tanh_obj);
MP_DECLARE_CONST_FUN_OBJ_KW(vector_vectorize_obj);
typedef struct _vectorized_function_obj_t {
mp_obj_base_t base;
@@ -53,12 +57,13 @@ typedef struct _vectorized_function_obj_t {
} vectorized_function_obj_t;
#if ULAB_HAS_FUNCTION_ITERATOR
#define ITERATE_VECTOR(type, array, source, sarray)\
#define ITERATE_VECTOR(type, array, source, sarray, shift)\
({\
size_t *scoords = ndarray_new_coords((source)->ndim);\
for(size_t i=0; i < (source)->len/(source)->shape[ULAB_MAX_DIMS -1]; i++) {\
for(size_t l=0; l < (source)->shape[ULAB_MAX_DIMS - 1]; l++) {\
*(array)++ = f(*((type *)(sarray)));\
*(array) = f(*((type *)(sarray)));\
(array) += (shift);\
(sarray) += (source)->strides[ULAB_MAX_DIMS - 1];\
}\
ndarray_rewind_array((source)->ndim, sarray, (source)->shape, (source)->strides, scoords);\
@@ -149,8 +154,8 @@ typedef struct _vectorized_function_obj_t {
#endif /* ULAB_HAS_FUNCTION_ITERATOR */
#define MATH_FUN_1(py_name, c_name) \
static mp_obj_t vectorise_ ## py_name(mp_obj_t x_obj) { \
return vectorise_generic_vector(x_obj, MICROPY_FLOAT_C_FUN(c_name)); \
static mp_obj_t vector_ ## py_name(mp_obj_t x_obj) { \
return vector_generic_vector(x_obj, MICROPY_FLOAT_C_FUN(c_name)); \
}
#endif /* _VECTOR_ */

4
python/port/mod/ulab/scipy/optimize/optimize.c
View File

@@ -121,7 +121,7 @@ MP_DEFINE_CONST_FUN_OBJ_KW(optimize_bisect_obj, 3, optimize_bisect);
//| 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 required."""
//| steps are requried."""
//| ...
//|
@@ -344,7 +344,7 @@ MP_DEFINE_CONST_FUN_OBJ_KW(optimize_curve_fit_obj, 2, optimize_curve_fit);
//|
//| 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 required."""
//| ``maxiter`` steps are requried."""
//| ...
//|

2
python/port/mod/ulab/scipy/scipy.c
View File

@@ -48,4 +48,4 @@ mp_obj_module_t ulab_scipy_module = {
.base = { &mp_type_module },
.globals = (mp_obj_dict_t*)&mp_module_ulab_scipy_globals,
};
#endif
#endif /* ULAB_HAS_SCIPY */

40
python/port/mod/ulab/scipy/signal/signal.c
View File

@@ -18,32 +18,9 @@
#include "../../ulab.h"
#include "../../ndarray.h"
#include "../../numpy/fft/fft_tools.h"
#include "../../numpy/carray/carray_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(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);
}
}
MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(signal_spectrogram_obj, 1, 2, signal_spectrogram);
#endif /* ULAB_SCIPY_SIGNAL_HAS_SPECTROGRAM */
#if ULAB_SCIPY_SIGNAL_HAS_SOSFILT
#if ULAB_SCIPY_SIGNAL_HAS_SOSFILT & ULAB_MAX_DIMS > 1
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;
@@ -68,6 +45,12 @@ mp_obj_t signal_sosfilt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
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;
@@ -102,7 +85,7 @@ mp_obj_t signal_sosfilt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
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)) {
if((zi->shape[ULAB_MAX_DIMS - 2] != 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) {
@@ -139,10 +122,7 @@ MP_DEFINE_CONST_FUN_OBJ_KW(signal_sosfilt_obj, 2, signal_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
#if ULAB_SCIPY_SIGNAL_HAS_SOSFILT & ULAB_MAX_DIMS > 1
{ MP_OBJ_NEW_QSTR(MP_QSTR_sosfilt), (mp_obj_t)&signal_sosfilt_obj },
#endif
};

1
python/port/mod/ulab/scipy/signal/signal.h
View File

@@ -18,7 +18,6 @@
extern 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_ */

8
python/port/mod/ulab/scipy/special/special.c
View File

@@ -21,16 +21,16 @@
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)&vectorise_erf_obj },
{ 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)&vectorise_erfc_obj },
{ 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)&vectorise_gamma_obj },
{ 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)&vectorise_lgamma_obj },
{ MP_OBJ_NEW_QSTR(MP_QSTR_gammaln), (mp_obj_t)&vector_lgamma_obj },
#endif
};

24
python/port/mod/ulab/ulab.c
View File

@@ -20,9 +20,9 @@
#include "py/objarray.h"
#include "ulab.h"
#include "ulab_create.h"
#include "ndarray.h"
#include "ndarray_properties.h"
#include "numpy/create.h"
#include "numpy/ndarray/ndarray_iter.h"
#include "numpy/numpy.h"
@@ -33,13 +33,21 @@
#include "user/user.h"
#include "utils/utils.h"
#define ULAB_VERSION 3.3.8
#define ULAB_VERSION 5.0.7
#define xstr(s) str(s)
#define str(s) #s
#if ULAB_SUPPORTS_COMPLEX
#define ULAB_VERSION_STRING xstr(ULAB_VERSION) xstr(-) xstr(ULAB_MAX_DIMS) xstr(D-c)
#else
#define ULAB_VERSION_STRING xstr(ULAB_VERSION) xstr(-) xstr(ULAB_MAX_DIMS) xstr(D)
#endif
STATIC MP_DEFINE_STR_OBJ(ulab_version_obj, ULAB_VERSION_STRING);
#ifdef ULAB_HASH
STATIC MP_DEFINE_STR_OBJ(ulab_sha_obj, xstr(ULAB_HASH));
#endif
STATIC const mp_rom_map_elem_t ulab_ndarray_locals_dict_table[] = {
#if ULAB_MAX_DIMS > 1
@@ -62,6 +70,9 @@ STATIC const mp_rom_map_elem_t ulab_ndarray_locals_dict_table[] = {
#if NDARRAY_HAS_TOBYTES
{ MP_ROM_QSTR(MP_QSTR_tobytes), MP_ROM_PTR(&ndarray_tobytes_obj) },
#endif
#if NDARRAY_HAS_TOLIST
{ MP_ROM_QSTR(MP_QSTR_tolist), MP_ROM_PTR(&ndarray_tolist_obj) },
#endif
#if NDARRAY_HAS_SORT
{ MP_ROM_QSTR(MP_QSTR_sort), MP_ROM_PTR(&numerical_sort_inplace_obj) },
#endif
@@ -141,6 +152,9 @@ const mp_obj_type_t ndarray_flatiter_type = {
STATIC const mp_map_elem_t ulab_globals_table[] = {
{ MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_ulab) },
{ MP_ROM_QSTR(MP_QSTR___version__), MP_ROM_PTR(&ulab_version_obj) },
#ifdef ULAB_HASH
{ MP_ROM_QSTR(MP_QSTR___sha__), MP_ROM_PTR(&ulab_sha_obj) },
#endif
#if ULAB_HAS_DTYPE_OBJECT
{ MP_OBJ_NEW_QSTR(MP_QSTR_dtype), (mp_obj_t)&ulab_dtype_type },
#else
@@ -174,4 +188,10 @@ const mp_obj_module_t ulab_user_cmodule = {
.globals = (mp_obj_dict_t*)&mp_module_ulab_globals,
};
// Use old three-argument MP_REGISTER_MODULE for
// MicroPython <= v1.18.0: (1 << 16) | (18 << 8) | 0
#if MICROPY_VERSION <= 70144
MP_REGISTER_MODULE(MP_QSTR_ulab, ulab_user_cmodule, MODULE_ULAB_ENABLED);
#else
MP_REGISTER_MODULE(MP_QSTR_ulab, ulab_user_cmodule);
#endif

93
python/port/mod/ulab/ulab.h
View File

@@ -6,7 +6,7 @@
*
* The MIT License (MIT)
*
* Copyright (c) 2019-2021 Zoltán Vörös
* Copyright (c) 2019-2022 Zoltán Vörös
*/
#ifndef __ULAB__
@@ -18,9 +18,9 @@
//
// - how many dimensions ulab can handle
// - which functions are included in the compiled firmware
// - whether the python syntax is numpy-like, or modular
// - whether arrays can be sliced and iterated over
// - which binary/unary operators are supported
// - whether ulab can deal with complex numbers
//
// A considerable amount of flash space can be saved by removing (setting
// the corresponding constants to 0) the unnecessary functions and features.
@@ -31,6 +31,10 @@
#include ULAB_CONFIG_FILE
#endif
// Adds support for complex ndarrays
#ifndef ULAB_SUPPORTS_COMPLEX
#define ULAB_SUPPORTS_COMPLEX (0)
#endif
// Determines, whether scipy is defined in ulab. The sub-modules and functions
// of scipy have to be defined separately
@@ -228,6 +232,10 @@
#define NDARRAY_HAS_TOBYTES (1)
#endif
#ifndef NDARRAY_HAS_TOLIST
#define NDARRAY_HAS_TOLIST (1)
#endif
#ifndef NDARRAY_HAS_TRANSPOSE
#define NDARRAY_HAS_TRANSPOSE (1)
#endif
@@ -385,6 +393,15 @@
#define ULAB_NUMPY_HAS_FFT_MODULE (1)
#endif
// By setting this constant to 1, the FFT routine will behave in a
// numpy-compatible way, i.e., it will output a complex array
// This setting has no effect, if ULAB_SUPPORTS_COMPLEX is 0
// Note that in this case, the input also must be numpythonic,
// i.e., the real an imaginary parts cannot be passed as two arguments
#ifndef ULAB_FFT_IS_NUMPY_COMPATIBLE
#define ULAB_FFT_IS_NUMPY_COMPATIBLE (0)
#endif
#ifndef ULAB_FFT_HAS_FFT
#define ULAB_FFT_HAS_FFT (1)
#endif
@@ -409,6 +426,14 @@
#define ULAB_NUMPY_HAS_ARGSORT (1)
#endif
#ifndef ULAB_NUMPY_HAS_ASARRAY
#define ULAB_NUMPY_HAS_ASARRAY (0)
#endif
#ifndef ULAB_NUMPY_HAS_COMPRESS
#define ULAB_NUMPY_HAS_COMPRESS (1)
#endif
#ifndef ULAB_NUMPY_HAS_CONVOLVE
#define ULAB_NUMPY_HAS_CONVOLVE (1)
#endif
@@ -417,6 +442,10 @@
#define ULAB_NUMPY_HAS_CROSS (1)
#endif
#ifndef ULAB_NUMPY_HAS_DELETE
#define ULAB_NUMPY_HAS_DELETE (1)
#endif
#ifndef ULAB_NUMPY_HAS_DIFF
#define ULAB_NUMPY_HAS_DIFF (1)
#endif
@@ -433,6 +462,14 @@
#define ULAB_NUMPY_HAS_INTERP (1)
#endif
#ifndef ULAB_NUMPY_HAS_LOAD
#define ULAB_NUMPY_HAS_LOAD (0)
#endif
#ifndef ULAB_NUMPY_HAS_LOADTXT
#define ULAB_NUMPY_HAS_LOADTXT (0)
#endif
#ifndef ULAB_NUMPY_HAS_MEAN
#define ULAB_NUMPY_HAS_MEAN (1)
#endif
@@ -457,6 +494,18 @@
#define ULAB_NUMPY_HAS_ROLL (1)
#endif
#ifndef ULAB_NUMPY_HAS_SAVE
#define ULAB_NUMPY_HAS_SAVE (0)
#endif
#ifndef ULAB_NUMPY_HAS_SAVETXT
#define ULAB_NUMPY_HAS_SAVETXT (0)
#endif
#ifndef ULAB_NUMPY_HAS_SIZE
#define ULAB_NUMPY_HAS_SIZE (1)
#endif
#ifndef ULAB_NUMPY_HAS_SORT
#define ULAB_NUMPY_HAS_SORT (1)
#endif
@@ -579,6 +628,25 @@
#define ULAB_NUMPY_HAS_VECTORIZE (1)
#endif
// Complex functions. The implementations are compiled into
// the firmware, only if ULAB_SUPPORTS_COMPLEX is set to 1
#ifndef ULAB_NUMPY_HAS_CONJUGATE
#define ULAB_NUMPY_HAS_CONJUGATE (1)
#endif
#ifndef ULAB_NUMPY_HAS_IMAG
#define ULAB_NUMPY_HAS_IMAG (1)
#endif
#ifndef ULAB_NUMPY_HAS_REAL
#define ULAB_NUMPY_HAS_REAL (1)
#endif
#ifndef ULAB_NUMPY_HAS_SORT_COMPLEX
#define ULAB_NUMPY_HAS_SORT_COMPLEX (0)
#endif
// scipy modules
#ifndef ULAB_SCIPY_HAS_LINALG_MODULE
#define ULAB_SCIPY_HAS_LINALG_MODULE (1)
#endif
@@ -595,10 +663,6 @@
#define ULAB_SCIPY_HAS_SIGNAL_MODULE (1)
#endif
#ifndef ULAB_SCIPY_SIGNAL_HAS_SPECTROGRAM
#define ULAB_SCIPY_SIGNAL_HAS_SPECTROGRAM (1)
#endif
#ifndef ULAB_SCIPY_SIGNAL_HAS_SOSFILT
#define ULAB_SCIPY_SIGNAL_HAS_SOSFILT (1)
#endif
@@ -643,12 +707,7 @@
#define ULAB_SCIPY_SPECIAL_HAS_GAMMALN (1)
#endif
// user-defined module; source of the module and
// its sub-modules should be placed in code/user/
#ifndef ULAB_HAS_USER_MODULE
#define ULAB_HAS_USER_MODULE (0)
#endif
// functions of the utils module
#ifndef ULAB_HAS_UTILS_MODULE
#define ULAB_HAS_UTILS_MODULE (1)
#endif
@@ -669,4 +728,14 @@
#define ULAB_UTILS_HAS_FROM_UINT32_BUFFER (1)
#endif
#ifndef ULAB_UTILS_HAS_SPECTROGRAM
#define ULAB_UTILS_HAS_SPECTROGRAM (1)
#endif
// user-defined module; source of the module and
// its sub-modules should be placed in code/user/
#ifndef ULAB_HAS_USER_MODULE
#define ULAB_HAS_USER_MODULE (0)
#endif
#endif

45
python/port/mod/ulab/ulab_tools.c
View File

@@ -5,7 +5,7 @@
*
* The MIT License (MIT)
*
* Copyright (c) 2020-2021 Zoltán Vörös
* Copyright (c) 2020-2022 Zoltán Vörös
*/
@@ -216,6 +216,14 @@ shape_strides tools_reduce_axes(ndarray_obj_t *ndarray, mp_obj_t axis) {
return _shape_strides;
}
int8_t tools_get_axis(mp_obj_t axis, uint8_t ndim) {
int8_t ax = mp_obj_get_int(axis);
if(ax < 0) ax += ndim;
if((ax < 0) || (ax > ndim - 1)) {
mp_raise_ValueError(translate("axis is out of bounds"));
}
return ax;
}
#if ULAB_MAX_DIMS > 1
ndarray_obj_t *tools_object_is_square(mp_obj_t obj) {
@@ -231,3 +239,38 @@ ndarray_obj_t *tools_object_is_square(mp_obj_t obj) {
return ndarray;
}
#endif
uint8_t ulab_binary_get_size(uint8_t dtype) {
#if ULAB_SUPPORTS_COMPLEX
if(dtype == NDARRAY_COMPLEX) {
return 2 * (uint8_t)sizeof(mp_float_t);
}
#endif
return dtype == NDARRAY_BOOL ? 1 : mp_binary_get_size('@', dtype, NULL);
}
#if ULAB_SUPPORTS_COMPLEX
void ulab_rescale_float_strides(int32_t *strides) {
// re-scale the strides, so that we can work with floats, when iterating
uint8_t sz = sizeof(mp_float_t);
for(uint8_t i = 0; i < ULAB_MAX_DIMS; i++) {
strides[i] /= sz;
}
}
#endif
bool ulab_tools_mp_obj_is_scalar(mp_obj_t obj) {
#if ULAB_SUPPORTS_COMPLEX
if(mp_obj_is_int(obj) || mp_obj_is_float(obj) || mp_obj_is_type(obj, &mp_type_complex)) {
return true;
} else {
return false;
}
#else
if(mp_obj_is_int(obj) || mp_obj_is_float(obj)) {
return true;
} else {
return false;
}
#endif
}

11
python/port/mod/ulab/ulab_tools.h
View File

@@ -5,7 +5,7 @@
*
* The MIT License (MIT)
*
* Copyright (c) 2020-2021 Zoltán Vörös
* Copyright (c) 2020-2022 Zoltán Vörös
*/
#ifndef _TOOLS_
@@ -33,5 +33,14 @@ uint8_t ndarray_upcast_dtype(uint8_t , uint8_t );
void *ndarray_set_float_function(uint8_t );
shape_strides tools_reduce_axes(ndarray_obj_t *, mp_obj_t );
int8_t tools_get_axis(mp_obj_t , uint8_t );
ndarray_obj_t *tools_object_is_square(mp_obj_t );
uint8_t ulab_binary_get_size(uint8_t );
#if ULAB_SUPPORTS_COMPLEX
void ulab_rescale_float_strides(int32_t *);
#endif
bool ulab_tools_mp_obj_is_scalar(mp_obj_t );
#endif

2
python/port/mod/ulab/user/user.c
View File

@@ -74,7 +74,7 @@ static mp_obj_t user_square(mp_obj_t arg) {
*rarray++ = (*array) * (*array);
}
}
// at the end, return a micropython object
// at the end, return a micrppython object
return MP_OBJ_FROM_PTR(results);
}

40
python/port/mod/ulab/utils/utils.c
View File

@@ -16,6 +16,8 @@
#include "py/misc.h"
#include "utils.h"
#include "../numpy/fft/fft_tools.h"
#if ULAB_HAS_UTILS_MODULE
enum UTILS_BUFFER_TYPE {
@@ -187,8 +189,41 @@ static mp_obj_t utils_from_uint32_buffer(size_t n_args, const mp_obj_t *pos_args
MP_DEFINE_CONST_FUN_OBJ_KW(utils_from_uint32_buffer_obj, 1, utils_from_uint32_buffer);
#endif
#endif /* ULAB_UTILS_HAS_FROM_INT16_BUFFER | ULAB_UTILS_HAS_FROM_UINT16_BUFFER | ULAB_UTILS_HAS_FROM_INT32_BUFFER | ULAB_UTILS_HAS_FROM_UINT32_BUFFER */
#if ULAB_UTILS_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.welch`` https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html."""
//| ...
//|
mp_obj_t utils_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(utils_spectrogram_obj, 1, 1, utils_spectrogram);
#else
MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(utils_spectrogram_obj, 1, 2, utils_spectrogram);
#endif
#endif /* ULAB_UTILS_HAS_SPECTROGRAM */
static const mp_rom_map_elem_t ulab_utils_globals_table[] = {
{ MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_utils) },
#if ULAB_UTILS_HAS_FROM_INT16_BUFFER
@@ -203,6 +238,9 @@ static const mp_rom_map_elem_t ulab_utils_globals_table[] = {
#if ULAB_UTILS_HAS_FROM_UINT32_BUFFER
{ MP_OBJ_NEW_QSTR(MP_QSTR_from_uint32_buffer), (mp_obj_t)&utils_from_uint32_buffer_obj },
#endif
#if ULAB_UTILS_HAS_SPECTROGRAM
{ MP_OBJ_NEW_QSTR(MP_QSTR_spectrogram), (mp_obj_t)&utils_spectrogram_obj },
#endif
};
static MP_DEFINE_CONST_DICT(mp_module_ulab_utils_globals, ulab_utils_globals_table);
@@ -212,4 +250,4 @@ mp_obj_module_t ulab_utils_module = {
.globals = (mp_obj_dict_t*)&mp_module_ulab_utils_globals,
};
#endif
#endif /* ULAB_HAS_UTILS_MODULE */

42
python/port/mpconfigport.h
View File

@@ -1,6 +1,9 @@
#include <stdint.h>
#include <alloca.h>
// Include helpers when this is not a MicroPython build.
#ifdef EPSILON_VERSION
#include "helpers.h"
#endif
/* MicroPython configuration options
* We're not listing the default options as defined in mpconfig.h */
@@ -16,6 +19,11 @@
* are therefore erased prematurely. */
#define MICROPY_ENABLE_PYSTACK (0)
// Whether to encode None/False/True as immediate objects instead of pointers to
// real objects. Reduces code size by a decent amount without hurting
// performance, for all representations except D on some architectures.
#define MICROPY_OBJ_IMMEDIATE_OBJS 0
// Maximum length of a path in the filesystem
#define MICROPY_ALLOC_PATH_MAX (32)
@@ -139,40 +147,6 @@ typedef long mp_off_t;
#define MP_STATE_PORT MP_STATE_VM
extern const struct _mp_obj_module_t modion_module;
extern const struct _mp_obj_module_t modkandinsky_module;
extern const struct _mp_obj_module_t modmatplotlib_module;
extern const struct _mp_obj_module_t modpyplot_module;
extern const struct _mp_obj_module_t modtime_module;
extern const struct _mp_obj_module_t modos_module;
extern const struct _mp_obj_module_t modturtle_module;
#if !defined(INCLUDE_ULAB)
#define MICROPY_PORT_BUILTIN_MODULES \
{ MP_ROM_QSTR(MP_QSTR_ion), MP_ROM_PTR(&modion_module) }, \
{ MP_ROM_QSTR(MP_QSTR_kandinsky), MP_ROM_PTR(&modkandinsky_module) }, \
{ MP_ROM_QSTR(MP_QSTR_matplotlib), MP_ROM_PTR(&modmatplotlib_module) }, \
{ MP_ROM_QSTR(MP_QSTR_matplotlib_dot_pyplot), MP_ROM_PTR(&modpyplot_module) }, \
{ MP_ROM_QSTR(MP_QSTR_time), MP_ROM_PTR(&modtime_module) }, \
{ MP_ROM_QSTR(MP_QSTR_os), MP_ROM_PTR(&modos_module) }, \
{ MP_ROM_QSTR(MP_QSTR_turtle), MP_ROM_PTR(&modturtle_module) }, \
#else
extern const struct _mp_obj_module_t ulab_user_cmodule;
#define MICROPY_PORT_BUILTIN_MODULES \
{ MP_ROM_QSTR(MP_QSTR_ion), MP_ROM_PTR(&modion_module) }, \
{ MP_ROM_QSTR(MP_QSTR_kandinsky), MP_ROM_PTR(&modkandinsky_module) }, \
{ MP_ROM_QSTR(MP_QSTR_matplotlib), MP_ROM_PTR(&modmatplotlib_module) }, \
{ MP_ROM_QSTR(MP_QSTR_matplotlib_dot_pyplot), MP_ROM_PTR(&modpyplot_module) }, \
{ MP_ROM_QSTR(MP_QSTR_time), MP_ROM_PTR(&modtime_module) }, \
{ MP_ROM_QSTR(MP_QSTR_os), MP_ROM_PTR(&modos_module) }, \
{ MP_ROM_QSTR(MP_QSTR_turtle), MP_ROM_PTR(&modturtle_module) }, \
{ MP_ROM_QSTR(MP_QSTR_ulab), MP_ROM_PTR(&ulab_user_cmodule) }, \
#endif
// Enable setjmp in debug mode. This is to avoid some optimizations done
// specifically for x86_64 using inline assembly, which makes the debug binary