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

fixed glitch in std/mean functions

Browse source
This commit is contained in:
Zoltán Vörös 2020-01-07 21:53:23 +01:00
parent 992e48b84e
commit ac1111c251
5 changed files with 505 additions and 424 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

382
code/numerical.c
View file

@ -7,7 +7,7 @@
*
* Copyright (c) 2019-2020 Zoltán Vörös
*/
#include <math.h>
#include <stdlib.h>
#include <string.h>
@ -77,7 +77,33 @@ mp_obj_t numerical_linspace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *k
}
}
mp_obj_t numerical_sum_mean_std_array(mp_obj_t oin, uint8_t optype) {
void axis_sorter(ndarray_obj_t *ndarray, mp_obj_t axis, size_t *m, size_t *n, size_t *N,
size_t *increment, size_t *len, size_t *start_inc) {
if(axis == mp_const_none) { // flatten the array
*m = 1;
*n = 1;
*len = ndarray->array->len;
*N = 1;
*increment = 1;
*start_inc = ndarray->array->len;
} else if((mp_obj_get_int(axis) == 1)) { // along the horizontal axis
*m = ndarray->m;
*n = 1;
*len = ndarray->n;
*N = ndarray->m;
*increment = 1;
*start_inc = ndarray->n;
} else { // along vertical axis
*m = 1;
*n = ndarray->n;
*len = ndarray->m;
*N = ndarray->n;
*increment = ndarray->n;
*start_inc = 1;
}
}
mp_obj_t numerical_sum_mean_std_iterable(mp_obj_t oin, uint8_t optype, size_t ddof) {
mp_float_t value, sum = 0.0, sq_sum = 0.0;
mp_obj_iter_buf_t iter_buf;
mp_obj_t item, iterable = mp_getiter(oin, &iter_buf);
@ -85,194 +111,153 @@ mp_obj_t numerical_sum_mean_std_array(mp_obj_t oin, uint8_t optype) {
while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
value = mp_obj_get_float(item);
sum += value;
if(optype == NUMERICAL_STD) {
sq_sum += value*value;
}
}
if(optype == NUMERICAL_SUM) {
return mp_obj_new_float(sum);
} else if(optype == NUMERICAL_MEAN) {
return mp_obj_new_float(sum/len);
} else {
} else { // this should be the case of the standard deviation
// TODO: note that we could get away with a single pass, if we used the Weldorf algorithm
// That should save a fair amount of time, because we would have to extract the values only once
iterable = mp_getiter(oin, &iter_buf);
sum /= len; // this is now the mean!
return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/len-sum*sum));
}
}
STATIC mp_float_t numerical_sum_mean_std_single_line(void *data, size_t start, size_t stop,
size_t stride, uint8_t typecode, uint8_t optype) {
mp_float_t sum = 0.0, sq_sum = 0.0, value;
size_t len = 0;
for(size_t i=start; i < stop; i+=stride, len++) {
value = ndarray_get_float_value(data, typecode, i);
sum += value;
if(optype == NUMERICAL_STD) {
sq_sum += value*value;
}
}
if(len == 0) {
mp_raise_ValueError("data length is 0!");
}
if(optype == NUMERICAL_SUM) {
return sum;
} else if(optype == NUMERICAL_MEAN) {
return sum/len;
} else {
sum /= len; // this is now the mean!
return MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/len-sum*sum);
}
}
STATIC mp_obj_t numerical_sum_mean_std_matrix(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {
ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
if((axis == mp_const_none) || (in->m == 1) || (in->n == 1)) {
// return the value for the flattened array
return mp_obj_new_float(numerical_sum_mean_std_single_line(in->array->items, 0,
in->array->len, 1, in->array->typecode, optype));
} else {
uint8_t _axis = mp_obj_get_int(axis);
size_t m = (_axis == 0) ? 1 : in->m;
size_t n = (_axis == 0) ? in->n : 1;
size_t len = in->array->len;
mp_float_t sms;
// TODO: pass in->array->typcode to create_new_ndarray
ndarray_obj_t *out = create_new_ndarray(m, n, NDARRAY_FLOAT);
// TODO: these two cases could probably be combined in a more elegant fashion...
if(_axis == 0) { // vertical
for(size_t i=0; i < n; i++) {
sms = numerical_sum_mean_std_single_line(in->array->items, i, len,
n, in->array->typecode, optype);
((float_t *)out->array->items)[i] = sms;
}
} else { // horizontal
for(size_t i=0; i < m; i++) {
sms = numerical_sum_mean_std_single_line(in->array->items, i*in->n,
(i+1)*in->n, 1, in->array->typecode, optype);
((float_t *)out->array->items)[i] = sms;
}
}
return MP_OBJ_FROM_PTR(out);
}
}
size_t numerical_argmin_argmax_array(ndarray_obj_t *in, size_t start,
size_t stop, size_t stride, uint8_t op) {
size_t best_idx = start;
if(in->array->typecode == NDARRAY_UINT8) {
ARG_MIN_LOOP(in, uint8_t, start, stop, stride, op);
} else if(in->array->typecode == NDARRAY_INT8) {
ARG_MIN_LOOP(in, int8_t, start, stop, stride, op);
} else if(in->array->typecode == NDARRAY_UINT16) {
ARG_MIN_LOOP(in, uint16_t, start, stop, stride, op);
} else if(in->array->typecode == NDARRAY_INT16) {
ARG_MIN_LOOP(in, uint16_t, start, stop, stride, op);
} else if(in->array->typecode == NDARRAY_FLOAT) {
ARG_MIN_LOOP(in, mp_float_t, start, stop, stride, op);
}
return best_idx;
}
void copy_value_into_ndarray(ndarray_obj_t *target, ndarray_obj_t *source, size_t target_idx, size_t source_idx) {
// since we are simply copying, it doesn't matter, whether the arrays are signed or unsigned,
// we can cast them in any way we like
// This could also be done with byte copies. I don't know, whether that would have any benefits
if((target->array->typecode == NDARRAY_UINT8) || (target->array->typecode == NDARRAY_INT8)) {
((uint8_t *)target->array->items)[target_idx] = ((uint8_t *)source->array->items)[source_idx];
} else if((target->array->typecode == NDARRAY_UINT16) || (target->array->typecode == NDARRAY_INT16)) {
((uint16_t *)target->array->items)[target_idx] = ((uint16_t *)source->array->items)[source_idx];
} else {
((float *)target->array->items)[target_idx] = ((float *)source->array->items)[source_idx];
}
}
STATIC mp_obj_t numerical_argmin_argmax(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {
if(mp_obj_is_type(oin, &mp_type_tuple) || mp_obj_is_type(oin, &mp_type_list) ||
mp_obj_is_type(oin, &mp_type_range)) {
// This case will work for single iterables only
size_t idx = 0, best_idx = 0;
mp_obj_iter_buf_t iter_buf;
mp_obj_t iterable = mp_getiter(oin, &iter_buf);
mp_obj_t best_obj = MP_OBJ_NULL;
mp_obj_t item;
mp_uint_t op = MP_BINARY_OP_LESS;
if((optype == NUMERICAL_ARGMAX) || (optype == NUMERICAL_MAX)) op = MP_BINARY_OP_MORE;
while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
if ((best_obj == MP_OBJ_NULL) || (mp_binary_op(op, item, best_obj) == mp_const_true)) {
best_obj = item;
best_idx = idx;
}
idx++;
value = mp_obj_get_float(item) - sum;
sq_sum += value * value;
}
if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {
return MP_OBJ_NEW_SMALL_INT(best_idx);
} else {
return best_obj;
}
} else if(mp_obj_is_type(oin, &ulab_ndarray_type)) {
ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
size_t best_idx;
if((axis == mp_const_none) || (in->m == 1) || (in->n == 1)) {
// return the value for the flattened array
best_idx = numerical_argmin_argmax_array(in, 0, in->array->len, 1, optype);
if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {
return MP_OBJ_NEW_SMALL_INT(best_idx);
} else {
if(in->array->typecode == NDARRAY_FLOAT) {
return mp_obj_new_float(ndarray_get_float_value(in->array->items, in->array->typecode, best_idx));
} else {
return mp_binary_get_val_array(in->array->typecode, in->array->items, best_idx);
}
}
} else { // we have to work with a full matrix here
uint8_t _axis = mp_obj_get_int(axis);
size_t m = (_axis == 0) ? 1 : in->m;
size_t n = (_axis == 0) ? in->n : 1;
size_t len = in->array->len;
ndarray_obj_t *ndarray = NULL;
if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {
ndarray = create_new_ndarray(m, n, in->array->typecode);
} else { // argmin/argmax
// TODO: one might get away with uint8_t, if both m, and n < 255
ndarray = create_new_ndarray(m, n, NDARRAY_UINT16);
}
// TODO: these two cases could probably be combined in a more elegant fashion...
if(_axis == 0) { // vertical
for(size_t i=0; i < n; i++) {
best_idx = numerical_argmin_argmax_array(in, i, len, n, optype);
if((optype == NUMERICAL_MIN) || (optype == NUMERICAL_MAX)) {
copy_value_into_ndarray(ndarray, in, i, best_idx);
} else {
((uint16_t *)ndarray->array->items)[i] = (uint16_t)(best_idx / n);
}
}
} else { // horizontal
for(size_t i=0; i < m; i++) {
best_idx = numerical_argmin_argmax_array(in, i*in->n, (i+1)*in->n, 1, optype);
if((optype == NUMERICAL_MIN) || (optype == NUMERICAL_MAX)) {
copy_value_into_ndarray(ndarray, in, i, best_idx);
} else {
((uint16_t *)ndarray->array->items)[i] = (uint16_t)(best_idx - i*in->n);
}
}
}
return MP_OBJ_FROM_PTR(ndarray);
}
return mp_const_none;
}
mp_raise_TypeError("input type is not supported");
return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/(len-ddof)));
}
}
STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t type) {
STATIC mp_obj_t numerical_sum_mean_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {
size_t m, n, increment, start, start_inc, N, len;
axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);
ndarray_obj_t *results = create_new_ndarray(m, n, NDARRAY_FLOAT);
mp_float_t sum, sq_sum;
mp_float_t *farray = (mp_float_t *)results->array->items;
for(size_t j=0; j < N; j++) { // result index
start = j * start_inc;
sum = sq_sum = 0.0;
if(ndarray->array->typecode == NDARRAY_UINT8) {
RUN_SUM(ndarray, uint8_t, optype, len, start, increment);
} else if(ndarray->array->typecode == NDARRAY_INT8) {
RUN_SUM(ndarray, int8_t, optype, len, start, increment);
} else if(ndarray->array->typecode == NDARRAY_UINT16) {
RUN_SUM(ndarray, uint16_t, optype, len, start, increment);
} else if(ndarray->array->typecode == NDARRAY_INT16) {
RUN_SUM(ndarray, int16_t, optype, len, start, increment);
} else { // this will be mp_float_t, no need to check
RUN_SUM(ndarray, mp_float_t, optype, len, start, increment);
}
if(optype == NUMERICAL_SUM) {
farray[j] = sum;
} else { // this is the case of the mean
farray[j] = sum / len;
}
}
if(results->array->len == 1) {
return mp_obj_new_float(farray[0]);
}
return MP_OBJ_FROM_PTR(results);
}
mp_obj_t numerical_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, size_t ddof) {
size_t m, n, increment, start, start_inc, N, len;
mp_float_t sum, sum_sq;
axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);
if(ddof > len) {
mp_raise_ValueError("ddof must be smaller than length of data set");
}
ndarray_obj_t *results = create_new_ndarray(m, n, NDARRAY_FLOAT);
mp_float_t *farray = (mp_float_t *)results->array->items;
for(size_t j=0; j < N; j++) { // result index
start = j * start_inc;
sum = 0.0;
sum_sq = 0.0;
if(ndarray->array->typecode == NDARRAY_UINT8) {
RUN_STD(ndarray, uint8_t, len, start, increment);
} else if(ndarray->array->typecode == NDARRAY_INT8) {
RUN_STD(ndarray, int8_t, len, start, increment);
} else if(ndarray->array->typecode == NDARRAY_UINT16) {
RUN_STD(ndarray, uint16_t, len, start, increment);
} else if(ndarray->array->typecode == NDARRAY_INT16) {
RUN_STD(ndarray, int16_t, len, start, increment);
} else { // this will be mp_float_t, no need to check
RUN_STD(ndarray, mp_float_t, len, start, increment);
}
farray[j] = MICROPY_FLOAT_C_FUN(sqrt)(sum_sq/(len - ddof));
}
if(results->array->len == 1) {
return mp_obj_new_float(farray[0]);
}
return MP_OBJ_FROM_PTR(results);
}
mp_obj_t numerical_argmin_argmax_iterable(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {
size_t idx = 0, best_idx = 0;
mp_obj_iter_buf_t iter_buf;
mp_obj_t iterable = mp_getiter(oin, &iter_buf);
mp_obj_t best_obj = MP_OBJ_NULL;
mp_obj_t item;
mp_uint_t op = MP_BINARY_OP_LESS;
if((optype == NUMERICAL_ARGMAX) || (optype == NUMERICAL_MAX)) op = MP_BINARY_OP_MORE;
while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
if ((best_obj == MP_OBJ_NULL) || (mp_binary_op(op, item, best_obj) == mp_const_true)) {
best_obj = item;
best_idx = idx;
}
idx++;
}
if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {
return MP_OBJ_NEW_SMALL_INT(best_idx);
} else {
return best_obj;
}
}
mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {
size_t m, n, increment, start, start_inc, N, len;
axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);
ndarray_obj_t *results;
if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {
// we could save some RAM by taking NDARRAY_UINT8, if the dimensions
// are smaller than 256, but the code would become more involving
// (we would also need extra flash space)
results = create_new_ndarray(m, n, NDARRAY_UINT16);
} else {
results = create_new_ndarray(m, n, ndarray->array->typecode);
}
for(size_t j=0; j < N; j++) { // result index
start = j * start_inc;
if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {
if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {
RUN_ARGMIN(ndarray, results, uint8_t, uint8_t, len, start, increment, optype, j);
} else {
RUN_ARGMIN(ndarray, results, uint8_t, uint16_t, len, start, increment, optype, j);
}
} else if((ndarray->array->typecode == NDARRAY_UINT16) || (ndarray->array->typecode == NDARRAY_INT16)) {
RUN_ARGMIN(ndarray, results, uint16_t, uint16_t, len, start, increment, optype, j);
} else {
if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {
RUN_ARGMIN(ndarray, results, mp_float_t, mp_float_t, len, start, increment, optype, j);
} else {
RUN_ARGMIN(ndarray, results, mp_float_t, uint16_t, len, start, increment, optype, j);
}
}
}
return MP_OBJ_FROM_PTR(results);
}
STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t optype) {
static const mp_arg_t allowed_args[] = {
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} } ,
{ MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },
{ MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },
};
mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
mp_obj_t oin = args[0].u_obj;
mp_obj_t axis = args[1].u_obj;
@ -281,32 +266,31 @@ STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_m
mp_raise_ValueError("axis must be None, 0, or 1");
}
if(mp_obj_is_type(oin, &mp_type_tuple) || mp_obj_is_type(oin, &mp_type_list) ||
mp_obj_is_type(oin, &mp_type_range)) {
switch(type) {
if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) ||
MP_OBJ_IS_TYPE(oin, &mp_type_range)) {
switch(optype) {
case NUMERICAL_MIN:
case NUMERICAL_ARGMIN:
case NUMERICAL_MAX:
case NUMERICAL_ARGMAX:
return numerical_argmin_argmax(oin, axis, type);
return numerical_argmin_argmax_iterable(oin, axis, optype);
case NUMERICAL_SUM:
case NUMERICAL_MEAN:
case NUMERICAL_STD:
return numerical_sum_mean_std_array(oin, type);
return numerical_sum_mean_std_iterable(oin, optype, 0);
default: // we should never reach this point, but whatever
return mp_const_none;
}
} else if(mp_obj_is_type(oin, &ulab_ndarray_type)) {
switch(type) {
} else if(MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);
switch(optype) {
case NUMERICAL_MIN:
case NUMERICAL_MAX:
case NUMERICAL_ARGMIN:
case NUMERICAL_ARGMAX:
return numerical_argmin_argmax(oin, axis, type);
return numerical_argmin_argmax_ndarray(ndarray, axis, optype);
case NUMERICAL_SUM:
case NUMERICAL_MEAN:
case NUMERICAL_STD:
return numerical_sum_mean_std_matrix(oin, axis, type);
return numerical_sum_mean_ndarray(ndarray, axis, optype);
default:
mp_raise_NotImplementedError("operation is not implemented on ndarrays");
}
@ -341,7 +325,31 @@ mp_obj_t numerical_mean(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
}
mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
return numerical_function(n_args, pos_args, kw_args, NUMERICAL_STD);
static const mp_arg_t allowed_args[] = {
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} } ,
{ MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },
{ MP_QSTR_ddof, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} },
};
mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
mp_obj_t oin = args[0].u_obj;
mp_obj_t axis = args[1].u_obj;
size_t ddof = args[2].u_int;
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("axis must be None, 0, or 1");
}
if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) || MP_OBJ_IS_TYPE(oin, &mp_type_range)) {
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_std_ndarray(ndarray, axis, ddof);
} else {
mp_raise_TypeError("input must be tuple, list, range, or ndarray");
}
return mp_const_none;
}
mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {

45
code/numerical.h
View file

@ -5,7 +5,7 @@
*
* The MIT License (MIT)
*
* Copyright (c) 2019 Zoltán Vörös
* Copyright (c) 2019-2020 Zoltán Vörös
*/
#ifndef _NUMERICAL_
@ -33,19 +33,44 @@ mp_obj_t numerical_sort(size_t , const mp_obj_t *, mp_map_t *);
mp_obj_t numerical_sort_inplace(size_t , const mp_obj_t *, mp_map_t *);
mp_obj_t numerical_argsort(size_t , const mp_obj_t *, mp_map_t *);
// this macro could be tighter, if we moved the ifs to the argmin function, assigned <, as well as >
#define ARG_MIN_LOOP(in, type, start, stop, stride, op) do {\
type *array = (type *)(in)->array->items;\
#define RUN_ARGMIN(in, out, typein, typeout, len, start, increment, op, pos) do {\
typein *array = (typein *)(in)->array->items;\
typeout *outarray = (typeout *)(out)->array->items;\
size_t best_index = 0;\
if(((op) == NUMERICAL_MAX) || ((op) == NUMERICAL_ARGMAX)) {\
for(size_t i=(start)+(stride); i < (stop); i+=(stride)) {\
if((array)[i] > (array)[best_idx]) {\
best_idx = i;\
}\
for(size_t i=1; i < (len); i++) {\
if(array[(start)+i*(increment)] > array[(start)+best_index*(increment)]) best_index = i;\
}\
if((op) == NUMERICAL_MAX) outarray[(pos)] = array[(start)+best_index*(increment)];\
else outarray[(pos)] = best_index;\
} else{\
for(size_t i=(start)+(stride); i < (stop); i+=(stride)) {\
if((array)[i] < (array)[best_idx]) best_idx = i;\
for(size_t i=1; i < (len); i++) {\
if(array[(start)+i*(increment)] < array[(start)+best_index*(increment)]) best_index = i;\
}\
if((op) == NUMERICAL_MIN) outarray[(pos)] = array[(start)+best_index*(increment)];\
else outarray[(pos)] = best_index;\
}\
} while(0)
#define RUN_SUM(ndarray, type, optype, len, start, increment) do {\
type *array = (type *)(ndarray)->array->items;\
type value;\
for(size_t j=0; j < (len); j++) {\
value = array[(start)+j*(increment)];\
sum += value;\
}\
} while(0)
#define RUN_STD(ndarray, type, len, start, increment) do {\
type *array = (type *)(ndarray)->array->items;\
mp_float_t value;\
for(size_t j=0; j < (len); j++) {\
sum += array[(start)+j*(increment)];\
}\
sum /= (len);\
for(size_t j=0; j < (len); j++) {\
value = (array[(start)+j*(increment)] - sum);\
sum_sq += value * value;\
}\
} while(0)

