jepler/micropython-ulab
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Compare commits

merge into: jepler:keepdims
Branches Tags
jepler:master
jepler:bincount
jepler:transpose
jepler:modulo
jepler:keepdims-fix
jepler:keepdims
jepler:take
jepler:buffer
jepler:svd
jepler:pid2
jepler:spectrum
jepler:random
jepler:circuitpython9
jepler:logical2
jepler:logical
jepler:vector-out
jepler:bitwise
jepler:bool-slice
jepler:floordiv
jepler:sort_complex-fix
jepler:arange-fix
jepler:vectorize-fix
jepler:ifft-fix
jepler:version4
jepler:doc0546
jepler:doc0361
jepler:legacy
jepler:dtype-mod
jepler:testing
jepler:6.7.3
jepler:6.7.2
jepler:6.6.1
jepler:6.5.5
jepler:6.5.2
jepler:6.0.12
jepler:5.1.1
jepler:4.0.0
jepler:3.3.7
jepler:0.54.6
jepler:0.36.1
jepler:1.7.9
jepler:3.3.1
jepler:3.0.0
jepler:1.7.8
jepler:2.7.1
jepler:2.6.1
jepler:2.6.0
jepler:1.7.2
jepler:2.1.5
jepler:2.1.3
jepler:2.1.2
jepler:1.7.1
jepler:1.7.0
jepler:1.6.1
jepler:1.1.0
jepler:0.54.5
jepler:0.54
jepler:0.51.1
jepler:0.50.2
jepler:0.36.0
jepler:0.24
jepler:0.23
...
pull from: jepler:master
Branches Tags
jepler:bincount
jepler:transpose
jepler:master
jepler:modulo
jepler:keepdims-fix
jepler:keepdims
jepler:take
jepler:buffer
jepler:svd
jepler:pid2
jepler:spectrum
jepler:random
jepler:circuitpython9
jepler:logical2
jepler:logical
jepler:vector-out
jepler:bitwise
jepler:bool-slice
jepler:floordiv
jepler:sort_complex-fix
jepler:arange-fix
jepler:vectorize-fix
jepler:ifft-fix
jepler:version4
jepler:doc0546
jepler:doc0361
jepler:legacy
jepler:dtype-mod
jepler:testing
jepler:6.7.3
jepler:6.7.2
jepler:6.6.1
jepler:6.5.5
jepler:6.5.2
jepler:6.0.12
jepler:5.1.1
jepler:4.0.0
jepler:3.3.7
jepler:0.54.6
jepler:0.36.1
jepler:1.7.9
jepler:3.3.1
jepler:3.0.0
jepler:1.7.8
jepler:2.7.1
jepler:2.6.1
jepler:2.6.0
jepler:1.7.2
jepler:2.1.5
jepler:2.1.3
jepler:2.1.2
jepler:1.7.1
jepler:1.7.0
jepler:1.6.1
jepler:1.1.0
jepler:0.54.5
jepler:0.54
jepler:0.51.1
jepler:0.50.2
jepler:0.36.0
jepler:0.24
jepler:0.23

12 commits
keepdims ... master

Author SHA1 Message Date
Zoltán Vörös
068da5fc96
Modulo (#734) 
* add modulo operator

* fix module loops

* add in-place modulo operator

* update readme

* add test files, update documentation
 2025-08-05 20:40:40 +02:00
Zoltán Vörös
a0999aba79
add modulo operator (#733) 
* add modulo operator

* fix modulo loops

* add in-place modulo operator

* update readme
 2025-08-04 22:55:02 +02:00
Zoltán Vörös
8eb8eaf5a1 update changelog; will push directly into master 2025-06-06 21:47:50 +02:00
Zoltán Vörös
66e9eb3ed3
add ndim property (#725) 2025-06-06 21:46:21 +02:00
Zoltán Vörös
844f6e5e1e
fix type inference 
* fix type inference

* make checks inclusive for negative values of integers

* move to newer version of ubuntu

* fix iterator macros

* fix iterator macros
 2025-06-06 21:29:41 +02:00
Zoltán Vörös
1398a8606f
loadtxt can deal with multi-line comments (#720) 
* loadtxt can deal with multi-line comments

* multiline headers/footers are treated correctly

* add macro to traverse arrays
 2025-06-06 20:58:31 +02:00
Zoltán Vörös
88ef893540
Random fix (#718) 
* fix typo and shape in radnom module

* add random test files

* fix expected file
 2025-05-29 19:16:43 +02:00
Zoltán Vörös
825ec2b143
re-name integration constants (#714) 2025-03-16 20:37:00 +01:00
Brian Whitman
66daa8960a
Adding tulipcc to ulab's users (#711) 2025-02-16 19:36:33 +01:00
Zoltán Vörös
20f7259a47
ix keepdims for min, max, argmin, argmax (#707) 2025-01-26 21:04:20 +01:00
Zoltán Vörös
be15d62632
fix keepdims for std, update docs (#704) 2025-01-19 22:51:02 +01:00
Zoltán Vörös
903506ca9a
Add keepdims keyword argument (#701) 
* add function to deal with keepdims=True

* preliminary keepdims fix

* fux keepdims code

* remove out-commented code
 2025-01-15 10:29:27 +01:00
34 changed files with 1414 additions and 1074 deletions
Show stats Expand all files Collapse all files

4
.github/workflows/build.yml vendored
View file

@ -20,7 +20,7 @@ jobs:
strategy:
matrix:
os:
- ubuntu-20.04
- ubuntu-24.04
- macOS-latest
dims: [1, 2, 3, 4]
runs-on: ${{ matrix.os }}
@ -61,7 +61,7 @@ jobs:
strategy:
matrix:
os:
- ubuntu-20.04
- ubuntu-24.04
- macOS-latest
dims: [1, 2, 3, 4]
runs-on: ${{ matrix.os }}

3
README.md
View file

@ -36,7 +36,7 @@ detected and handled.
## ndarray methods
`ulab` implements `numpy`'s `ndarray` with the `==`, `!=`, `<`, `<=`, `>`, `>=`, `+`, `-`, `/`, `*`, `**`,
`+=`, `-=`, `*=`, `/=`, `**=` binary operators, and the `len`, `~`, `-`, `+`, `abs` unary operators that
`%`, `+=`, `-=`, `*=`, `/=`, `**=`, `%=` binary operators, and the `len`, `~`, `-`, `+`, `abs` unary operators that
operate element-wise. Type-aware `ndarray`s can be initialised from any `micropython` iterable, lists of
iterables via the `array` constructor, or by means of the `arange`, `concatenate`, `diag`, `eye`,
`frombuffer`, `full`, `linspace`, `logspace`, `ones`, or `zeros` functions.
@ -112,6 +112,7 @@ of the user manual.
1. `MaixPy` https://github.com/sipeed/MaixPy
1. `OpenMV` https://github.com/openmv/openmv
1. `pimoroni-pico` https://github.com/pimoroni/pimoroni-pico
1. `Tulip Creative Computer` https://github.com/shorepine/tulipcc
## Compiling

174
code/ndarray.c
View file

@ -563,6 +563,9 @@ ndarray_obj_t *ndarray_new_dense_ndarray(uint8_t ndim, size_t *shape, uint8_t dt
ndarray_obj_t *ndarray_new_ndarray_from_tuple(mp_obj_tuple_t *_shape, uint8_t dtype) {
// creates a dense array from a tuple
// the function should work in the general n-dimensional case
if(_shape->len > ULAB_MAX_DIMS) {
mp_raise_ValueError(MP_ERROR_TEXT("maximum number of dimensions is " MP_STRINGIFY(ULAB_MAX_DIMS)));
}
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
for(size_t i = 0; i < _shape->len; i++) {
shape[ULAB_MAX_DIMS - 1 - i] = mp_obj_get_int(_shape->items[_shape->len - 1 - i]);
@ -583,43 +586,10 @@ void ndarray_copy_array(ndarray_obj_t *source, ndarray_obj_t *target, uint8_t sh
}
#endif
#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 {
ITERATOR_HEAD();
memcpy(tarray, sarray, target->itemsize);
tarray += target->itemsize;
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
ITERATOR_TAIL(source, sarray);
}
ndarray_obj_t *ndarray_new_view(ndarray_obj_t *source, uint8_t ndim, size_t *shape, int32_t *strides, int32_t offset) {
@ -673,20 +643,7 @@ ndarray_obj_t *ndarray_copy_view_convert_type(ndarray_obj_t *source, uint8_t dty
uint8_t complex_size = 2 * sizeof(mp_float_t);
#endif
#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 {
ITERATOR_HEAD()
mp_obj_t item;
#if ULAB_SUPPORTS_COMPLEX
if(source->dtype == NDARRAY_COMPLEX) {
@ -715,27 +672,7 @@ ndarray_obj_t *ndarray_copy_view_convert_type(ndarray_obj_t *source, uint8_t dty
ndarray_set_value(dtype, array, 0, item);
#endif
array += ndarray->itemsize;
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
ITERATOR_TAIL(source, sarray);
return ndarray;
}
@ -762,20 +699,7 @@ mp_obj_t ndarray_byteswap(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_
return MP_OBJ_FROM_PTR(ndarray);
} else {
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 {
ITERATOR_HEAD();
if(self->dtype == NDARRAY_FLOAT) {
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
SWAP(uint8_t, array[0], array[3]);
@ -789,27 +713,7 @@ mp_obj_t ndarray_byteswap(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_
} else {
SWAP(uint8_t, array[0], array[1]);
}
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
ITERATOR_TAIL(ndarray, array);
}
return MP_OBJ_FROM_PTR(ndarray);
}
@ -1438,43 +1342,10 @@ mp_obj_t ndarray_flatten(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_a
uint8_t *array = (uint8_t *)ndarray->array;
if(memcmp(order, "C", 1) == 0) { // C-type ordering
#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 {
ITERATOR_HEAD();
memcpy(array, sarray, self->itemsize);
array += ndarray->strides[ULAB_MAX_DIMS - 1];
sarray += self->strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < self->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
sarray -= self->strides[ULAB_MAX_DIMS - 1] * self->shape[ULAB_MAX_DIMS-1];
sarray += self->strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < self->shape[ULAB_MAX_DIMS - 2]);
#endif
#if ULAB_MAX_DIMS > 2
sarray -= self->strides[ULAB_MAX_DIMS - 2] * self->shape[ULAB_MAX_DIMS-2];
sarray += self->strides[ULAB_MAX_DIMS - 3];
j++;
} while(j < self->shape[ULAB_MAX_DIMS - 3]);
#endif
#if ULAB_MAX_DIMS > 3
sarray -= self->strides[ULAB_MAX_DIMS - 3] * self->shape[ULAB_MAX_DIMS-3];
sarray += self->strides[ULAB_MAX_DIMS - 4];
i++;
} while(i < self->shape[ULAB_MAX_DIMS - 4]);
#endif
ITERATOR_TAIL(self, sarray);
} else { // 'F', Fortran-type ordering
#if ULAB_MAX_DIMS > 3
size_t i = 0;
@ -1527,6 +1398,13 @@ mp_obj_t ndarray_itemsize(mp_obj_t self_in) {
}
#endif
#if NDARRAY_HAS_NDIM
mp_obj_t ndarray_ndim(mp_obj_t self_in) {
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
return MP_OBJ_NEW_SMALL_INT(self->ndim);
}
#endif
#if NDARRAY_HAS_SHAPE
mp_obj_t ndarray_shape(mp_obj_t self_in) {
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
@ -1611,7 +1489,7 @@ ndarray_obj_t *ndarray_from_mp_obj(mp_obj_t obj, uint8_t other_type) {
if(mp_obj_is_int(obj)) {
int32_t ivalue = mp_obj_get_int(obj);
if((ivalue < -32767) || (ivalue > 32767)) {
if((ivalue < -32768) || (ivalue > 65535)) {
// the integer value clearly does not fit the ulab integer types, so move on to float
ndarray = ndarray_new_linear_array(1, NDARRAY_FLOAT);
mp_float_t *array = (mp_float_t *)ndarray->array;
@ -1619,7 +1497,7 @@ ndarray_obj_t *ndarray_from_mp_obj(mp_obj_t obj, uint8_t other_type) {
} else {
uint8_t dtype;
if(ivalue < 0) {
if(ivalue > -128) {
if(ivalue >= -128) {
dtype = NDARRAY_INT8;
} else {
dtype = NDARRAY_INT16;
@ -1770,6 +1648,12 @@ mp_obj_t ndarray_binary_op(mp_binary_op_t _op, mp_obj_t lobj, mp_obj_t robj) {
return ndarray_inplace_ams(lhs, rhs, rstrides, op);
break;
#endif
#if NDARRAY_HAS_INPLACE_MODULO
case MP_BINARY_OP_INPLACE_MODULO:
COMPLEX_DTYPE_NOT_IMPLEMENTED(lhs->dtype);
return ndarray_inplace_modulo(lhs, rhs, rstrides);
break;
#endif
#if NDARRAY_HAS_INPLACE_MULTIPLY
case MP_BINARY_OP_INPLACE_MULTIPLY:
COMPLEX_DTYPE_NOT_IMPLEMENTED(lhs->dtype);
@ -1825,6 +1709,12 @@ mp_obj_t ndarray_binary_op(mp_binary_op_t _op, mp_obj_t lobj, mp_obj_t robj) {
return ndarray_binary_add(lhs, rhs, ndim, shape, lstrides, rstrides);
break;
#endif
#if NDARRAY_HAS_BINARY_OP_MODULO
case MP_BINARY_OP_MODULO:
COMPLEX_DTYPE_NOT_IMPLEMENTED(lhs->dtype);
return ndarray_binary_modulo(lhs, rhs, ndim, shape, lstrides, rstrides);
break;
#endif
#if NDARRAY_HAS_BINARY_OP_MULTIPLY
case MP_BINARY_OP_MULTIPLY:
return ndarray_binary_multiply(lhs, rhs, ndim, shape, lstrides, rstrides);

89
code/ndarray.h
View file

@ -232,6 +232,10 @@ mp_obj_t ndarray_dtype(mp_obj_t );
mp_obj_t ndarray_itemsize(mp_obj_t );
#endif
#if NDARRAY_HAS_NDIM
mp_obj_t ndarray_ndim(mp_obj_t );
#endif
#if NDARRAY_HAS_SIZE
mp_obj_t ndarray_size(mp_obj_t );
#endif
@ -709,4 +713,89 @@ ndarray_obj_t *ndarray_from_mp_obj(mp_obj_t , uint8_t );
#endif /* ULAB_MAX_DIMS == 4 */
#endif /* ULAB_HAS_FUNCTION_ITERATOR */
// iterator macro for traversing arrays over all dimensions
#if ULAB_MAX_DIMS == 1
#define ITERATOR_HEAD()\
size_t _l_ = 0;\
do {
#define ITERATOR_TAIL(_source_, _source_array_)\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 1];\
_l_++;\
} while(_l_ < (_source_)->shape[ULAB_MAX_DIMS - 1]);
#endif /* ULAB_MAX_DIMS == 1 */
#if ULAB_MAX_DIMS == 2
#define ITERATOR_HEAD()\
size_t _k_ = 0;\
do {\
size_t _l_ = 0;\
do {
#define ITERATOR_TAIL(_source_, _source_array_)\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 1];\
_l_++;\
} while(_l_ < (_source_)->shape[ULAB_MAX_DIMS - 1]);\
(_source_array_) -= (_source_)->strides[ULAB_MAX_DIMS - 1] * (_source_)->shape[ULAB_MAX_DIMS - 1];\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 2];\
_k_++;\
} while(_k_ < (_source_)->shape[ULAB_MAX_DIMS - 2]);
#endif /* ULAB_MAX_DIMS == 2 */
#if ULAB_MAX_DIMS == 3
#define ITERATOR_HEAD()\
size_t _j_ = 0;\
do {\
size_t _k_ = 0;\
do {\
size_t _l_ = 0;\
do {
#define ITERATOR_TAIL(_source_, _source_array_)\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 1];\
_l_++;\
} while(_l_ < (_source_)->shape[ULAB_MAX_DIMS - 1]);\
(_source_array_) -= (_source_)->strides[ULAB_MAX_DIMS - 1] * (_source_)->shape[ULAB_MAX_DIMS - 1];\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 2];\
_k_++;\
} while(_k_ < (_source_)->shape[ULAB_MAX_DIMS - 2]);\
(_source_array_) -= (_source_)->strides[ULAB_MAX_DIMS - 2] * (_source_)->shape[ULAB_MAX_DIMS - 2];\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 3];\
_j_++;\
} while(_j_ < (_source_)->shape[ULAB_MAX_DIMS - 3]);
#endif /* ULAB_MAX_DIMS == 3 */
#if ULAB_MAX_DIMS == 4
#define ITERATOR_HEAD()\
size_t _i_ = 0;\
do {\
size_t _j_ = 0;\
do {\
size_t _k_ = 0;\
do {\
size_t _l_ = 0;\
do {
#define ITERATOR_TAIL(_source_, _source_array_)\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 1];\
_l_++;\
} while(_l_ < (_source_)->shape[ULAB_MAX_DIMS - 1]);\
(_source_array_) -= (_source_)->strides[ULAB_MAX_DIMS - 1] * (_source_)->shape[ULAB_MAX_DIMS - 1];\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 2];\
_k_++;\
} while(_k_ < (_source_)->shape[ULAB_MAX_DIMS - 2]);\
(_source_array_) -= (_source_)->strides[ULAB_MAX_DIMS - 2] * (_source_)->shape[ULAB_MAX_DIMS - 2];\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 3];\
_j_++;\
} while(_j_ < (_source_)->shape[ULAB_MAX_DIMS - 3]);\
(_source_array_) -= (_source_)->strides[ULAB_MAX_DIMS - 3] * (_source_)->shape[ULAB_MAX_DIMS - 3];\
(_source_array_) += (_source_)->strides[ULAB_MAX_DIMS - 4];\
_i_++;\
} while(_i_ < (_source_)->shape[ULAB_MAX_DIMS - 4]);
#endif /* ULAB_MAX_DIMS == 4 */
#endif

122
code/ndarray_operators.c
View file

@ -248,6 +248,105 @@ mp_obj_t ndarray_binary_add(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
}
#endif /* NDARRAY_HAS_BINARY_OP_ADD */
#if NDARRAY_HAS_BINARY_OP_MODULO
mp_obj_t ndarray_binary_modulo(ndarray_obj_t *lhs, ndarray_obj_t *rhs,
uint8_t ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {
ndarray_obj_t *results = NULL;
uint8_t *larray = (uint8_t *)lhs->array;
uint8_t *rarray = (uint8_t *)rhs->array;
if(lhs->dtype == NDARRAY_UINT8) {
if(rhs->dtype == NDARRAY_UINT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_UINT8);
BINARY_LOOP(results, uint8_t, uint8_t, uint8_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_INT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_INT16);
BINARY_LOOP(results, int16_t, uint8_t, int8_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_UINT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_UINT16);
BINARY_LOOP(results, uint16_t, uint8_t, uint16_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_INT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_INT16);
BINARY_LOOP(results, int16_t, uint8_t, int16_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_FLOAT) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
MODULO_FLOAT_LOOP(results, mp_float_t, uint8_t, mp_float_t, larray, lstrides, rarray, rstrides);
}
} else if(lhs->dtype == NDARRAY_INT8) {
if(rhs->dtype == NDARRAY_UINT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_INT16);
BINARY_LOOP(results, int16_t, int8_t, uint8_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_INT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_INT8);
BINARY_LOOP(results, int8_t, int8_t, int8_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_UINT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_INT16);
BINARY_LOOP(results, int16_t, int8_t, int16_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_INT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_INT16);
BINARY_LOOP(results, int16_t, int8_t, int16_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_FLOAT) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
MODULO_FLOAT_LOOP(results, mp_float_t, int8_t, mp_float_t, larray, lstrides, rarray, rstrides);
}
} else if(lhs->dtype == NDARRAY_UINT16) {
if(rhs->dtype == NDARRAY_UINT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_UINT8);
BINARY_LOOP(results, uint16_t, uint16_t, uint8_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_INT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
BINARY_LOOP(results, mp_float_t, uint16_t, int8_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_UINT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_UINT16);
BINARY_LOOP(results, uint16_t, uint16_t, uint16_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_INT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
BINARY_LOOP(results, mp_float_t, uint16_t, int16_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_FLOAT) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
MODULO_FLOAT_LOOP(results, mp_float_t, uint16_t, mp_float_t, larray, lstrides, rarray, rstrides);
}
} else if(lhs->dtype == NDARRAY_INT16) {
if(rhs->dtype == NDARRAY_UINT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_INT16);
BINARY_LOOP(results, int16_t, int16_t, uint8_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_INT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_INT16);
BINARY_LOOP(results, int16_t, int16_t, int8_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_UINT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
BINARY_LOOP(results, mp_float_t, int16_t, uint16_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_INT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_INT16);
BINARY_LOOP(results, int16_t, int16_t, int16_t, larray, lstrides, rarray, rstrides, %);
} else if(rhs->dtype == NDARRAY_FLOAT) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
MODULO_FLOAT_LOOP(results, mp_float_t, int16_t, mp_float_t, larray, lstrides, rarray, rstrides);
}
} else if(lhs->dtype == NDARRAY_FLOAT) {
if(rhs->dtype == NDARRAY_UINT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
MODULO_FLOAT_LOOP(results, mp_float_t, mp_float_t, uint8_t, larray, lstrides, rarray, rstrides);
} else if(rhs->dtype == NDARRAY_INT8) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
MODULO_FLOAT_LOOP(results, mp_float_t, mp_float_t, int8_t, larray, lstrides, rarray, rstrides);
} else if(rhs->dtype == NDARRAY_UINT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
MODULO_FLOAT_LOOP(results, mp_float_t, mp_float_t, uint16_t, larray, lstrides, rarray, rstrides);
} else if(rhs->dtype == NDARRAY_INT16) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
MODULO_FLOAT_LOOP(results, mp_float_t, mp_float_t, int16_t, larray, lstrides, rarray, rstrides);
} else if(rhs->dtype == NDARRAY_FLOAT) {
results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
MODULO_FLOAT_LOOP(results, mp_float_t, mp_float_t, mp_float_t, larray, lstrides, rarray, rstrides);
}
}
return MP_OBJ_FROM_PTR(results);
}
#endif /* NDARRAY_HAS_BINARY_OP_MODULO */
#if NDARRAY_HAS_BINARY_OP_MULTIPLY
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) {
@ -1074,6 +1173,29 @@ mp_obj_t ndarray_inplace_ams(ndarray_obj_t *lhs, ndarray_obj_t *rhs, int32_t *rs
}
#endif /* NDARRAY_HAS_INPLACE_ADD || NDARRAY_HAS_INPLACE_MULTIPLY || NDARRAY_HAS_INPLACE_SUBTRACT */
#if NDARRAY_HAS_INPLACE_MODULO
mp_obj_t ndarray_inplace_modulo(ndarray_obj_t *lhs, ndarray_obj_t *rhs, int32_t *rstrides) {
if((lhs->dtype != NDARRAY_FLOAT) && (rhs->dtype == NDARRAY_FLOAT)) {
mp_raise_TypeError(MP_ERROR_TEXT("results cannot be cast to specified type"));
}
if(lhs->dtype == NDARRAY_FLOAT) {
if(rhs->dtype == NDARRAY_UINT8) {
INLINE_MODULO_FLOAT_LOOP(lhs, uint8_t, larray, rarray, rstrides);
} else if(rhs->dtype == NDARRAY_UINT8) {
INLINE_MODULO_FLOAT_LOOP(lhs, int8_t, larray, rarray, rstrides);
} else if(rhs->dtype == NDARRAY_UINT16) {
INLINE_MODULO_FLOAT_LOOP(lhs, uint16_t, larray, rarray, rstrides);
} else if(rhs->dtype == NDARRAY_INT16) {
INLINE_MODULO_FLOAT_LOOP(lhs, int16_t, larray, rarray, rstrides);
} else {
INLINE_MODULO_FLOAT_LOOP(lhs, mp_float_t, larray, rarray, rstrides);
}
}
return MP_OBJ_FROM_PTR(lhs);
}
#endif /* NDARRAY_HAS_INPLACE_MODULO */
#if NDARRAY_HAS_INPLACE_TRUE_DIVIDE
mp_obj_t ndarray_inplace_divide(ndarray_obj_t *lhs, ndarray_obj_t *rhs, int32_t *rstrides) {

175
code/ndarray_operators.h
View file

@ -12,6 +12,7 @@
mp_obj_t ndarray_binary_equality(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *, mp_binary_op_t );
mp_obj_t ndarray_binary_add(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *);
mp_obj_t ndarray_binary_modulo(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *);
mp_obj_t ndarray_binary_multiply(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *);
mp_obj_t ndarray_binary_more(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *, mp_binary_op_t );
mp_obj_t ndarray_binary_power(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *);
@ -21,6 +22,7 @@ mp_obj_t ndarray_binary_logical(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size
mp_obj_t ndarray_binary_floor_divide(ndarray_obj_t *, ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t *);
mp_obj_t ndarray_inplace_ams(ndarray_obj_t *, ndarray_obj_t *, int32_t *, uint8_t );
mp_obj_t ndarray_inplace_modulo(ndarray_obj_t *, ndarray_obj_t *, int32_t *);
mp_obj_t ndarray_inplace_power(ndarray_obj_t *, ndarray_obj_t *, int32_t *);
mp_obj_t ndarray_inplace_divide(ndarray_obj_t *, ndarray_obj_t *, int32_t *);
@ -537,3 +539,176 @@ mp_obj_t ndarray_inplace_divide(ndarray_obj_t *, ndarray_obj_t *, int32_t *);
} while(0)
#endif /* ULAB_MAX_DIMS == 4 */
#define MODULO_FLOAT1(results, array, type_left, type_right, larray, lstrides, rarray, rstrides)\
({\
size_t l = 0;\
do {\
*(array)++ = MICROPY_FLOAT_C_FUN(fmod)(*((type_left *)(larray)), *((type_right *)(rarray)));\
(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
#define MODULO_FLOAT_LOOP(results, type_out, type_left, type_right, larray, lstrides, rarray, rstrides) do {\
type_out *array = (type_out *)(results)->array;\
MODULO_FLOAT1((results), (array), type_left, type_right, (larray), (lstrides), (rarray), (rstrides));\
} while(0)
#endif /* ULAB_MAX_DIMS == 1 */
#if ULAB_MAX_DIMS == 2
#define MODULO_FLOAT_LOOP(results, type_out, type_left, type_right, larray, lstrides, rarray, rstrides) do {\
type_out *array = (type_out *)(results)->array;\
size_t l = 0;\
do {\
MODULO_FLOAT1((results), (array), type_left, 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];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 2]);\
} while(0)
#endif /* ULAB_MAX_DIMS == 2 */
#if ULAB_MAX_DIMS == 3
#define MODULO_FLOAT_LOOP(results, type_out, type_left, type_right, larray, lstrides, rarray, rstrides) do {\
type_out *array = (type_out *)(results)->array;\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
MODULO_FLOAT1((results), (array), type_left, 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];\
l++;\
} while(l < (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];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 3]);\
} while(0)
#endif /* ULAB_MAX_DIMS == 3 */
#if ULAB_MAX_DIMS == 4
#define MODULO_FLOAT_LOOP(results, type_out, type_left, type_right, larray, lstrides, rarray, rstrides) do {\
type_out *array = (type_out *)(results)->array;\
size_t j = 0;\
do {\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
MODULO_FLOAT1((results), (array), type_left, 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];\
l++;\
} while(l < (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];\
k++;\
} while(k < (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];\
j++;\
} while(j < (results)->shape[ULAB_MAX_DIMS - 4]);\
} while(0)
#endif /* ULAB_MAX_DIMS == 4 */
#define INPLACE_MODULO_FLOAT1(results, type_right, larray, rarray, rstrides)\
({\
size_t l = 0;\
do {\
*((mp_float_t *)larray) = MICROPY_FLOAT_C_FUN(fmod)(*((mp_float_t *)(larray)), *((type_right *)(rarray)));\
(larray) += (results)->strides[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 1]);\
})
#if ULAB_MAX_DIMS == 1
#define INPLACE_MODULO_FLOAT_LOOP(results, type_right, larray, rarray, rstrides) do {\
INPLACE_MODULO_FLOAT1((results), type_right, (larray), (rarray), (rstrides));\
} while(0)
#endif /* ULAB_MAX_DIMS == 1 */
#if ULAB_MAX_DIMS == 2
#define INLINE_MODULO_FLOAT_LOOP(results, type_right, larray, rarray, rstrides) do {\
size_t l = 0;\
do {\
INPLACE_MODULO_FLOAT1((results), type_right, (larray), (rarray), (rstrides));\
(larray) -= (results)->strides[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(larray) += (results)->strides[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 2]);\
} while(0)
#endif /* ULAB_MAX_DIMS == 2 */
#if ULAB_MAX_DIMS == 3
#define INLINE_MODULO_FLOAT_LOOP(results, type_right, larray, rarray, rstrides) do {\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
INPLACE_MODULO_FLOAT1((results), type_right, (larray), (rarray), (rstrides));\
(larray) -= (results)->strides[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(larray) += (results)->strides[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 2]);\
(larray) -= (results)->strides[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(larray) += (results)->strides[ULAB_MAX_DIMS - 3];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 3];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 3]);\
} while(0)
#endif /* ULAB_MAX_DIMS == 3 */
#if ULAB_MAX_DIMS == 4
#define INLINE_MODULO_FLOAT_LOOP(results, type_right, larray, rarray, rstrides) do {\
size_t j = 0;\
do {\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
INPLACE_MODULO_FLOAT1((results), type_right, (larray), (rarray), (rstrides));\
(larray) -= (results)->strides[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(larray) += (results)->strides[ULAB_MAX_DIMS - 2];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 1] * (results)->shape[ULAB_MAX_DIMS - 1];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 2];\
l++;\
} while(l < (results)->shape[ULAB_MAX_DIMS - 2]);\
(larray) -= (results)->strides[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(larray) += (results)->strides[ULAB_MAX_DIMS - 3];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 2] * (results)->shape[ULAB_MAX_DIMS - 2];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 3];\
k++;\
} while(k < (results)->shape[ULAB_MAX_DIMS - 3]);\
(larray) -= (results)->strides[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(larray) += (results)->strides[ULAB_MAX_DIMS - 4];\
(rarray) -= (rstrides)[ULAB_MAX_DIMS - 3] * (results)->shape[ULAB_MAX_DIMS - 3];\
(rarray) += (rstrides)[ULAB_MAX_DIMS - 4];\
j++;\
} while(j < (results)->shape[ULAB_MAX_DIMS - 4]);\
} while(0)
#endif /* ULAB_MAX_DIMS == 4 */

7
code/ndarray_properties.c
View file

@ -6,7 +6,7 @@
*
* The MIT License (MIT)
*
* Copyright (c) 2021 Zoltán Vörös
* Copyright (c) 2021-2025 Zoltán Vörös
*
*/
@ -42,6 +42,11 @@ void ndarray_properties_attr(mp_obj_t self_in, qstr attr, mp_obj_t *dest) {
dest[0] = ndarray_itemsize(self_in);
break;
#endif
#if NDARRAY_HAS_NDIM
case MP_QSTR_ndim:
dest[0] = ndarray_ndim(self_in);
break;
#endif
#if NDARRAY_HAS_SHAPE
case MP_QSTR_shape:
dest[0] = ndarray_shape(self_in);

6
code/ndarray_properties.h
View file

@ -7,7 +7,7 @@
* The MIT License (MIT)
*
* Copyright (c) 2020 Jeff Epler for Adafruit Industries
* 2020-2021 Zoltán Vörös
* 2020-2025 Zoltán Vörös
*/
#ifndef _NDARRAY_PROPERTIES_
@ -36,6 +36,10 @@ MP_DEFINE_CONST_FUN_OBJ_1(ndarray_flatiter_make_new_obj, ndarray_flatiter_make_n
MP_DEFINE_CONST_FUN_OBJ_1(ndarray_itemsize_obj, ndarray_itemsize);
#endif
#if NDARRAY_HAS_NDIM
MP_DEFINE_CONST_FUN_OBJ_1(ndarray_ndim_obj, ndarray_ndim);
#endif
#if NDARRAY_HAS_SHAPE
MP_DEFINE_CONST_FUN_OBJ_1(ndarray_shape_obj, ndarray_shape);
#endif

38
code/numpy/compare.c
View file

@ -260,48 +260,14 @@ static mp_obj_t compare_isinf_isfinite(mp_obj_t _x, uint8_t mask) {
}
uint8_t *xarray = (uint8_t *)x->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 {
ITERATOR_HEAD();
mp_float_t value = *(mp_float_t *)xarray;
if(isnan(value)) {
*rarray++ = 0;
} else {
*rarray++ = isinf(value) ? mask : 1 - mask;
}
xarray += x->strides[ULAB_MAX_DIMS - 1];
l++;
} while(l < x->shape[ULAB_MAX_DIMS - 1]);
#if ULAB_MAX_DIMS > 1
xarray -= x->strides[ULAB_MAX_DIMS - 1] * x->shape[ULAB_MAX_DIMS-1];
xarray += x->strides[ULAB_MAX_DIMS - 2];
k++;
} while(k < x->shape[ULAB_MAX_DIMS - 2]);
#endif
#if ULAB_MAX_DIMS > 2
xarray -= x->strides[ULAB_MAX_DIMS - 2] * x->shape[ULAB_MAX_DIMS-2];
xarray += x->strides[ULAB_MAX_DIMS - 3];
j++;
} while(j < x->shape[ULAB_MAX_DIMS - 3]);
#endif
#if ULAB_MAX_DIMS > 3
xarray -= x->strides[ULAB_MAX_DIMS - 3] * x->shape[ULAB_MAX_DIMS-3];
xarray += x->strides[ULAB_MAX_DIMS - 4];
i++;
} while(i < x->shape[ULAB_MAX_DIMS - 4]);
#endif
ITERATOR_TAIL(x, xarray);
return MP_OBJ_FROM_PTR(results);
} else {
mp_raise_TypeError(MP_ERROR_TEXT("wrong input type"));

91
code/numpy/io/io.c
View file

@ -338,7 +338,7 @@ static mp_obj_t io_loadtxt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
buffer[read] = '\0';
offset = buffer;
while(*offset != '\0') {
if(*offset == comment_char) {
while(*offset == comment_char) {
// clear the line till the end, or the buffer's end
while((*offset != '\0')) {
offset++;
@ -425,7 +425,7 @@ static mp_obj_t io_loadtxt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
offset = buffer;
while(*offset != '\0') {
if(*offset == comment_char) {
while(*offset == comment_char) {
// clear the line till the end, or the buffer's end
while((*offset != '\0')) {
offset++;
@ -619,48 +619,14 @@ static mp_obj_t io_save(mp_obj_t file, mp_obj_t ndarray_) {
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 {
ITERATOR_HEAD();
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
ITERATOR_TAIL(ndarray, array);
stream_p->write(stream, buffer, offset, &error);
stream_p->ioctl(stream, MP_STREAM_CLOSE, 0, &error);
@ -751,16 +717,32 @@ static mp_obj_t io_savetxt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
char *buffer = m_new(char, ULAB_IO_BUFFER_SIZE);
int error;
size_t len_comment;
char *comments;
if(mp_obj_is_str(args[5].u_obj)) {
const char *_comments = mp_obj_str_get_data(args[5].u_obj, &len_comment);
comments = (char *)_comments;
} else {
len_comment = 2;
comments = m_new(char, len_comment);
comments[0] = '#';
comments[1] = ' ';
}
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, comments, len_comment, &error);
// We can't write the header in the single chunk, for it might contain line breaks
for(size_t i = 0; i < _len; header++, i++) {
stream_p->write(stream, header, 1, &error);
if((*header == '\n') && (i < _len)) {
stream_p->write(stream, comments, len_comment, &error);
}
}
stream_p->write(stream, "\n", 1, &error);
}
@ -799,16 +781,19 @@ static mp_obj_t io_savetxt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
} while(k < ndarray->shape[ULAB_MAX_DIMS - 2]);
#endif
if(mp_obj_is_str(args[4].u_obj)) {
if(mp_obj_is_str(args[4].u_obj)) { // footer string
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, comments, len_comment, &error);
// We can't write the header in the single chunk, for it might contain line breaks
for(size_t i = 0; i < _len; footer++, i++) {
stream_p->write(stream, footer, 1, &error);
if((*footer == '\n') && (i < _len)) {
stream_p->write(stream, comments, len_comment, &error);
}
}
stream_p->write(stream, "\n", 1, &error);
}

38
code/numpy/numerical.c
View file

@ -274,7 +274,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) {
static mp_obj_t numerical_sum_mean_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, mp_obj_t keepdims, 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);
@ -380,7 +380,7 @@ static mp_obj_t numerical_sum_mean_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t
bool isStd = optype == NUMERICAL_STD ? 1 : 0;
results = ndarray_new_dense_ndarray(_shape_strides.ndim, _shape_strides.shape, NDARRAY_FLOAT);
farray = (mp_float_t *)results->array;
// we can return the 0 array here, if the degrees of freedom is larger than the length of the axis
// we can return the 0 array here, if the degrees of freedom are larger than the length of the axis
if((optype == NUMERICAL_STD) && (_shape_strides.shape[0] <= ddof)) {
return MP_OBJ_FROM_PTR(results);
}
@ -397,11 +397,9 @@ static mp_obj_t numerical_sum_mean_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t
RUN_MEAN_STD(mp_float_t, array, farray, _shape_strides, div, isStd);
}
}
if(results->ndim == 0) { // return a scalar here
return mp_binary_get_val_array(results->dtype, results->array, 0);
}
return MP_OBJ_FROM_PTR(results);
return ulab_tools_restore_dims(ndarray, results, keepdims, _shape_strides);
}
// we should never get to this point
return mp_const_none;
}
#endif
@ -441,7 +439,7 @@ static mp_obj_t numerical_argmin_argmax_iterable(mp_obj_t oin, uint8_t optype) {
}
}
static mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {
static mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t keepdims, mp_obj_t axis, uint8_t optype) {
// TODO: treat the flattened array
if(ndarray->len == 0) {
mp_raise_ValueError(MP_ERROR_TEXT("attempt to get (arg)min/(arg)max of empty sequence"));
@ -521,7 +519,9 @@ static mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t
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;
shape_strides _shape_strides = tools_reduce_axes(ndarray, axis);
uint8_t index = _shape_strides.axis;
ndarray_obj_t *results = NULL;
@ -546,12 +546,9 @@ static mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t
}
m_del(int32_t, strides, ULAB_MAX_DIMS);
if(results->len == 1) {
return mp_binary_get_val_array(results->dtype, results->array, 0);
}
return MP_OBJ_FROM_PTR(results);
return ulab_tools_restore_dims(ndarray, results, keepdims, _shape_strides);
}
// we should never get to this point
return mp_const_none;
}
#endif
@ -560,6 +557,7 @@ static mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_m
static const mp_arg_t allowed_args[] = {
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_NONE} } ,
{ MP_QSTR_axis, MP_ARG_OBJ, { .u_rom_obj = MP_ROM_NONE } },
{ MP_QSTR_keepdims, MP_ARG_OBJ, { .u_rom_obj = MP_ROM_FALSE } },
};
mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
@ -567,6 +565,8 @@ static mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_m
mp_obj_t oin = args[0].u_obj;
mp_obj_t axis = args[1].u_obj;
mp_obj_t keepdims = args[2].u_obj;
if((axis != mp_const_none) && (!mp_obj_is_int(axis))) {
mp_raise_TypeError(MP_ERROR_TEXT("axis must be None, or an integer"));
}
@ -598,11 +598,12 @@ static mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_m
case NUMERICAL_ARGMIN:
case NUMERICAL_ARGMAX:
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
return numerical_argmin_argmax_ndarray(ndarray, axis, optype);
case NUMERICAL_SUM:
return numerical_argmin_argmax_ndarray(ndarray, keepdims, axis, optype);
case NUMERICAL_MEAN:
case NUMERICAL_STD:
case NUMERICAL_SUM:
COMPLEX_DTYPE_NOT_IMPLEMENTED(ndarray->dtype)
return numerical_sum_mean_std_ndarray(ndarray, axis, optype, 0);
return numerical_sum_mean_std_ndarray(ndarray, axis, keepdims, optype, 0);
default:
mp_raise_NotImplementedError(MP_ERROR_TEXT("operation is not implemented on ndarrays"));
}
@ -1385,6 +1386,7 @@ mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_arg
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } } ,
{ MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
{ MP_QSTR_ddof, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} },
{ MP_QSTR_keepdims, MP_ARG_OBJ, { .u_rom_obj = MP_ROM_FALSE } },
};
mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
@ -1393,6 +1395,8 @@ mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_arg
mp_obj_t oin = args[0].u_obj;
mp_obj_t axis = args[1].u_obj;
size_t ddof = args[2].u_int;
mp_obj_t keepdims = args[3].u_obj;
if((axis != mp_const_none) && (mp_obj_get_int(axis) != 0) && (mp_obj_get_int(axis) != 1)) {
// this seems to pass with False, and True...
mp_raise_ValueError(MP_ERROR_TEXT("axis must be None, or an integer"));
@ -1401,7 +1405,7 @@ mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_arg
return numerical_sum_mean_std_iterable(oin, NUMERICAL_STD, ddof);
} else if(mp_obj_is_type(oin, &ulab_ndarray_type)) {
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);
return numerical_sum_mean_std_ndarray(ndarray, axis, NUMERICAL_STD, ddof);
return numerical_sum_mean_std_ndarray(ndarray, axis, keepdims, NUMERICAL_STD, ddof);
} else {
mp_raise_TypeError(MP_ERROR_TEXT("input must be tuple, list, range, or ndarray"));
}