2
code/ulab.c
View file

@ -24,7 +24,7 @@
#include "fft.h"
#include "numerical.h"
STATIC MP_DEFINE_STR_OBJ(ulab_version_obj, "0.26.4");
STATIC MP_DEFINE_STR_OBJ(ulab_version_obj, "0.26.5");
MP_DEFINE_CONST_FUN_OBJ_1(ndarray_shape_obj, ndarray_shape);
MP_DEFINE_CONST_FUN_OBJ_1(ndarray_rawsize_obj, ndarray_rawsize);

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

@ -1,4 +1,11 @@
Tue, 7 Jan 2020
version 0.26.5
fixed glitch in numerical.c, numerical.h
Mon, 6 Jan 2020
version 0.26.4

493
docs/ulab.ipynb
View file

@ -68,8 +68,8 @@
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-05T20:18:02.700773Z",
"start_time": "2020-01-05T20:18:02.695726Z"
"end_time": "2020-01-07T20:29:29.327256Z",
"start_time": "2020-01-07T20:29:29.211644Z"
}
},
"outputs": [
@ -87,11 +87,11 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-05T20:18:30.766901Z",
"start_time": "2020-01-05T20:18:30.759404Z"
"end_time": "2020-01-07T20:29:31.701926Z",
"start_time": "2020-01-07T20:29:31.698817Z"
}
},
"outputs": [],
@ -134,11 +134,11 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-05T20:21:58.371317Z",
"start_time": "2020-01-05T20:21:58.357586Z"
"end_time": "2020-01-07T20:29:35.150872Z",
"start_time": "2020-01-07T20:29:35.124059Z"
}
},
"outputs": [],
@ -328,8 +328,8 @@
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-05T20:18:37.351688Z",
"start_time": "2020-01-05T20:18:37.343794Z"
"end_time": "2020-01-07T20:29:38.889453Z",
"start_time": "2020-01-07T20:29:38.883666Z"
}
},
"outputs": [],
@ -6965,11 +6965,11 @@
},
{
"cell_type": "code",
"execution_count": 157,
"execution_count": 7,
"metadata": {
"ExecuteTime": {
"end_time": "2019-11-06T06:20:29.556126Z",
"start_time": "2019-11-06T06:20:29.547543Z"
"end_time": "2020-01-07T20:46:49.563652Z",
"start_time": "2020-01-07T20:46:49.544742Z"
}
},
"outputs": [
@ -6977,7 +6977,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 4779 bytes to numerical.h\n"
"written 5721 bytes to numerical.h\n"
]
}
],
@ -7009,19 +7009,44 @@
"mp_obj_t numerical_sort_inplace(size_t , const mp_obj_t *, mp_map_t *);\n",
"mp_obj_t numerical_argsort(size_t , const mp_obj_t *, mp_map_t *);\n",
"\n",
"// this macro could be tighter, if we moved the ifs to the argmin function, assigned <, as well as >\n",
"#define ARG_MIN_LOOP(in, type, start, stop, stride, op) do {\\\n",
" type *array = (type *)(in)->array->items;\\\n",
"#define RUN_ARGMIN(in, out, typein, typeout, len, start, increment, op, pos) do {\\\n",
" typein *array = (typein *)(in)->array->items;\\\n",
" typeout *outarray = (typeout *)(out)->array->items;\\\n",
" size_t best_index = 0;\\\n",
" if(((op) == NUMERICAL_MAX) || ((op) == NUMERICAL_ARGMAX)) {\\\n",
" for(size_t i=(start)+(stride); i < (stop); i+=(stride)) {\\\n",
" if((array)[i] > (array)[best_idx]) {\\\n",
" best_idx = i;\\\n",
" }\\\n",
" for(size_t i=1; i < (len); i++) {\\\n",
" if(array[(start)+i*(increment)] > array[(start)+best_index*(increment)]) best_index = i;\\\n",
" }\\\n",
" if((op) == NUMERICAL_MAX) outarray[(pos)] = array[(start)+best_index*(increment)];\\\n",
" else outarray[(pos)] = best_index;\\\n",
" } else{\\\n",
" for(size_t i=(start)+(stride); i < (stop); i+=(stride)) {\\\n",
" if((array)[i] < (array)[best_idx]) best_idx = i;\\\n",
" for(size_t i=1; i < (len); i++) {\\\n",
" if(array[(start)+i*(increment)] < array[(start)+best_index*(increment)]) best_index = i;\\\n",
" }\\\n",
" if((op) == NUMERICAL_MIN) outarray[(pos)] = array[(start)+best_index*(increment)];\\\n",
" else outarray[(pos)] = best_index;\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"#define RUN_SUM(ndarray, type, optype, len, start, increment) do {\\\n",
" type *array = (type *)(ndarray)->array->items;\\\n",
" type value;\\\n",
" for(size_t j=0; j < (len); j++) {\\\n",
" value = array[(start)+j*(increment)];\\\n",
" sum += value;\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"#define RUN_STD(ndarray, type, len, start, increment) do {\\\n",
" type *array = (type *)(ndarray)->array->items;\\\n",
" mp_float_t value;\\\n",
" for(size_t j=0; j < (len); j++) {\\\n",
" sum += array[(start)+j*(increment)];\\\n",
" }\\\n",
" sum /= (len);\\\n",
" for(size_t j=0; j < (len); j++) {\\\n",
" value = (array[(start)+j*(increment)] - sum);\\\n",
" sum_sq += value * value;\\\n",
" }\\\n",
"} while(0)\n",
"\n",
@ -7119,16 +7144,17 @@
"\n",
"### Argument parsing\n",
"\n",
"\n",
"Since most of these functions operate on matrices along an axis, it might make sense to factor out the parsing of arguments and keyword arguments. The void function `numerical_parse_args` fills in the pointer for the matrix/array, and the axis."
]
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 8,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-01T19:31:51.602146Z",
"start_time": "2020-01-01T19:31:51.594464Z"
"end_time": "2020-01-07T20:49:34.247209Z",
"start_time": "2020-01-07T20:49:34.223304Z"
}
},
"outputs": [
@ -7136,13 +7162,13 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 28604 bytes to numerical.c\n"
"written 28402 bytes to numerical.c\n"
]
}
],
"source": [
"%%ccode numerical.c\n",
"\n",
" \n",
"#include <math.h>\n",
"#include <stdlib.h>\n",
"#include <string.h>\n",
@ -7212,7 +7238,33 @@
" }\n",
"}\n",
"\n",
"mp_obj_t numerical_sum_mean_std_array(mp_obj_t oin, uint8_t optype) {\n",
"void axis_sorter(ndarray_obj_t *ndarray, mp_obj_t axis, size_t *m, size_t *n, size_t *N, \n",
" size_t *increment, size_t *len, size_t *start_inc) {\n",
" if(axis == mp_const_none) { // flatten the array\n",
" *m = 1;\n",
" *n = 1;\n",
" *len = ndarray->array->len;\n",
" *N = 1;\n",
" *increment = 1;\n",
" *start_inc = ndarray->array->len;\n",
" } else if((mp_obj_get_int(axis) == 1)) { // along the horizontal axis\n",
" *m = ndarray->m;\n",
" *n = 1;\n",
" *len = ndarray->n;\n",
" *N = ndarray->m;\n",
" *increment = 1;\n",
" *start_inc = ndarray->n;\n",
" } else { // along vertical axis\n",
" *m = 1;\n",
" *n = ndarray->n;\n",
" *len = ndarray->m;\n",
" *N = ndarray->n;\n",
" *increment = ndarray->n;\n",
" *start_inc = 1;\n",
" } \n",
"}\n",
"\n",
"mp_obj_t numerical_sum_mean_std_iterable(mp_obj_t oin, uint8_t optype, size_t ddof) {\n",
" mp_float_t value, sum = 0.0, sq_sum = 0.0;\n",
" mp_obj_iter_buf_t iter_buf;\n",
" mp_obj_t item, iterable = mp_getiter(oin, &iter_buf);\n",
@ -7220,194 +7272,153 @@
" while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" value = mp_obj_get_float(item);\n",
" sum += value;\n",
" if(optype == NUMERICAL_STD) {\n",
" sq_sum += value*value;\n",
" }\n",
" }\n",
" if(optype == NUMERICAL_SUM) {\n",
" return mp_obj_new_float(sum);\n",
" } else if(optype == NUMERICAL_MEAN) {\n",
" return mp_obj_new_float(sum/len);\n",
" } else {\n",
" } else { // this should be the case of the standard deviation\n",
" // TODO: note that we could get away with a single pass, if we used the Weldorf algorithm\n",
" // That should save a fair amount of time, because we would have to extract the values only once\n",
" iterable = mp_getiter(oin, &iter_buf);\n",
" sum /= len; // this is now the mean!\n",
" return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/len-sum*sum));\n",
" }\n",
"}\n",
"\n",
"STATIC mp_float_t numerical_sum_mean_std_single_line(void *data, size_t start, size_t stop, \n",
" size_t stride, uint8_t typecode, uint8_t optype) {\n",
" \n",
" mp_float_t sum = 0.0, sq_sum = 0.0, value;\n",
" size_t len = 0;\n",
" for(size_t i=start; i < stop; i+=stride, len++) {\n",
" value = ndarray_get_float_value(data, typecode, i); \n",
" sum += value;\n",
" if(optype == NUMERICAL_STD) {\n",
" sq_sum += value*value;\n",
" }\n",
" }\n",
" if(len == 0) {\n",
" mp_raise_ValueError(\"data length is 0!\");\n",
" }\n",
" if(optype == NUMERICAL_SUM) {\n",
" return sum;\n",
" } else if(optype == NUMERICAL_MEAN) {\n",
" return sum/len;\n",
" } else {\n",
" sum /= len; // this is now the mean!\n",
" return MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/len-sum*sum);\n",
" }\n",
"}\n",
"\n",
"STATIC mp_obj_t numerical_sum_mean_std_matrix(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {\n",
" ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);\n",
" if((axis == mp_const_none) || (in->m == 1) || (in->n == 1)) { \n",
" // return the value for the flattened array\n",
" return mp_obj_new_float(numerical_sum_mean_std_single_line(in->array->items, 0, \n",
" in->array->len, 1, in->array->typecode, optype));\n",
" } else {\n",
" uint8_t _axis = mp_obj_get_int(axis);\n",
" size_t m = (_axis == 0) ? 1 : in->m;\n",
" size_t n = (_axis == 0) ? in->n : 1;\n",
" size_t len = in->array->len;\n",
" mp_float_t sms;\n",
" // TODO: pass in->array->typcode to create_new_ndarray\n",
" ndarray_obj_t *out = create_new_ndarray(m, n, NDARRAY_FLOAT);\n",
"\n",
" // TODO: these two cases could probably be combined in a more elegant fashion...\n",
" if(_axis == 0) { // vertical\n",
" for(size_t i=0; i < n; i++) {\n",
" sms = numerical_sum_mean_std_single_line(in->array->items, i, len, \n",
" n, in->array->typecode, optype);\n",
" ((float_t *)out->array->items)[i] = sms;\n",
" }\n",
" } else { // horizontal\n",
" for(size_t i=0; i < m; i++) {\n",
" sms = numerical_sum_mean_std_single_line(in->array->items, i*in->n, \n",
" (i+1)*in->n, 1, in->array->typecode, optype);\n",
" ((float_t *)out->array->items)[i] = sms;\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(out);\n",
" }\n",
"}\n",
"\n",
"size_t numerical_argmin_argmax_array(ndarray_obj_t *in, size_t start, \n",
" size_t stop, size_t stride, uint8_t op) {\n",
" size_t best_idx = start;\n",
" if(in->array->typecode == NDARRAY_UINT8) {\n",
" ARG_MIN_LOOP(in, uint8_t, start, stop, stride, op);\n",
" } else if(in->array->typecode == NDARRAY_INT8) {\n",
" ARG_MIN_LOOP(in, int8_t, start, stop, stride, op);\n",
" } else if(in->array->typecode == NDARRAY_UINT16) {\n",
" ARG_MIN_LOOP(in, uint16_t, start, stop, stride, op);\n",
" } else if(in->array->typecode == NDARRAY_INT16) {\n",
" ARG_MIN_LOOP(in, uint16_t, start, stop, stride, op);\n",
" } else if(in->array->typecode == NDARRAY_FLOAT) {\n",
" ARG_MIN_LOOP(in, mp_float_t, start, stop, stride, op);\n",
" }\n",
" return best_idx;\n",
"}\n",
"\n",
"void copy_value_into_ndarray(ndarray_obj_t *target, ndarray_obj_t *source, size_t target_idx, size_t source_idx) {\n",
" // since we are simply copying, it doesn't matter, whether the arrays are signed or unsigned, \n",
" // we can cast them in any way we like\n",
" // This could also be done with byte copies. I don't know, whether that would have any benefits\n",
" if((target->array->typecode == NDARRAY_UINT8) || (target->array->typecode == NDARRAY_INT8)) {\n",
" ((uint8_t *)target->array->items)[target_idx] = ((uint8_t *)source->array->items)[source_idx];\n",
" } else if((target->array->typecode == NDARRAY_UINT16) || (target->array->typecode == NDARRAY_INT16)) {\n",
" ((uint16_t *)target->array->items)[target_idx] = ((uint16_t *)source->array->items)[source_idx];\n",
" } else { \n",
" ((float *)target->array->items)[target_idx] = ((float *)source->array->items)[source_idx];\n",
" }\n",
"}\n",
" \n",
"STATIC mp_obj_t numerical_argmin_argmax(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {\n",
" if(mp_obj_is_type(oin, &mp_type_tuple) || mp_obj_is_type(oin, &mp_type_list) || \n",
" mp_obj_is_type(oin, &mp_type_range)) {\n",
" // This case will work for single iterables only \n",
" size_t idx = 0, best_idx = 0;\n",
" mp_obj_iter_buf_t iter_buf;\n",
" mp_obj_t iterable = mp_getiter(oin, &iter_buf);\n",
" mp_obj_t best_obj = MP_OBJ_NULL;\n",
" mp_obj_t item;\n",
" mp_uint_t op = MP_BINARY_OP_LESS;\n",
" if((optype == NUMERICAL_ARGMAX) || (optype == NUMERICAL_MAX)) op = MP_BINARY_OP_MORE;\n",
" while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if ((best_obj == MP_OBJ_NULL) || (mp_binary_op(op, item, best_obj) == mp_const_true)) {\n",
" best_obj = item;\n",
" best_idx = idx;\n",
" }\n",
" idx++;\n",
" value = mp_obj_get_float(item) - sum;\n",
" sq_sum += value * value;\n",
" }\n",
" if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {\n",
" return MP_OBJ_NEW_SMALL_INT(best_idx);\n",
" } else {\n",
" return best_obj;\n",
" }\n",
" } else if(mp_obj_is_type(oin, &ulab_ndarray_type)) {\n",
" ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);\n",
" size_t best_idx;\n",
" if((axis == mp_const_none) || (in->m == 1) || (in->n == 1)) {\n",
" // return the value for the flattened array \n",
" best_idx = numerical_argmin_argmax_array(in, 0, in->array->len, 1, optype);\n",
" if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {\n",
" return MP_OBJ_NEW_SMALL_INT(best_idx);\n",
" } else {\n",
" if(in->array->typecode == NDARRAY_FLOAT) {\n",
" return mp_obj_new_float(ndarray_get_float_value(in->array->items, in->array->typecode, best_idx));\n",
" } else {\n",
" return mp_binary_get_val_array(in->array->typecode, in->array->items, best_idx);\n",
" }\n",
" }\n",
" } else { // we have to work with a full matrix here\n",
" uint8_t _axis = mp_obj_get_int(axis);\n",
" size_t m = (_axis == 0) ? 1 : in->m;\n",
" size_t n = (_axis == 0) ? in->n : 1;\n",
" size_t len = in->array->len;\n",
" ndarray_obj_t *ndarray = NULL;\n",
" if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {\n",
" ndarray = create_new_ndarray(m, n, in->array->typecode);\n",
" } else { // argmin/argmax\n",
" // TODO: one might get away with uint8_t, if both m, and n < 255\n",
" ndarray = create_new_ndarray(m, n, NDARRAY_UINT16);\n",
" }\n",
"\n",
" // TODO: these two cases could probably be combined in a more elegant fashion...\n",
" if(_axis == 0) { // vertical\n",
" for(size_t i=0; i < n; i++) {\n",
" best_idx = numerical_argmin_argmax_array(in, i, len, n, optype);\n",
" if((optype == NUMERICAL_MIN) || (optype == NUMERICAL_MAX)) {\n",
" copy_value_into_ndarray(ndarray, in, i, best_idx);\n",
" } else {\n",
" ((uint16_t *)ndarray->array->items)[i] = (uint16_t)(best_idx / n);\n",
" }\n",
" }\n",
" } else { // horizontal\n",
" for(size_t i=0; i < m; i++) {\n",
" best_idx = numerical_argmin_argmax_array(in, i*in->n, (i+1)*in->n, 1, optype);\n",
" if((optype == NUMERICAL_MIN) || (optype == NUMERICAL_MAX)) {\n",
" copy_value_into_ndarray(ndarray, in, i, best_idx);\n",
" } else {\n",
" ((uint16_t *)ndarray->array->items)[i] = (uint16_t)(best_idx - i*in->n);\n",
" }\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
" }\n",
" return mp_const_none;\n",
" }\n",
" mp_raise_TypeError(\"input type is not supported\");\n",
" return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/(len-ddof)));\n",
" }\n",
"}\n",
"\n",
"STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t type) {\n",
"STATIC mp_obj_t numerical_sum_mean_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {\n",
" size_t m, n, increment, start, start_inc, N, len; \n",
" axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);\n",
" ndarray_obj_t *results = create_new_ndarray(m, n, NDARRAY_FLOAT);\n",
" mp_float_t sum, sq_sum;\n",
" mp_float_t *farray = (mp_float_t *)results->array->items;\n",
" for(size_t j=0; j < N; j++) { // result index\n",
" start = j * start_inc;\n",
" sum = sq_sum = 0.0;\n",
" if(ndarray->array->typecode == NDARRAY_UINT8) {\n",
" RUN_SUM(ndarray, uint8_t, optype, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_INT8) {\n",
" RUN_SUM(ndarray, int8_t, optype, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_UINT16) {\n",
" RUN_SUM(ndarray, uint16_t, optype, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_INT16) {\n",
" RUN_SUM(ndarray, int16_t, optype, len, start, increment);\n",
" } else { // this will be mp_float_t, no need to check\n",
" RUN_SUM(ndarray, mp_float_t, optype, len, start, increment);\n",
" }\n",
" if(optype == NUMERICAL_SUM) {\n",
" farray[j] = sum;\n",
" } else { // this is the case of the mean\n",
" farray[j] = sum / len;\n",
" }\n",
" }\n",
" if(results->array->len == 1) {\n",
" return mp_obj_new_float(farray[0]);\n",
" }\n",
" return MP_OBJ_FROM_PTR(results);\n",
"}\n",
"\n",
"mp_obj_t numerical_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, size_t ddof) {\n",
" size_t m, n, increment, start, start_inc, N, len; \n",
" mp_float_t sum, sum_sq;\n",
" \n",
" axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);\n",
" if(ddof > len) {\n",
" mp_raise_ValueError(\"ddof must be smaller than length of data set\");\n",
" }\n",
" ndarray_obj_t *results = create_new_ndarray(m, n, NDARRAY_FLOAT);\n",
" mp_float_t *farray = (mp_float_t *)results->array->items;\n",
" for(size_t j=0; j < N; j++) { // result index\n",
" start = j * start_inc;\n",
" sum = 0.0;\n",
" sum_sq = 0.0;\n",
" if(ndarray->array->typecode == NDARRAY_UINT8) {\n",
" RUN_STD(ndarray, uint8_t, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_INT8) {\n",
" RUN_STD(ndarray, int8_t, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_UINT16) {\n",
" RUN_STD(ndarray, uint16_t, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_INT16) {\n",
" RUN_STD(ndarray, int16_t, len, start, increment);\n",
" } else { // this will be mp_float_t, no need to check\n",
" RUN_STD(ndarray, mp_float_t, len, start, increment);\n",
" }\n",
" farray[j] = MICROPY_FLOAT_C_FUN(sqrt)(sum_sq/(len - ddof));\n",
" }\n",
" if(results->array->len == 1) {\n",
" return mp_obj_new_float(farray[0]);\n",
" }\n",
" return MP_OBJ_FROM_PTR(results);\n",
"}\n",
"\n",
"mp_obj_t numerical_argmin_argmax_iterable(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {\n",
" size_t idx = 0, best_idx = 0;\n",
" mp_obj_iter_buf_t iter_buf;\n",
" mp_obj_t iterable = mp_getiter(oin, &iter_buf);\n",
" mp_obj_t best_obj = MP_OBJ_NULL;\n",
" mp_obj_t item;\n",
" mp_uint_t op = MP_BINARY_OP_LESS;\n",
" if((optype == NUMERICAL_ARGMAX) || (optype == NUMERICAL_MAX)) op = MP_BINARY_OP_MORE;\n",
" while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if ((best_obj == MP_OBJ_NULL) || (mp_binary_op(op, item, best_obj) == mp_const_true)) {\n",
" best_obj = item;\n",
" best_idx = idx;\n",
" }\n",
" idx++;\n",
" }\n",
" if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {\n",
" return MP_OBJ_NEW_SMALL_INT(best_idx);\n",
" } else {\n",
" return best_obj;\n",
" } \n",
"}\n",
"\n",
"mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {\n",
" size_t m, n, increment, start, start_inc, N, len;\n",
" axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);\n",
" ndarray_obj_t *results;\n",
" if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {\n",
" // we could save some RAM by taking NDARRAY_UINT8, if the dimensions \n",
" // are smaller than 256, but the code would become more involving \n",
" // (we would also need extra flash space)\n",
" results = create_new_ndarray(m, n, NDARRAY_UINT16);\n",
" } else {\n",
" results = create_new_ndarray(m, n, ndarray->array->typecode);\n",
" }\n",
" \n",
" for(size_t j=0; j < N; j++) { // result index\n",
" start = j * start_inc;\n",
" if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {\n",
" if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {\n",
" RUN_ARGMIN(ndarray, results, uint8_t, uint8_t, len, start, increment, optype, j);\n",
" } else {\n",
" RUN_ARGMIN(ndarray, results, uint8_t, uint16_t, len, start, increment, optype, j); \n",
" }\n",
" } else if((ndarray->array->typecode == NDARRAY_UINT16) || (ndarray->array->typecode == NDARRAY_INT16)) {\n",
" RUN_ARGMIN(ndarray, results, uint16_t, uint16_t, len, start, increment, optype, j);\n",
" } else {\n",
" if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {\n",
" RUN_ARGMIN(ndarray, results, mp_float_t, mp_float_t, len, start, increment, optype, j);\n",
" } else {\n",
" RUN_ARGMIN(ndarray, results, mp_float_t, uint16_t, len, start, increment, optype, j); \n",
" }\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(results);\n",
"}\n",
"\n",
"STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t optype) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} } ,\n",
" { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },\n",
" { MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" \n",
" mp_obj_t oin = args[0].u_obj;\n",
" mp_obj_t axis = args[1].u_obj;\n",
@ -7416,32 +7427,31 @@
" mp_raise_ValueError(\"axis must be None, 0, or 1\");\n",
" }\n",
" \n",
" if(mp_obj_is_type(oin, &mp_type_tuple) || mp_obj_is_type(oin, &mp_type_list) || \n",
" mp_obj_is_type(oin, &mp_type_range)) {\n",
" switch(type) {\n",
" if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) || \n",
" MP_OBJ_IS_TYPE(oin, &mp_type_range)) {\n",
" switch(optype) {\n",
" case NUMERICAL_MIN:\n",
" case NUMERICAL_ARGMIN:\n",
" case NUMERICAL_MAX:\n",
" case NUMERICAL_ARGMAX:\n",
" return numerical_argmin_argmax(oin, axis, type);\n",
" return numerical_argmin_argmax_iterable(oin, axis, optype);\n",
" case NUMERICAL_SUM:\n",
" case NUMERICAL_MEAN:\n",
" case NUMERICAL_STD:\n",
" return numerical_sum_mean_std_array(oin, type);\n",
" return numerical_sum_mean_std_iterable(oin, optype, 0);\n",
" default: // we should never reach this point, but whatever\n",
" return mp_const_none;\n",
" }\n",
" } else if(mp_obj_is_type(oin, &ulab_ndarray_type)) {\n",
" switch(type) {\n",
" } else if(MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {\n",
" ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);\n",
" switch(optype) {\n",
" case NUMERICAL_MIN:\n",
" case NUMERICAL_MAX:\n",
" case NUMERICAL_ARGMIN:\n",
" case NUMERICAL_ARGMAX:\n",
" return numerical_argmin_argmax(oin, axis, type);\n",
" return numerical_argmin_argmax_ndarray(ndarray, axis, optype);\n",
" case NUMERICAL_SUM:\n",
" case NUMERICAL_MEAN:\n",
" case NUMERICAL_STD:\n",
" return numerical_sum_mean_std_matrix(oin, axis, type); \n",
" return numerical_sum_mean_ndarray(ndarray, axis, optype);\n",
" default:\n",
" mp_raise_NotImplementedError(\"operation is not implemented on ndarrays\");\n",
" }\n",
@ -7476,7 +7486,31 @@
"}\n",
"\n",
"mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" return numerical_function(n_args, pos_args, kw_args, NUMERICAL_STD);\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} } ,\n",
" { MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },\n",
" { MP_QSTR_ddof, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" \n",
" mp_obj_t oin = args[0].u_obj;\n",
" mp_obj_t axis = args[1].u_obj;\n",
" size_t ddof = args[2].u_int;\n",
" if((axis != mp_const_none) && (mp_obj_get_int(axis) != 0) && (mp_obj_get_int(axis) != 1)) {\n",
" // this seems to pass with False, and True...\n",
" mp_raise_ValueError(\"axis must be None, 0, or 1\");\n",
" }\n",
" if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) || MP_OBJ_IS_TYPE(oin, &mp_type_range)) {\n",
" return numerical_sum_mean_std_iterable(oin, NUMERICAL_STD, ddof);\n",
" } else if(MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {\n",
" ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);\n",
" return numerical_std_ndarray(ndarray, axis, ddof);\n",
" } else {\n",
" mp_raise_TypeError(\"input must be tuple, list, range, or ndarray\");\n",
" }\n",
" return mp_const_none;\n",
"}\n",
"\n",
"mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
@ -7844,11 +7878,11 @@
},
{
"cell_type": "code",
"execution_count": 105,
"execution_count": 12,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-07T19:59:33.514774Z",
"start_time": "2020-01-07T19:59:33.506709Z"
"end_time": "2020-01-07T20:51:32.722886Z",
"start_time": "2020-01-07T20:51:32.713985Z"
}
},
"outputs": [
@ -7879,7 +7913,7 @@
"#include \"fft.h\"\n",
"#include \"numerical.h\"\n",
"\n",
"STATIC MP_DEFINE_STR_OBJ(ulab_version_obj, \"0.26.4\");\n",
"STATIC MP_DEFINE_STR_OBJ(ulab_version_obj, \"0.26.5\");\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_1(ndarray_shape_obj, ndarray_shape);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(ndarray_rawsize_obj, ndarray_rawsize);\n",
@ -8123,11 +8157,11 @@
},
{
"cell_type": "code",
"execution_count": 102,
"execution_count": 9,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-07T19:53:37.999138Z",
"start_time": "2020-01-07T19:53:37.196162Z"
"end_time": "2020-01-07T20:49:53.299366Z",
"start_time": "2020-01-07T20:49:51.912441Z"
}
},
"outputs": [
@ -8151,8 +8185,8 @@
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-07T19:59:39.979093Z",
"start_time": "2020-01-07T19:59:38.833117Z"
"end_time": "2020-01-07T20:50:29.269290Z",
"start_time": "2020-01-07T20:49:55.946789Z"
},
"scrolled": false
},
@ -8214,11 +8248,11 @@
},
{
"cell_type": "code",
"execution_count": 108,
"execution_count": 13,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-07T20:00:33.493890Z",
"start_time": "2020-01-07T20:00:33.436113Z"
"end_time": "2020-01-07T20:52:26.547974Z",
"start_time": "2020-01-07T20:52:26.542000Z"
}
},
"outputs": [
@ -8233,6 +8267,13 @@
"source": [
"%%writefile ../../../ulab/docs/ulab-change-log.md\n",
"\n",
"Tue, 7 Jan 2020\n",
"\n",
"version 0.26.5\n",
"\n",
" fixed glitch in numerical.c, numerical.h\n",
"\n",
"Mon, 6 Jan 2020\n",
"\n",
"version 0.26.4\n",
"\n",