130
code/numpy/numerical.h
View file

@ -57,35 +57,9 @@
(rarray) += (results)->itemsize;\
})
// The mean could be calculated by simply dividing the sum by
// the number of elements, but that method is numerically unstable
#define RUN_MEAN1(type, array, rarray, ss)\
({\
mp_float_t M = 0.0;\
for(size_t i=0; i < (ss).shape[0]; i++) {\
mp_float_t value = (mp_float_t)(*(type *)(array));\
M = M + (value - M) / (mp_float_t)(i+1);\
(array) += (ss).strides[0];\
}\
*(rarray)++ = M;\
})
// Instead of the straightforward implementation of the definition,
// we take the numerically stable Welford algorithm here
// https://www.johndcook.com/blog/2008/09/26/comparing-three-methods-of-computing-standard-deviation/
#define RUN_STD1(type, array, rarray, ss, div)\
({\
mp_float_t M = 0.0, m = 0.0, S = 0.0;\
for(size_t i=0; i < (ss).shape[0]; i++) {\
mp_float_t value = (mp_float_t)(*(type *)(array));\
m = M + (value - M) / (mp_float_t)(i+1);\
S = S + (value - M) * (value - m);\
M = m;\
(array) += (ss).strides[0];\
}\
*(rarray)++ = MICROPY_FLOAT_C_FUN(sqrt)(S / (div));\
})
#define RUN_MEAN_STD1(type, array, rarray, ss, div, isStd)\
({\
mp_float_t M = 0.0, m = 0.0, S = 0.0;\
@ -193,14 +167,6 @@
RUN_SUM1(type, (array), (results), (rarray), (ss));\
} while(0)
#define RUN_MEAN(type, array, rarray, ss) do {\
RUN_MEAN1(type, (array), (rarray), (ss));\
} while(0)
#define RUN_STD(type, array, rarray, ss, div) do {\
RUN_STD1(type, (array), (results), (rarray), (ss), (div));\
} while(0)
#define RUN_MEAN_STD(type, array, rarray, ss, div, isStd) do {\
RUN_MEAN_STD1(type, (array), (rarray), (ss), (div), (isStd));\
} while(0)
@ -234,26 +200,6 @@
} while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\
} while(0)
#define RUN_MEAN(type, array, rarray, ss) do {\
size_t l = 0;\
do {\
RUN_MEAN1(type, (array), (rarray), (ss));\
(array) -= (ss).strides[0] * (ss).shape[0];\
(array) += (ss).strides[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\
} while(0)
#define RUN_STD(type, array, rarray, ss, div) do {\
size_t l = 0;\
do {\
RUN_STD1(type, (array), (rarray), (ss), (div));\
(array) -= (ss).strides[0] * (ss).shape[0];\
(array) += (ss).strides[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\
} while(0)
#define RUN_MEAN_STD(type, array, rarray, ss, div, isStd) do {\
size_t l = 0;\
do {\
@ -325,38 +271,6 @@
} while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\
} while(0)
#define RUN_MEAN(type, array, rarray, ss) do {\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
RUN_MEAN1(type, (array), (rarray), (ss));\
(array) -= (ss).strides[0] * (ss).shape[0];\
(array) += (ss).strides[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\
(array) -= (ss).strides[ULAB_MAX_DIMS - 1] * (ss).shape[ULAB_MAX_DIMS - 1];\
(array) += (ss).strides[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\
} while(0)
#define RUN_STD(type, array, rarray, ss, div) do {\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
RUN_STD1(type, (array), (rarray), (ss), (div));\
(array) -= (ss).strides[0] * (ss).shape[0];\
(array) += (ss).strides[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\
(array) -= (ss).strides[ULAB_MAX_DIMS - 1] * (ss).shape[ULAB_MAX_DIMS - 1];\
(array) += (ss).strides[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\
} while(0)
#define RUN_MEAN_STD(type, array, rarray, ss, div, isStd) do {\
size_t k = 0;\
do {\
@ -467,50 +381,6 @@
} while(j < (ss).shape[ULAB_MAX_DIMS - 3]);\
} while(0)
#define RUN_MEAN(type, array, rarray, ss) do {\
size_t j = 0;\
do {\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
RUN_MEAN1(type, (array), (rarray), (ss));\
(array) -= (ss).strides[0] * (ss).shape[0];\
(array) += (ss).strides[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\
(array) -= (ss).strides[ULAB_MAX_DIMS - 1] * (ss).shape[ULAB_MAX_DIMS - 1];\
(array) += (ss).strides[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\
(array) -= (ss).strides[ULAB_MAX_DIMS - 2] * (ss).shape[ULAB_MAX_DIMS - 2];\
(array) += (ss).strides[ULAB_MAX_DIMS - 3];\
j++;\
} while(j < (ss).shape[ULAB_MAX_DIMS - 3]);\
} while(0)
#define RUN_STD(type, array, rarray, ss, div) do {\
size_t j = 0;\
do {\
size_t k = 0;\
do {\
size_t l = 0;\
do {\
RUN_STD1(type, (array), (rarray), (ss), (div));\
(array) -= (ss).strides[0] * (ss).shape[0];\
(array) += (ss).strides[ULAB_MAX_DIMS - 1];\
l++;\
} while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\
(array) -= (ss).strides[ULAB_MAX_DIMS - 1] * (ss).shape[ULAB_MAX_DIMS - 1];\
(array) += (ss).strides[ULAB_MAX_DIMS - 2];\
k++;\
} while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\
(array) -= (ss).strides[ULAB_MAX_DIMS - 2] * (ss).shape[ULAB_MAX_DIMS - 2];\
(array) += (ss).strides[ULAB_MAX_DIMS - 3];\
j++;\
} while(j < (ss).shape[ULAB_MAX_DIMS - 3]);\
} while(0)
#define RUN_MEAN_STD(type, array, rarray, ss, div, isStd) do {\
size_t j = 0;\
do {\

37
code/numpy/poly.c
View file

@ -197,42 +197,9 @@ mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {
// 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;
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 {
ITERATOR_HEAD();
*array++ = poly_eval(func(sarray), p, plen);
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
ITERATOR_TAIL(source, sarray);
} else {
// o_x had better be a one-dimensional standard iterable
ndarray = ndarray_new_linear_array(mp_obj_get_int(mp_obj_len_maybe(o_x)), NDARRAY_FLOAT);

33
code/numpy/random/random.c
View file

@ -57,7 +57,7 @@ const mp_obj_type_t random_generator_type = {
void random_generator_print(const mp_print_t *print, mp_obj_t self_in, mp_print_kind_t kind) {
(void)kind;
random_generator_obj_t *self = MP_OBJ_TO_PTR(self_in);
mp_printf(MP_PYTHON_PRINTER, "Gnerator() at 0x%p", self);
mp_printf(MP_PYTHON_PRINTER, "Generator() at 0x%p", self);
}
mp_obj_t random_generator_make_new(const mp_obj_type_t *type, size_t n_args, size_t n_kw, const mp_obj_t *args) {
@ -149,12 +149,9 @@ static mp_obj_t random_normal(size_t n_args, const mp_obj_t *pos_args, mp_map_t
ndarray = ndarray_new_linear_array((size_t)mp_obj_get_int(size), NDARRAY_FLOAT);
} else if(mp_obj_is_type(size, &mp_type_tuple)) {
mp_obj_tuple_t *_shape = MP_OBJ_TO_PTR(size);
if(_shape->len > ULAB_MAX_DIMS) {
mp_raise_ValueError(MP_ERROR_TEXT("maximum number of dimensions is " MP_STRINGIFY(ULAB_MAX_DIMS)));
}
ndarray = ndarray_new_ndarray_from_tuple(_shape, NDARRAY_FLOAT);
} else { // input type not supported
mp_raise_TypeError(MP_ERROR_TEXT("shape must be None, and integer or a tuple of integers"));
mp_raise_TypeError(MP_ERROR_TEXT("shape must be None, an integer or a tuple of integers"));
}
} else {
// return single value
@ -221,27 +218,16 @@ static mp_obj_t random_random(size_t n_args, const mp_obj_t *pos_args, mp_map_t
mp_obj_t out = args[2].u_obj;
ndarray_obj_t *ndarray = NULL;
size_t *shape = m_new(size_t, ULAB_MAX_DIMS);
uint8_t ndim = 1;
size_t *shape = m_new0(size_t, ULAB_MAX_DIMS);
if(size != mp_const_none) {
if(mp_obj_is_int(size)) {
shape[ULAB_MAX_DIMS - 1] = (size_t)mp_obj_get_int(size);
} else if(mp_obj_is_type(size, &mp_type_tuple)) {
mp_obj_tuple_t *_shape = MP_OBJ_TO_PTR(size);
if(_shape->len > ULAB_MAX_DIMS) {
mp_raise_ValueError(MP_ERROR_TEXT("maximum number of dimensions is " MP_STRINGIFY(ULAB_MAX_DIMS)));
}
ndim = _shape->len;
for(size_t i = 0; i < ULAB_MAX_DIMS; i++) {
if(i >= ndim) {
shape[ULAB_MAX_DIMS - 1 - i] = 0;
} else {
shape[ULAB_MAX_DIMS - 1 - i] = mp_obj_get_int(_shape->items[i]);
}
}
ndarray = ndarray_new_ndarray_from_tuple(_shape, NDARRAY_FLOAT);
} else { // input type not supported
mp_raise_TypeError(MP_ERROR_TEXT("shape must be None, and integer or a tuple of integers"));
mp_raise_TypeError(MP_ERROR_TEXT("shape must be None, an integer or a tuple of integers"));
}
}
@ -267,7 +253,8 @@ static mp_obj_t random_random(size_t n_args, const mp_obj_t *pos_args, mp_map_t
}
} else { // out == None
if(size != mp_const_none) {
ndarray = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
mp_obj_tuple_t *_shape = MP_OBJ_TO_PTR(size);
ndarray = ndarray_new_ndarray_from_tuple(_shape, NDARRAY_FLOAT);
} else {
// return single value
mp_float_t value;
@ -336,13 +323,9 @@ static mp_obj_t random_uniform(size_t n_args, const mp_obj_t *pos_args, mp_map_t
return mp_obj_new_float(value);
} else if(mp_obj_is_type(size, &mp_type_tuple)) {
mp_obj_tuple_t *_shape = MP_OBJ_TO_PTR(size);
// TODO: this could be reduced, if the inspection was in the ndarray_new_ndarray_from_tuple function
if(_shape->len > ULAB_MAX_DIMS) {
mp_raise_ValueError(MP_ERROR_TEXT("maximum number of dimensions is " MP_STRINGIFY(ULAB_MAX_DIMS)));
}
ndarray = ndarray_new_ndarray_from_tuple(_shape, NDARRAY_FLOAT);
} else { // input type not supported
mp_raise_TypeError(MP_ERROR_TEXT("shape must be None, and integer or a tuple of integers"));
mp_raise_TypeError(MP_ERROR_TEXT("shape must be None, an integer or a tuple of integers"));
}
mp_float_t *array = (mp_float_t *)ndarray->array;

262
code/numpy/vector.c
View file

@ -91,43 +91,10 @@ static mp_obj_t vector_generic_vector(size_t n_args, const mp_obj_t *pos_args, m
mp_float_t (*func)(void *) = ndarray_get_float_function(source->dtype);
#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 {
ITERATOR_HEAD();
mp_float_t value = func(sarray);
*tarray++ = f(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 */
ITERATOR_TAIL(source, sarray);
#else
if(source->dtype == NDARRAY_UINT8) {
ITERATE_VECTOR(uint8_t, target, tarray, tstrides, source, sarray);
@ -171,43 +138,10 @@ static mp_obj_t vector_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)
mp_float_t (*func)(void *) = ndarray_get_float_function(source->dtype);
#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 {
ITERATOR_HEAD();
mp_float_t value = func(sarray);
*array++ = f(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 */
ITERATOR_TAIL(source, sarray);
#else
if(source->dtype == NDARRAY_UINT8) {
ITERATE_VECTOR(uint8_t, array, source, sarray);
@ -327,43 +261,11 @@ mp_obj_t vector_around(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_arg
mp_float_t (*func)(void *) = ndarray_get_float_function(source->dtype);
#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 {
ITERATOR_HEAD();
mp_float_t f = func(sarray);
*narray++ = MICROPY_FLOAT_C_FUN(round)(f * mul) / mul;
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
ITERATOR_TAIL(source, sarray);
return MP_OBJ_FROM_PTR(ndarray);
}
@ -631,46 +533,13 @@ static mp_obj_t vector_exp(mp_obj_t o_in) {
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 {
ITERATOR_HEAD();
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 */
ITERATOR_TAIL(source, sarray);
return MP_OBJ_FROM_PTR(ndarray);
}
}
@ -921,20 +790,7 @@ mp_obj_t vector_sqrt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args)
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 {
ITERATOR_HEAD();
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);
@ -942,47 +798,15 @@ mp_obj_t vector_sqrt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args)
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 */
ITERATOR_TAIL(source, sarray);
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 {
ITERATOR_HEAD();
mp_float_t value = *(mp_float_t *)sarray;
if(value >= MICROPY_FLOAT_CONST(0.0)) {
*array++ = MICROPY_FLOAT_C_FUN(sqrt)(value);
@ -991,27 +815,8 @@ mp_obj_t vector_sqrt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args)
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 */
ITERATOR_TAIL(source, sarray);
return MP_OBJ_FROM_PTR(ndarray);
} else {
mp_raise_TypeError(MP_ERROR_TEXT("input dtype must be float or complex"));
@ -1071,45 +876,12 @@ static mp_obj_t vector_vectorized_function_call(mp_obj_t self_in, size_t n_args,
uint8_t *sarray = (uint8_t *)source->array;
uint8_t *narray = (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 {
ITERATOR_HEAD();
avalue[0] = mp_binary_get_val_array(source->dtype, sarray, 0);
fvalue = MP_OBJ_TYPE_GET_SLOT(self->type, call)(self->fun, 1, 0, avalue);
ndarray_set_value(self->otypes, narray, 0, fvalue);
sarray += source->strides[ULAB_MAX_DIMS - 1];
narray += ndarray->itemsize;
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 */
ITERATOR_TAIL(source, sarray);
return MP_OBJ_FROM_PTR(ndarray);
} else if(mp_obj_is_type(args[0], &mp_type_tuple) || mp_obj_is_type(args[0], &mp_type_list) ||

102
code/scipy/integrate/integrate.c
View file

@ -41,21 +41,21 @@
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
ULAB_DEFINE_FLOAT_CONST(etolerance, MICROPY_FLOAT_CONST(1e-14), 0x283424dcUL, 0x3e901b2b29a4692bULL);
#define MACHEPS MICROPY_FLOAT_CONST(1e-17)
#define ULAB_MACHEPS MICROPY_FLOAT_CONST(1e-17)
#else
ULAB_DEFINE_FLOAT_CONST(etolerance, MICROPY_FLOAT_CONST(1e-8), 0x358637cfUL, 0x3e7010c6f7d42d18ULL);
#define MACHEPS MICROPY_FLOAT_CONST(1e-8)
#define ULAB_MACHEPS MICROPY_FLOAT_CONST(1e-8)
#endif
#define ZERO MICROPY_FLOAT_CONST(0.0)
#define POINT_TWO_FIVE MICROPY_FLOAT_CONST(0.25)
#define ONE MICROPY_FLOAT_CONST(1.0)
#define TWO MICROPY_FLOAT_CONST(2.0)
#define FOUR MICROPY_FLOAT_CONST(4.0)
#define SIX MICROPY_FLOAT_CONST(6.0)
#define TEN MICROPY_FLOAT_CONST(10.0)
#define FIFTEEN MICROPY_FLOAT_CONST(15.0)
#define EPS_5 MICROPY_FLOAT_CONST(1e-5)
#define ULAB_ZERO MICROPY_FLOAT_CONST(0.0)
#define ULAB_POINT_TWO_FIVE MICROPY_FLOAT_CONST(0.25)
#define ULAB_ONE MICROPY_FLOAT_CONST(1.0)
#define ULAB_TWO MICROPY_FLOAT_CONST(2.0)
#define ULAB_FOUR MICROPY_FLOAT_CONST(4.0)
#define ULAB_SIX MICROPY_FLOAT_CONST(6.0)
#define ULAB_TEN MICROPY_FLOAT_CONST(10.0)
#define ULAB_FIFTEEN MICROPY_FLOAT_CONST(15.0)
#define ULAB_EPSILON_5 MICROPY_FLOAT_CONST(1e-5)
static mp_float_t integrate_python_call(const mp_obj_type_t *type, mp_obj_t fun, mp_float_t x, mp_obj_t *fargs, uint8_t nparams) {
@ -68,7 +68,7 @@ static mp_float_t integrate_python_call(const mp_obj_type_t *type, mp_obj_t fun,
// sign helper function
int sign(mp_float_t x) {
if (x >= ZERO)
if (x >= ULAB_ZERO)
return 1;
else
return -1;
@ -85,7 +85,7 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_
mp_obj_t fargs[1];
mp_float_t h2 = integrate_python_call(type, fun, a + d/2, fargs, 0) - integrate_python_call(type, fun, (a + d*2)*4, fargs, 0);
int i = 1, j = 32; // j=32 is optimal to find r
if (isfinite(h2) && MICROPY_FLOAT_C_FUN(fabs)(h2) > EPS_5) { // if |h2| > 2^-16
if (isfinite(h2) && MICROPY_FLOAT_C_FUN(fabs)(h2) > ULAB_EPSILON_5) { // if |h2| > 2^-16
mp_float_t r, fl, fr, h, s = 0, lfl, lfr, lr = 2;
do { // find max j such that fl and fr are finite
j /= 2;
@ -118,8 +118,8 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_
if (s > eps) { // if sum of |h| > eps
h = lfl - lfr; // use last fl and fr before the sign change
r = lr; // use last r before the sign change
if (h != ZERO) // if last diff != 0, back up r by one step
r /= TWO;
if (h != ULAB_ZERO) // if last diff != 0, back up r by one step
r /= ULAB_TWO;
if (MICROPY_FLOAT_C_FUN(fabs)(lfl) < MICROPY_FLOAT_C_FUN(fabs)(lfr))
d /= r; // move d closer to the finite endpoint
else
@ -135,8 +135,8 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_
mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint16_t n, mp_float_t eps, mp_float_t *e) {
const mp_obj_type_t *type = mp_obj_get_type(fun);
mp_obj_t fargs[1];
const mp_float_t tol = TEN*eps;
mp_float_t c = ZERO, d = ONE, s, sign = ONE, v, h = TWO;
const mp_float_t tol = ULAB_TEN * eps;
mp_float_t c = ULAB_ZERO, d = ULAB_ONE, s, sign = ULAB_ONE, v, h = ULAB_TWO;
int k = 0, mode = 0; // Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
if (b < a) { // swap bounds
v = b;
@ -145,8 +145,8 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
sign = -1;
}
if (isfinite(a) && isfinite(b)) {
c = (a+b)/TWO;
d = (b-a)/TWO;
c = (a+b) / ULAB_TWO;
d = (b-a) / ULAB_TWO;
v = c;
}
else if (isfinite(a)) {
@ -165,20 +165,20 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
}
else {
mode = 2; // Sinh-Sinh
v = ZERO;
v = ULAB_ZERO;
}
s = integrate_python_call(type, fun, v, fargs, 0);
do {
mp_float_t p = ZERO, q, fp = ZERO, fm = ZERO, t, eh;
h /= TWO;
mp_float_t p = ULAB_ZERO, q, fp = ULAB_ZERO, fm = ULAB_ZERO, t, eh;
h /= ULAB_TWO;
t = eh = MICROPY_FLOAT_C_FUN(exp)(h);
if (k > ZERO)
if (k > ULAB_ZERO)
eh *= eh;
if (mode == 0) { // Tanh-Sinh
do {
mp_float_t u = MICROPY_FLOAT_C_FUN(exp)(ONE / t - t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
mp_float_t r = TWO * u / (ONE + u); // = 1 - tanh(sinh(j*h))
mp_float_t w = (t + ONE / t) * r / (ONE + u); // = cosh(j*h)/cosh(sinh(j*h))^2
mp_float_t u = MICROPY_FLOAT_C_FUN(exp)(ULAB_ONE / t - t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
mp_float_t r = ULAB_TWO * u / (ULAB_ONE + u); // = 1 - tanh(sinh(j*h))
mp_float_t w = (t + ULAB_ONE / t) * r / (ULAB_ONE + u); // = cosh(j*h)/cosh(sinh(j*h))^2
mp_float_t x = d*r;
if (a+x > a) { // if too close to a then reuse previous fp
mp_float_t y = integrate_python_call(type, fun, a+x, fargs, 0);
@ -196,11 +196,11 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
} while (MICROPY_FLOAT_C_FUN(fabs)(q) > eps*MICROPY_FLOAT_C_FUN(fabs)(p));
}
else {
t /= TWO;
t /= ULAB_TWO;
do {
mp_float_t r = MICROPY_FLOAT_C_FUN(exp)(t - POINT_TWO_FIVE / t); // = exp(sinh(j*h))
mp_float_t r = MICROPY_FLOAT_C_FUN(exp)(t - ULAB_POINT_TWO_FIVE / t); // = exp(sinh(j*h))
mp_float_t x, y, w = r;
q = ZERO;
q = ULAB_ZERO;
if (mode == 1) { // Exp-Sinh
x = c + d/r;
if (x == c) // if x hit the finite endpoint then break
@ -210,8 +210,8 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
q += y/w;
}
else { // Sinh-Sinh
r = (r - ONE / r) / TWO; // = sinh(sinh(j*h))
w = (w + ONE / w) / TWO; // = cosh(sinh(j*h))
r = (r - ULAB_ONE / r) / ULAB_TWO; // = sinh(sinh(j*h))
w = (w + ULAB_ONE / w) / ULAB_TWO; // = cosh(sinh(j*h))
x = c - d*r;
y = integrate_python_call(type, fun, x, fargs, 0);
if (isfinite(y)) // if f(x) is finite, add to local sum
@ -221,7 +221,7 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
y = integrate_python_call(type, fun, x, fargs, 0);
if (isfinite(y)) // if f(x) is finite, add to local sum
q += y*w;
q *= t + POINT_TWO_FIVE / t; // q *= cosh(j*h)
q *= t + ULAB_POINT_TWO_FIVE / t; // q *= cosh(j*h)
p += q;
t *= eh;
} while (MICROPY_FLOAT_C_FUN(fabs)(q) > eps*MICROPY_FLOAT_C_FUN(fabs)(p));
@ -319,12 +319,12 @@ mp_float_t qromb(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint
for (i = 1; i < n; ++i) {
unsigned long long k = 1UL << i;
unsigned long long s = 1;
mp_float_t sum = ZERO;
mp_float_t sum = ULAB_ZERO;
mp_float_t *Rt;
h /= TWO;
h /= ULAB_TWO;
for (j = 1; j < k; j += 2)
sum += integrate_python_call(type, fun, a+j*h, fargs, 0);
Ru[0] = h*sum + Ro[0] / TWO;
Ru[0] = h*sum + Ro[0] / ULAB_TWO;
for (j = 1; j <= i; ++j) {
s <<= 2;
Ru[j] = (s*Ru[j-1] - Ro[j-1])/(s-1);
@ -408,17 +408,17 @@ mp_float_t as(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, mp_floa
mp_float_t fb, mp_float_t v, mp_float_t eps, int n, mp_float_t t) {
const mp_obj_type_t *type = mp_obj_get_type(fun);
mp_obj_t fargs[1];
mp_float_t h = (b-a) / TWO;
mp_float_t f1 = integrate_python_call(type, fun, a + h / TWO, fargs, 0);
mp_float_t f2 = integrate_python_call(type, fun, b - h / TWO, fargs, 0);
mp_float_t sl = h*(fa + FOUR * f1 + fm) / SIX;
mp_float_t sr = h*(fm + FOUR * f2 + fb) / SIX;
mp_float_t h = (b-a) / ULAB_TWO;
mp_float_t f1 = integrate_python_call(type, fun, a + h / ULAB_TWO, fargs, 0);
mp_float_t f2 = integrate_python_call(type, fun, b - h / ULAB_TWO, fargs, 0);
mp_float_t sl = h*(fa + ULAB_FOUR * f1 + fm) / ULAB_SIX;
mp_float_t sr = h*(fm + ULAB_FOUR * f2 + fb) / ULAB_SIX;
mp_float_t s = sl+sr;
mp_float_t d = (s-v) / FIFTEEN;
mp_float_t d = (s-v) / ULAB_FIFTEEN;
mp_float_t m = a+h;
if (n <= 0 || MICROPY_FLOAT_C_FUN(fabs)(d) < eps)
return t + s + d; // note: fabs(d) can be used as error estimate
eps /= TWO;
eps /= ULAB_TWO;
--n;
t = as(fun, a, m, fa, f1, fm, sl, eps, n, t);
return as(fun, m, b, fm, f2, fb, sr, eps, n, t);
@ -430,7 +430,7 @@ mp_float_t qasi(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n
mp_float_t fa = integrate_python_call(type, fun, a, fargs, 0);
mp_float_t fm = integrate_python_call(type, fun, (a+b)/2, fargs, 0);
mp_float_t fb = integrate_python_call(type, fun, b, fargs, 0);
mp_float_t v = (fa + FOUR * fm + fb) * (b-a) / SIX;
mp_float_t v = (fa + ULAB_FOUR * fm + fb) * (b-a) / ULAB_SIX;
return as(fun, a, b, fa, fm, fb, v, eps, n, 0);
}
@ -562,8 +562,8 @@ mp_float_t gk(mp_float_t (*fun)(mp_float_t), mp_float_t c, mp_float_t d, mp_floa
};
const mp_obj_type_t *type = mp_obj_get_type(fun);
mp_obj_t fargs[1];
mp_float_t p = ZERO; // kronrod quadrature sum
mp_float_t q = ZERO; // gauss quadrature sum
mp_float_t p = ULAB_ZERO; // kronrod quadrature sum
mp_float_t q = ULAB_ZERO; // gauss quadrature sum
mp_float_t fp, fm;
mp_float_t e;
int i;
@ -581,24 +581,24 @@ mp_float_t gk(mp_float_t (*fun)(mp_float_t), mp_float_t c, mp_float_t d, mp_floa
p += (fp + fm) * weights[i];
}
*err = MICROPY_FLOAT_C_FUN(fabs)(p - q);
e = MICROPY_FLOAT_C_FUN(fabs)(2 * p * MACHEPS); // optional, to take 1e-17 MachEps prec. into account
e = MICROPY_FLOAT_C_FUN(fabs)(2 * p * ULAB_MACHEPS); // optional, to take 1e-17 MachEps prec. into account
if (*err < e)
*err = e;
return p;
}
mp_float_t qakro(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n, mp_float_t tol, mp_float_t eps, mp_float_t *err) {
mp_float_t c = (a+b) / TWO;
mp_float_t d = (b-a) / TWO;
mp_float_t c = (a+b) / ULAB_TWO;
mp_float_t d = (b-a) / ULAB_TWO;
mp_float_t e;
mp_float_t r = gk(fun, c, d, &e);
mp_float_t s = d*r;
mp_float_t t = MICROPY_FLOAT_C_FUN(fabs)(s*tol);
if (tol == ZERO)
if (tol == ULAB_ZERO)
tol = t;
if (n > 0 && t < e && tol < e) {
s = qakro(fun, a, c, n-1, t / TWO, eps, err);
s += qakro(fun, c, b, n-1, t / TWO, eps, &e);
s = qakro(fun, a, c, n-1, t / ULAB_TWO, eps, err);
s += qakro(fun, c, b, n-1, t / ULAB_TWO, eps, &e);
*err += e;
return s;
}

2
code/ulab.c
View file

@ -33,7 +33,7 @@
#include "user/user.h"
#include "utils/utils.h"
#define ULAB_VERSION 6.7.0
#define ULAB_VERSION 6.9.0
#define xstr(s) str(s)
#define str(s) #s

12
code/ulab.h
View file

@ -117,6 +117,10 @@
#define NDARRAY_HAS_BINARY_OP_LESS_EQUAL (1)
#endif
#ifndef NDARRAY_HAS_BINARY_OP_MODULO
#define NDARRAY_HAS_BINARY_OP_MODULO (1)
#endif
#ifndef NDARRAY_HAS_BINARY_OP_MORE
#define NDARRAY_HAS_BINARY_OP_MORE (1)
#endif
@ -161,6 +165,10 @@
#define NDARRAY_HAS_INPLACE_ADD (1)
#endif
#ifndef NDARRAY_HAS_INPLACE_MODULO
#define NDARRAY_HAS_INPLACE_MODU (1)
#endif
#ifndef NDARRAY_HAS_INPLACE_MULTIPLY
#define NDARRAY_HAS_INPLACE_MULTIPLY (1)
#endif
@ -245,6 +253,10 @@
#define NDARRAY_HAS_ITEMSIZE (1)
#endif
#ifndef NDARRAY_HAS_NDIM
#define NDARRAY_HAS_NDIM (1)
#endif
#ifndef NDARRAY_HAS_RESHAPE
#define NDARRAY_HAS_RESHAPE (1)
#endif

76
code/ulab_tools.c
View file

@ -162,6 +162,15 @@ void *ndarray_set_float_function(uint8_t dtype) {
}
#endif /* NDARRAY_BINARY_USES_FUN_POINTER */
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(MP_ERROR_TEXT("axis is out of bounds"));
}
return ax;
}
shape_strides tools_reduce_axes(ndarray_obj_t *ndarray, mp_obj_t axis) {
// TODO: replace numerical_reduce_axes with this function, wherever applicable
// This function should be used, whenever a tensor is contracted;
@ -172,38 +181,36 @@ shape_strides tools_reduce_axes(ndarray_obj_t *ndarray, mp_obj_t axis) {
}
shape_strides _shape_strides;
_shape_strides.increment = 0;
// this is the contracted dimension (won't be overwritten for axis == None)
_shape_strides.ndim = 0;
if(axis == mp_const_none) {
_shape_strides.shape = ndarray->shape;
_shape_strides.strides = ndarray->strides;
return _shape_strides;
}
size_t *shape = m_new(size_t, ULAB_MAX_DIMS + 1);
_shape_strides.shape = shape;
int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS + 1);
_shape_strides.strides = strides;
_shape_strides.increment = 0;
// this is the contracted dimension (won't be overwritten for axis == None)
_shape_strides.ndim = 0;
memcpy(_shape_strides.shape, ndarray->shape, sizeof(size_t) * ULAB_MAX_DIMS);
memcpy(_shape_strides.strides, ndarray->strides, sizeof(int32_t) * ULAB_MAX_DIMS);
if(axis == mp_const_none) {
return _shape_strides;
}
uint8_t index = ULAB_MAX_DIMS - 1; // value of index for axis == mp_const_none (won't be overwritten)
_shape_strides.axis = ULAB_MAX_DIMS - 1; // value of index for axis == mp_const_none (won't be overwritten)
if(axis != mp_const_none) { // i.e., axis is an integer
int8_t ax = mp_obj_get_int(axis);
if(ax < 0) ax += ndarray->ndim;
if((ax < 0) || (ax > ndarray->ndim - 1)) {
mp_raise_ValueError(MP_ERROR_TEXT("index out of range"));
}
index = ULAB_MAX_DIMS - ndarray->ndim + ax;
int8_t ax = tools_get_axis(axis, ndarray->ndim);
_shape_strides.axis = ULAB_MAX_DIMS - ndarray->ndim + ax;
_shape_strides.ndim = ndarray->ndim - 1;
}
// move the value stored at index to the leftmost position, and align everything else to the right
_shape_strides.shape[0] = ndarray->shape[index];
_shape_strides.strides[0] = ndarray->strides[index];
for(uint8_t i = 0; i < index; i++) {
_shape_strides.shape[0] = ndarray->shape[_shape_strides.axis];
_shape_strides.strides[0] = ndarray->strides[_shape_strides.axis];
for(uint8_t i = 0; i < _shape_strides.axis; i++) {
// entries to the right of index must be shifted by one position to the left
_shape_strides.shape[i + 1] = ndarray->shape[i];
_shape_strides.strides[i + 1] = ndarray->strides[i];
@ -213,16 +220,37 @@ shape_strides tools_reduce_axes(ndarray_obj_t *ndarray, mp_obj_t axis) {
_shape_strides.increment = 1;
}
if(_shape_strides.ndim == 0) {
_shape_strides.ndim = 1;
_shape_strides.shape[ULAB_MAX_DIMS - 1] = 1;
_shape_strides.strides[ULAB_MAX_DIMS - 1] = ndarray->itemsize;
}
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(MP_ERROR_TEXT("axis is out of bounds"));
mp_obj_t ulab_tools_restore_dims(ndarray_obj_t *ndarray, ndarray_obj_t *results, mp_obj_t keepdims, shape_strides _shape_strides) {
// restores the contracted dimension, if keepdims is True
if((ndarray->ndim == 1) && (keepdims != mp_const_true)) {
// since the original array has already been contracted and
// we don't want to keep the dimensions here, we have to return a scalar
return mp_binary_get_val_array(results->dtype, results->array, 0);
}
return ax;
if(keepdims == mp_const_true) {
results->ndim += 1;
for(int8_t i = 0; i < ULAB_MAX_DIMS; i++) {
results->shape[i] = ndarray->shape[i];
}
results->shape[_shape_strides.axis] = 1;
results->strides[ULAB_MAX_DIMS - 1] = ndarray->itemsize;
for(uint8_t i = ULAB_MAX_DIMS; i > 1; i--) {
results->strides[i - 2] = results->strides[i - 1] * results->shape[i - 1];
}
}
return MP_OBJ_FROM_PTR(results);
}
#if ULAB_MAX_DIMS > 1

2
code/ulab_tools.h
View file

@ -17,6 +17,7 @@
typedef struct _shape_strides_t {
uint8_t increment;
uint8_t axis;
uint8_t ndim;
size_t *shape;
int32_t *strides;
@ -34,6 +35,7 @@ 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 );
mp_obj_t ulab_tools_restore_dims(ndarray_obj_t * , ndarray_obj_t * , mp_obj_t , shape_strides );
ndarray_obj_t *tools_object_is_square(mp_obj_t );
uint8_t ulab_binary_get_size(uint8_t );

5
docs/manual/source/conf.py
View file

@ -23,11 +23,12 @@ from sphinx import addnodes
# -- Project information -----------------------------------------------------
project = 'The ulab book'
copyright = '2019-2024, Zoltán Vörös and contributors'
copyright = '2019-2025, Zoltán Vörös and contributors'
author = 'Zoltán Vörös'
# The full version, including alpha/beta/rc tags
release = '6.6.0'
release = '6.9.0'
# -- General configuration ---------------------------------------------------

1
docs/manual/source/index.rst
View file

@ -23,6 +23,7 @@ Welcome to the ulab book!
numpy-fft
numpy-linalg
numpy-random
scipy-integrate
scipy-linalg
scipy-optimize
scipy-signal

220
docs/manual/source/scipy-integrate.rst Normal file
View file

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

12
docs/manual/source/ulab-ndarray.rst
View file

@ -1814,12 +1814,12 @@ array.
Binary operators
================
``ulab`` implements the ``+``, ``-``, ``*``, ``/``, ``**``, ``<``,
``>``, ``<=``, ``>=``, ``==``, ``!=``, ``+=``, ``-=``, ``*=``, ``/=``,
``**=`` binary operators, as well as the ``AND``, ``OR``, ``XOR``
bit-wise operators that work element-wise. Note that the bit-wise
operators will raise an exception, if either of the operands is of
``float`` or ``complex`` type.
``ulab`` implements the ``+``, ``-``, ``*``, ``/``, ``**``, ``%``,
``<``, ``>``, ``<=``, ``>=``, ``==``, ``!=``, ``+=``, ``-=``, ``*=``,
``/=``, ``**=``, ``%=`` binary operators, as well as the ``AND``,
``OR``, ``XOR`` bit-wise operators that work element-wise. Note that the
bit-wise operators will raise an exception, if either of the operands is
of ``float`` or ``complex`` type.
Broadcasting is available, meaning that the two operands do not even
have to have the same shape. If the lengths along the respective axes

54
docs/ulab-change-log.md
View file

@ -1,3 +1,57 @@
Fri, 06 Jun 2025
version 6.8.0
expose ndim property
Fri, 06 Jun 2025
version 6.7.7
fix ndarray type inference for micropython objects
Thu, 29 May 2025
version 6.7.6
loadtxt can deal with multi-line comments
Thu, 29 May 2025
version 6.7.5
fix typo and shape in radnom module
Sun, 16 Mar 2025
version 6.7.4
re-name integration constants to avoid name clash with EPS ports
Sun, 26 Jan 2025
version 6.7.3
fix keepdims for min, max, argmin, argmax
Sun, 19 Jan 2025
version 6.7.2
fix keepdims for std, remove redundant macros from numerical.h, update documentation
Mon, 30 Dec 2024
version 6.7.1
add keepdims keyword argument to numerical functions
Sun, 15 Dec 2024
version 6.7.0
add scipy.integrate module
Sun, 24 Nov 2024
version 6.6.1

19
docs/ulab-convert.ipynb
View file

@ -14,7 +14,7 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2022-02-09T06:27:15.118699Z",
@ -57,11 +57,11 @@
"# -- Project information -----------------------------------------------------\n",
"\n",
"project = 'The ulab book'\n",
"copyright = '2019-2024, Zoltán Vörös and contributors'\n",
"copyright = '2019-2025, Zoltán Vörös and contributors'\n",
"author = 'Zoltán Vörös'\n",
"\n",
"# The full version, including alpha/beta/rc tags\n",
"release = '6.6.0'\n",
"release = '6.9.0'\n",
"\n",
"\n",
"# -- General configuration ---------------------------------------------------\n",
@ -191,6 +191,7 @@
" numpy-fft\n",
" numpy-linalg\n",
" numpy-random\n",
" scipy-integrate\n",
" scipy-linalg\n",
" scipy-optimize\n",
" scipy-signal\n",
@ -294,6 +295,10 @@
"/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n",
" _, nbc = validator.normalize(nbc)\n",
"/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n",
" _, nbc = validator.normalize(nbc)\n",
"/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n",
" _, nbc = validator.normalize(nbc)\n",
"/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n",
" _, nbc = validator.normalize(nbc)\n"
]
}
@ -306,6 +311,7 @@
" 'numpy-fft',\n",
" 'numpy-linalg',\n",
" 'numpy-random',\n",
" 'scipy-integrate',\n",
" 'scipy-linalg',\n",
" 'scipy-optimize',\n",
" 'scipy-signal',\n",
@ -470,7 +476,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.8.5 ('base')",
"display_name": "base",
"language": "python",
"name": "python3"
},
@ -532,11 +538,6 @@
"_Feature"
],
"window_display": false
},
"vscode": {
"interpreter": {
"hash": "9e4ec6f642f986afcc9e252c165e44859a62defc5c697cae6f82c2943465ec10"
}
}
},
"nbformat": 4,

2
docs/ulab-ndarray.ipynb
View file

@ -2599,7 +2599,7 @@
"source": [
"# Binary operators\n",
"\n",
"`ulab` implements the `+`, `-`, `*`, `/`, `**`, `<`, `>`, `<=`, `>=`, `==`, `!=`, `+=`, `-=`, `*=`, `/=`, `**=` binary operators, as well as the `AND`, `OR`, `XOR` bit-wise operators that work element-wise. Note that the bit-wise operators will raise an exception, if either of the operands is of `float` or `complex` type.\n",
"`ulab` implements the `+`, `-`, `*`, `/`, `**`, `%`, `<`, `>`, `<=`, `>=`, `==`, `!=`, `+=`, `-=`, `*=`, `/=`, `**=`, `%=` binary operators, as well as the `AND`, `OR`, `XOR` bit-wise operators that work element-wise. Note that the bit-wise operators will raise an exception, if either of the operands is of `float` or `complex` type.\n",
"\n",
"Broadcasting is available, meaning that the two operands do not even have to have the same shape. If the lengths along the respective axes are equal, or one of them is 1, or the axis is missing, the element-wise operation can still be carried out. \n",
"A thorough explanation of broadcasting can be found under https://numpy.org/doc/stable/user/basics.broadcasting.html. \n",

353
docs/ulab-numerical.ipynb
View file

@ -14,6 +14,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"%pylab is deprecated, use %matplotlib inline and import the required libraries.\n",
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
@ -31,7 +32,7 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-07T19:16:29.118001Z",
@ -77,7 +78,7 @@
" if args.unix: # tests the code on the unix port. Note that this works on unix only\n",
" with open('/dev/shm/micropython.py', 'w') as fout:\n",
" fout.write(cell)\n",
" proc = subprocess.Popen([\"../micropython/ports/unix/micropython-2\", \"/dev/shm/micropython.py\"], \n",
" proc = subprocess.Popen([\"../micropython/ports/unix/build-2/micropython-2\", \"/dev/shm/micropython.py\"], \n",
" stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
" print(proc.stdout.read().decode(\"utf-8\"))\n",
" print(proc.stderr.read().decode(\"utf-8\"))\n",
@ -182,7 +183,7 @@
"%%micropython -pyboard 1\n",
"\n",
"import utime\n",
"import ulab as np\n",
"from ulab import numpy as np\n",
"\n",
"def timeit(n=1000):\n",
" def wrapper(f, *args, **kwargs):\n",
@ -244,145 +245,14 @@
"\n",
"**WARNING:** Difference to `numpy`: the `out` keyword argument is not implemented.\n",
"\n",
"These functions follow the same pattern, and work with generic iterables, and `ndarray`s. `min`, and `max` return the minimum or maximum of a sequence. If the input array is two-dimensional, the `axis` keyword argument can be supplied, in which case the minimum/maximum along the given axis will be returned. If `axis=None` (this is also the default value), the minimum/maximum of the flattened array will be determined.\n",
"These functions follow the same pattern, and work with generic iterables, and `ndarray`s. `min`, and `max` return the minimum or maximum of a sequence. If the input array is two-dimensional, the `axis` keyword argument can be supplied, in which case the minimum/maximum along the given axis will be returned. If `axis=None` (this is also the default value), the minimum/maximum of the flattened array will be determined. The functions also accept the `keepdims=True` or `keepdims=False` keyword argument. The latter case is the default, while the former keeps the dimensions (the number of axes) of the supplied array. \n",
"\n",
"`argmin/argmax` return the position (index) of the minimum/maximum in the sequence."
]
},
{
"cell_type": "code",
"execution_count": 108,
"metadata": {
"ExecuteTime": {
"end_time": "2020-10-17T21:26:22.507996Z",
"start_time": "2020-10-17T21:26:22.492543Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([1.0, 2.0, 3.0], dtype=float)\n",
"array([], dtype=float)\n",
"[] 0\n",
"array([1.0, 2.0, 3.0], dtype=float)\n",
"array([], dtype=float)\n",
"\n",
"\n"
]
}
],
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"\n",
"a = np.array([1, 2, 3])\n",
"print(a)\n",
"print(a[-1:-1:-3])\n",
"try:\n",
" sa = list(a[-1:-1:-3])\n",
" la = len(sa)\n",
"except IndexError as e:\n",
" sa = str(e)\n",
" la = -1\n",
" \n",
"print(sa, la)\n",
"\n",
"a[-1:-1:-3] = np.ones(0)\n",
"print(a)\n",
"\n",
"b = np.ones(0) + 1\n",
"print(b)\n",
"# print('b', b.shape())"
]
},
{
"cell_type": "code",
"execution_count": 122,
"metadata": {
"ExecuteTime": {
"end_time": "2020-10-17T21:54:49.123748Z",
"start_time": "2020-10-17T21:54:49.093819Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0, 1, -3array([], dtype=float)\n",
"\n",
"\n"
]
}
],
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"a = np.array([1, 2, 3])\n",
"print(a[0:1:-3])"
]
},
{
"cell_type": "code",
"execution_count": 127,
"metadata": {
"ExecuteTime": {
"end_time": "2020-10-17T21:57:01.482277Z",
"start_time": "2020-10-17T21:57:01.477362Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"[0]"
]
},
"execution_count": 127,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"l = list(range(13))\n",
"\n",
"l[0:10:113]"
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {
"ExecuteTime": {
"end_time": "2020-10-17T20:59:58.285134Z",
"start_time": "2020-10-17T20:59:58.263605Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(0,)"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = np.array([1, 2, 3])\n",
"np.ones(0, dtype=uint8) / np.zeros(0, dtype=uint16)\n",
"np.ones(0).shape"
]
},
{
"cell_type": "code",
"execution_count": 375,
"execution_count": 10,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-18T13:08:28.113525Z",
@ -394,16 +264,16 @@
"name": "stdout",
"output_type": "stream",
"text": [
"a: array([1.0, 2.0, 0.0, 1.0, 10.0], dtype=float)\n",
"a: array([1.0, 2.0, 0.0, 1.0, 10.0], dtype=float64)\n",
"min of a: 0.0\n",
"argmin of a: 2\n",
"\n",
"b:\n",
" array([[1.0, 2.0, 0.0],\n",
"\t [1.0, 10.0, -1.0]], dtype=float)\n",
" [1.0, 10.0, -1.0]], dtype=float64)\n",
"min of b (flattened): -1.0\n",
"min of b (axis=0): array([1.0, 2.0, -1.0], dtype=float)\n",
"min of b (axis=1): array([0.0, -1.0], dtype=float)\n",
"min of b (axis=0): array([1.0, 2.0, -1.0], dtype=float64)\n",
"min of b (axis=1): array([0.0, -1.0], dtype=float64)\n",
"\n",
"\n"
]
@ -412,19 +282,55 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array([1, 2, 0, 1, 10])\n",
"print('a:', a)\n",
"print('min of a:', numerical.min(a))\n",
"print('argmin of a:', numerical.argmin(a))\n",
"print('min of a:', np.min(a))\n",
"print('argmin of a:', np.argmin(a))\n",
"\n",
"b = np.array([[1, 2, 0], [1, 10, -1]])\n",
"print('\\nb:\\n', b)\n",
"print('min of b (flattened):', numerical.min(b))\n",
"print('min of b (axis=0):', numerical.min(b, axis=0))\n",
"print('min of b (axis=1):', numerical.min(b, axis=1))"
"print('min of b (flattened):', np.min(b))\n",
"print('min of b (axis=0):', np.min(b, axis=0))\n",
"print('min of b (axis=1):', np.min(b, axis=1))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a: array([[0.0, 1.0, 2.0, 3.0],\n",
" [4.0, 5.0, 6.0, 7.0],\n",
" [8.0, 9.0, 10.0, 11.0]], dtype=float64)\n",
"\n",
"min of a (axis=1):\n",
" array([[0.0],\n",
" [4.0],\n",
" [8.0]], dtype=float64)\n",
"\n",
"min of a (axis=0):\n",
" array([[0.0, 1.0, 2.0, 3.0]], dtype=float64)\n",
"\n",
"\n"
]
}
],
"source": [
"%%micropython -unix 1\n",
"\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array(range(12)).reshape((3, 4))\n",
"\n",
"print('a:', a)\n",
"print('\\nmin of a (axis=1):\\n', np.min(a, axis=1, keepdims=True))\n",
"print('\\nmin of a (axis=0):\\n', np.min(a, axis=0, keepdims=True))"
]
},
{
@ -439,12 +345,14 @@
"\n",
"`numpy`: https://docs.scipy.org/doc/numpy/reference/generated/numpy.mean.html\n",
"\n",
"These three functions follow the same pattern: if the axis keyword is not specified, it assumes the default value of `None`, and returns the result of the computation for the flattened array. Otherwise, the calculation is along the given axis."
"These three functions follow the same pattern: if the `axis` keyword is not specified, they assume the default value of `None`, and return the result of the computation for the flattened array. Otherwise, the calculation is along the given axis. \n",
"\n",
"If the `axis` keyword argument is a number (this can also be negative to signify counting from the rightmost axis) the functions contract the arrays, i.e., the results will have one axis fewer than the input array. The only exception to this rule is when the `keepdims` keyword argument is supplied with a value `True`, in which case, the results will have the same number of axis as the input, but the axis specified in `axis` will have a length of 1. This is useful in cases, when the output is to be broadcast with the input in subsequent computations."
]
},
{
"cell_type": "code",
"execution_count": 527,
"execution_count": 12,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-20T06:51:58.845076Z",
@ -458,29 +366,73 @@
"text": [
"a: \n",
" array([[1.0, 2.0, 3.0],\n",
"\t [4.0, 5.0, 6.0],\n",
"\t [7.0, 8.0, 9.0]], dtype=float)\n",
" [4.0, 5.0, 6.0],\n",
" [7.0, 8.0, 9.0]], dtype=float64)\n",
"sum, flat array: 45.0\n",
"mean, horizontal: array([2.0, 5.0, 8.0], dtype=float)\n",
"std, vertical: array([2.44949, 2.44949, 2.44949], dtype=float)\n",
"mean, horizontal: array([2.0, 5.0, 8.0], dtype=float64)\n",
"std, vertical: array([2.449489742783178, 2.449489742783178, 2.449489742783178], dtype=float64)\n",
"\n",
"\n"
]
}
],
"source": [
"%%micropython -pyboard 1\n",
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])\n",
"print('a: \\n', a)\n",
"\n",
"print('sum, flat array: ', numerical.sum(a))\n",
"print('sum, flat array: ', np.sum(a))\n",
"\n",
"print('mean, horizontal: ', numerical.mean(a, axis=1))\n",
"print('mean, horizontal: ', np.mean(a, axis=1))\n",
"\n",
"print('std, vertical: ', numerical.std(a, axis=0))"
"print('std, vertical: ', np.std(a, axis=0))"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a: \n",
" array([[1.0, 2.0, 3.0],\n",
" [4.0, 5.0, 6.0],\n",
" [7.0, 8.0, 9.0]], dtype=float64)\n",
"\n",
"std, along 0th axis:\n",
" array([2.449489742783178, 2.449489742783178, 2.449489742783178], dtype=float64)\n",
"\n",
"a: \n",
" array([[1.0, 2.0, 3.0],\n",
" [4.0, 5.0, 6.0],\n",
" [7.0, 8.0, 9.0]], dtype=float64)\n",
"\n",
"std, along 1st axis, keeping dimensions:\n",
" array([[0.8164965809277261],\n",
" [0.8164965809277261],\n",
" [0.8164965809277261]], dtype=float64)\n",
"\n",
"\n"
]
}
],
"source": [
"%%micropython -unix 1\n",
"\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])\n",
"print('a: \\n', a)\n",
"print('\\nstd, along 0th axis:\\n', np.std(a, axis=0))\n",
"\n",
"print('\\na: \\n', a)\n",
"print('\\nstd, along 1st axis, keeping dimensions:\\n', np.std(a, axis=1, keepdims=True))"
]
},
{
@ -519,13 +471,12 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array([1, 2, 3, 4, 5, 6, 7, 8])\n",
"print(\"a:\\t\\t\\t\", a)\n",
"\n",
"numerical.roll(a, 2)\n",
"np.roll(a, 2)\n",
"print(\"a rolled to the left:\\t\", a)\n",
"\n",
"# this should be the original vector\n",
@ -581,19 +532,18 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])\n",
"print(\"a:\\n\", a)\n",
"\n",
"numerical.roll(a, 2)\n",
"np.roll(a, 2)\n",
"print(\"\\na rolled to the left:\\n\", a)\n",
"\n",
"numerical.roll(a, -1, axis=1)\n",
"np.roll(a, -1, axis=1)\n",
"print(\"\\na rolled up:\\n\", a)\n",
"\n",
"numerical.roll(a, 1, axis=None)\n",
"np.roll(a, 1, axis=None)\n",
"print(\"\\na rolled with None:\\n\", a)"
]
},
@ -649,9 +599,7 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import vector\n",
"from ulab import numpy as np\n",
"\n",
"def dummy_adc():\n",
" # dummy adc function, so that the results are reproducible\n",
@ -659,8 +607,8 @@
" \n",
"n = 10\n",
"# These are the normalised weights; the last entry is the most dominant\n",
"weight = vector.exp([1, 2, 3, 4, 5])\n",
"weight = weight/numerical.sum(weight)\n",
"weight = np.exp([1, 2, 3, 4, 5])\n",
"weight = weight/np.sum(weight)\n",
"\n",
"print(weight)\n",
"# initial array of samples\n",
@ -669,10 +617,10 @@
"for i in range(n):\n",
" # a new datum is inserted on the right hand side. This simply overwrites whatever was in the last slot\n",
" samples[-1] = dummy_adc()\n",
" print(numerical.mean(samples[-5:]*weight))\n",
" print(np.mean(samples[-5:]*weight))\n",
" print(samples[-5:])\n",
" # the data are shifted by one position to the left\n",
" numerical.roll(samples, 1)"
" numerical.np(samples, 1)"
]
},
{
@ -725,17 +673,16 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array([1, 2, 3, 4, 5])\n",
"print(\"a: \\t\", a)\n",
"print(\"a flipped:\\t\", np.flip(a))\n",
"\n",
"a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=np.uint8)\n",
"print(\"\\na flipped horizontally\\n\", numerical.flip(a, axis=1))\n",
"print(\"\\na flipped vertically\\n\", numerical.flip(a, axis=0))\n",
"print(\"\\na flipped horizontally+vertically\\n\", numerical.flip(a))"
"print(\"\\na flipped horizontally\\n\", np.flip(a, axis=1))\n",
"print(\"\\na flipped vertically\\n\", np.flip(a, axis=0))\n",
"print(\"\\na flipped horizontally+vertically\\n\", np.flip(a))"
]
},
{
@ -800,19 +747,18 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array(range(9), dtype=np.uint8)\n",
"print('a:\\n', a)\n",
"\n",
"print('\\nfirst derivative:\\n', numerical.diff(a, n=1))\n",
"print('\\nsecond derivative:\\n', numerical.diff(a, n=2))\n",
"print('\\nfirst derivative:\\n', np.diff(a, n=1))\n",
"print('\\nsecond derivative:\\n', np.diff(a, n=2))\n",
"\n",
"c = np.array([[1, 2, 3, 4], [4, 3, 2, 1], [1, 4, 9, 16], [0, 0, 0, 0]])\n",
"print('\\nc:\\n', c)\n",
"print('\\nfirst derivative, first axis:\\n', numerical.diff(c, axis=0))\n",
"print('\\nfirst derivative, second axis:\\n', numerical.diff(c, axis=1))"
"print('\\nfirst derivative, first axis:\\n', np.diff(c, axis=0))\n",
"print('\\nfirst derivative, second axis:\\n', np.diff(c, axis=1))"
]
},
{
@ -858,7 +804,7 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array(range(12), dtype=np.int8).reshape((3, 4))\n",
"print('a:\\n', a)\n",
@ -925,18 +871,17 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array([[1, 12, 3, 0], [5, 3, 4, 1], [9, 11, 1, 8], [7, 10, 0, 1]], dtype=np.float)\n",
"print('\\na:\\n', a)\n",
"b = numerical.sort(a, axis=0)\n",
"b = np.sort(a, axis=0)\n",
"print('\\na sorted along vertical axis:\\n', b)\n",
"\n",
"c = numerical.sort(a, axis=1)\n",
"c = np.sort(a, axis=1)\n",
"print('\\na sorted along horizontal axis:\\n', c)\n",
"\n",
"c = numerical.sort(a, axis=None)\n",
"c = np.sort(a, axis=None)\n",
"print('\\nflattened a sorted:\\n', c)"
]
},
@ -955,15 +900,13 @@
"source": [
"%%micropython -pyboard 1\n",
"\n",
"import ulab as np\n",
"from ulab import vector\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"@timeit\n",
"def sort_time(array):\n",
" return numerical.sort(array)\n",
" return np.sort(array)\n",
"\n",
"b = vector.sin(np.linspace(0, 6.28, num=1000))\n",
"b = np.sin(np.linspace(0, 6.28, num=1000))\n",
"print('b: ', b)\n",
"sort_time(b)\n",
"print('\\nb sorted:\\n', b)"
@ -1025,18 +968,17 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array([[1, 12, 3, 0], [5, 3, 4, 1], [9, 11, 1, 8], [7, 10, 0, 1]], dtype=np.float)\n",
"print('\\na:\\n', a)\n",
"b = numerical.argsort(a, axis=0)\n",
"b = np.argsort(a, axis=0)\n",
"print('\\na sorted along vertical axis:\\n', b)\n",
"\n",
"c = numerical.argsort(a, axis=1)\n",
"c = np.argsort(a, axis=1)\n",
"print('\\na sorted along horizontal axis:\\n', c)\n",
"\n",
"c = numerical.argsort(a, axis=None)\n",
"c = np.argsort(a, axis=None)\n",
"print('\\nflattened a sorted:\\n', c)"
]
},
@ -1078,12 +1020,11 @@
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"from ulab import numerical\n",
"from ulab import numpy as np\n",
"\n",
"a = np.array([0, 5, 1, 3, 2, 4], dtype=np.uint8)\n",
"print('\\na:\\n', a)\n",
"b = numerical.argsort(a, axis=1)\n",
"b = np.argsort(a, axis=1)\n",
"print('\\nsorting indices:\\n', b)\n",
"print('\\nthe original array:\\n', a)"
]
@ -1091,7 +1032,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "base",
"language": "python",
"name": "python3"
},
@ -1105,7 +1046,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
"version": "3.11.7"
},
"toc": {
"base_numbering": 1,

26
tests/2d/numpy/modulo.py Normal file
View file

@ -0,0 +1,26 @@
try:
from ulab import numpy as np
except:
import numpy as np
dtypes = (np.uint8, np.int8, np.uint16, np.int16, np.float)
for dtype1 in dtypes:
x1 = np.array(range(6), dtype=dtype1).reshape((2, 3))
for dtype2 in dtypes:
x2 = np.array(range(1, 4), dtype=dtype2)
print(x1 % x2)
print()
print('=' * 30)
print('inplace modulo')
print('=' * 30)
print()
for dtype1 in dtypes:
x1 = np.array(range(6), dtype=dtype1).reshape((2, 3))
for dtype2 in dtypes:
x2 = np.array(range(1, 4), dtype=dtype2)
x1 %= x2
print(x1)

105
tests/2d/numpy/modulo.py.exp Normal file
View file

@ -0,0 +1,105 @@
array([[0, 1, 2],
[0, 0, 2]], dtype=uint8)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0, 1, 2],
[0, 0, 2]], dtype=uint16)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0, 1, 2],
[0, 0, 2]], dtype=int8)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0, 0, 1],
[0, 2, 0]], dtype=uint8)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0, 1, 2],
[0, 0, 2]], dtype=uint16)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
==============================
inplace modulo
==============================
array([[0, 1, 2],
[0, 0, 2]], dtype=uint8)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0, 0, 1],
[0, 2, 0]], dtype=uint8)
array([[0, 0, 1],
[0, 0, 0]], dtype=int16)
array([[0.0, 0.0, 1.0],
[0.0, 0.0, 0.0]], dtype=float64)
array([[0.0, 0.0, 1.0],
[0.0, 0.0, 0.0]], dtype=float64)
array([[0.0, 0.0, 1.0],
[0.0, 0.0, 0.0]], dtype=float64)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0, 1, 2],
[0, 0, 2]], dtype=int16)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)
array([[0.0, 1.0, 2.0],
[0.0, 0.0, 2.0]], dtype=float64)

10
tests/2d/numpy/random.py Normal file
View file

@ -0,0 +1,10 @@
try:
from ulab import numpy as np
except ImportError:
import numpy as np
rng = np.random.Generator(1234)
for generator in (rng.normal, rng.random, rng.uniform):
random_array = generator(size=(1, 2))
print("array shape:", random_array.shape)

3
tests/2d/numpy/random.py.exp Normal file
View file

@ -0,0 +1,3 @@
array shape: (1, 2)
array shape: (1, 2)
array shape: (1, 2)

23
tests/2d/numpy/sum.py Normal file
View file

@ -0,0 +1,23 @@
try:
from ulab import numpy as np
except ImportError:
import numpy as np
for dtype in (np.uint8, np.int8, np.uint16, np.int8, np.float):
a = np.array(range(12), dtype=dtype)
b = a.reshape((3, 4))
print(a)
print(b)
print()
print(np.sum(a))
print(np.sum(a, axis=0))
print(np.sum(a, axis=0, keepdims=True))
print()
print(np.sum(b))
print(np.sum(b, axis=0))
print(np.sum(b, axis=1))
print(np.sum(b, axis=0, keepdims=True))
print(np.sum(b, axis=1, keepdims=True))

80
tests/2d/numpy/sum.py.exp Normal file
View file

@ -0,0 +1,80 @@
array([0, 1, 2, ..., 9, 10, 11], dtype=uint8)
array([[0, 1, 2, 3],
[4, 5, 6, 7],
[8, 9, 10, 11]], dtype=uint8)
66
66
array([66], dtype=uint8)
66
array([12, 15, 18, 21], dtype=uint8)
array([6, 22, 38], dtype=uint8)
array([[12, 15, 18, 21]], dtype=uint8)
array([[6],
[22],
[38]], dtype=uint8)
array([0, 1, 2, ..., 9, 10, 11], dtype=int8)
array([[0, 1, 2, 3],
[4, 5, 6, 7],
[8, 9, 10, 11]], dtype=int8)
66
66
array([66], dtype=int8)
66
array([12, 15, 18, 21], dtype=int8)
array([6, 22, 38], dtype=int8)
array([[12, 15, 18, 21]], dtype=int8)
array([[6],
[22],
[38]], dtype=int8)
array([0, 1, 2, ..., 9, 10, 11], dtype=uint16)
array([[0, 1, 2, 3],
[4, 5, 6, 7],
[8, 9, 10, 11]], dtype=uint16)
66
66
array([66], dtype=uint16)
66
array([12, 15, 18, 21], dtype=uint16)
array([6, 22, 38], dtype=uint16)
array([[12, 15, 18, 21]], dtype=uint16)
array([[6],
[22],
[38]], dtype=uint16)
array([0, 1, 2, ..., 9, 10, 11], dtype=int8)
array([[0, 1, 2, 3],
[4, 5, 6, 7],
[8, 9, 10, 11]], dtype=int8)
66
66
array([66], dtype=int8)
66
array([12, 15, 18, 21], dtype=int8)
array([6, 22, 38], dtype=int8)
array([[12, 15, 18, 21]], dtype=int8)
array([[6],
[22],
[38]], dtype=int8)
array([0.0, 1.0, 2.0, ..., 9.0, 10.0, 11.0], dtype=float64)
array([[0.0, 1.0, 2.0, 3.0],
[4.0, 5.0, 6.0, 7.0],
[8.0, 9.0, 10.0, 11.0]], dtype=float64)
66.0
66.0
array([66.0], dtype=float64)
66.0
array([12.0, 15.0, 18.0, 21.0], dtype=float64)
array([6.0, 22.0, 38.0], dtype=float64)
array([[12.0, 15.0, 18.0, 21.0]], dtype=float64)
array([[6.0],
[22.0],
[38.0]], dtype=float64)