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

replaced float with mp_float_t everywhere, so that the source can be compiled for all platforms

Browse source
This commit is contained in:
Zoltán Vörös 2019-11-01 14:28:29 +01:00
parent 8b6d03417b
commit 6e4911d96f
12 changed files with 566 additions and 307 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

20
code/fft.c
View file

@ -26,13 +26,13 @@ enum FFT_TYPE {
FFT_SPECTRUM,
};
void fft_kernel(float *real, float *imag, int n, int isign) {
void fft_kernel(mp_float_t *real, mp_float_t *imag, int n, int isign) {
// This is basically a modification of four1 from Numerical Recipes
// The main difference is that this function takes two arrays, one
// for the real, and one for the imaginary parts.
int j, m, mmax, istep;
float tempr, tempi;
float wtemp, wr, wpr, wpi, wi, theta;
mp_float_t tempr, tempi;
mp_float_t wtemp, wr, wpr, wpi, wi, theta;
j = 0;
for(int i = 0; i < n; i++) {
@ -52,9 +52,9 @@ void fft_kernel(float *real, float *imag, int n, int isign) {
while (n > mmax) {
istep = mmax << 1;
theta = -1.0*isign*6.28318530717959/istep;
wtemp = sinf(0.5 * theta);
wtemp = sin(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = sinf(theta);
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for(m = 0; m < mmax; m++) {
@ -92,19 +92,19 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
}
ndarray_obj_t *out_re = create_new_ndarray(1, len, NDARRAY_FLOAT);
float *data_re = (float *)out_re->array->items;
mp_float_t *data_re = (mp_float_t *)out_re->array->items;
if(re->array->typecode == NDARRAY_FLOAT) {
// By treating this case separately, we can save a bit of time.
// I don't know if it is worthwhile, though...
memcpy((float *)out_re->array->items, (float *)re->array->items, re->bytes);
memcpy((mp_float_t *)out_re->array->items, (mp_float_t *)re->array->items, re->bytes);
} else {
for(size_t i=0; i < len; i++) {
data_re[i] = ndarray_get_float_value(re->array->items, re->array->typecode, i);
}
}
ndarray_obj_t *out_im = create_new_ndarray(1, len, NDARRAY_FLOAT);
float *data_im = (float *)out_im->array->items;
mp_float_t *data_im = (mp_float_t *)out_im->array->items;
if(n_args == 2) {
ndarray_obj_t *im = MP_OBJ_TO_PTR(arg_im);
@ -112,7 +112,7 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
mp_raise_ValueError("real and imaginary parts must be of equal length");
}
if(im->array->typecode == NDARRAY_FLOAT) {
memcpy((float *)out_im->array->items, (float *)im->array->items, im->bytes);
memcpy((mp_float_t *)out_im->array->items, (mp_float_t *)im->array->items, im->bytes);
} else {
for(size_t i=0; i < len; i++) {
data_im[i] = ndarray_get_float_value(im->array->items, im->array->typecode, i);
@ -123,7 +123,7 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
fft_kernel(data_re, data_im, len, 1);
if(type == FFT_SPECTRUM) {
for(size_t i=0; i < len; i++) {
data_re[i] = sqrtf(data_re[i]*data_re[i] + data_im[i]*data_im[i]);
data_re[i] = sqrt(data_re[i]*data_re[i] + data_im[i]*data_im[i]);
}
}
} else { // inverse transform

56
code/linalg.c
View file

@ -108,24 +108,24 @@ mp_obj_t linalg_size(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args)
}
}
bool linalg_invert_matrix(float *data, size_t N) {
bool linalg_invert_matrix(mp_float_t *data, size_t N) {
// returns true, of the inversion was successful,
// false, if the matrix is singular
// initially, this is the unit matrix: the contents of this matrix is what
// will be returned after all the transformations
float *unit = m_new(float, N*N);
mp_float_t *unit = m_new(mp_float_t, N*N);
float elem = 1.0;
mp_float_t elem = 1.0;
// initialise the unit matrix
memset(unit, 0, sizeof(float)*N*N);
memset(unit, 0, sizeof(mp_float_t)*N*N);
for(size_t m=0; m < N; m++) {
memcpy(&unit[m*(N+1)], &elem, sizeof(float));
memcpy(&unit[m*(N+1)], &elem, sizeof(mp_float_t));
}
for(size_t m=0; m < N; m++){
// this could be faster with ((c < epsilon) && (c > -epsilon))
if(abs(data[m*(N+1)]) < epsilon) {
m_del(float, unit, N*N);
m_del(mp_float_t, unit, N*N);
return false;
}
for(size_t n=0; n < N; n++){
@ -145,8 +145,8 @@ bool linalg_invert_matrix(float *data, size_t N) {
unit[N*m+n] /= elem;
}
}
memcpy(data, unit, sizeof(float)*N*N);
m_del(float, unit, N*N);
memcpy(data, unit, sizeof(mp_float_t)*N*N);
m_del(mp_float_t, unit, N*N);
return true;
}
@ -163,21 +163,21 @@ mp_obj_t linalg_inv(mp_obj_t o_in) {
mp_raise_ValueError("only square matrices can be inverted");
}
ndarray_obj_t *inverted = create_new_ndarray(o->m, o->n, NDARRAY_FLOAT);
float *data = (float *)inverted->array->items;
mp_float_t *data = (mp_float_t *)inverted->array->items;
mp_obj_t elem;
for(size_t m=0; m < o->m; m++) { // rows first
for(size_t n=0; n < o->n; n++) { // columns next
// this could, perhaps, be done in single line...
// On the other hand, we probably spend little time here
elem = mp_binary_get_val_array(o->array->typecode, o->array->items, m*o->n+n);
data[m*o->n+n] = (float)mp_obj_get_float(elem);
data[m*o->n+n] = (mp_float_t)mp_obj_get_float(elem);
}
}
if(!linalg_invert_matrix(data, o->m)) {
// TODO: I am not sure this is needed here. Otherwise,
// how should we free up the unused RAM of inverted?
m_del(float, inverted->array->items, o->n*o->n);
m_del(mp_float_t, inverted->array->items, o->n*o->n);
mp_raise_ValueError("input matrix is singular");
}
return MP_OBJ_FROM_PTR(inverted);
@ -192,8 +192,8 @@ mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) {
}
// TODO: numpy uses upcasting here
ndarray_obj_t *out = create_new_ndarray(m1->m, m2->n, NDARRAY_FLOAT);
float *outdata = (float *)out->array->items;
float sum, v1, v2;
mp_float_t *outdata = (mp_float_t *)out->array->items;
mp_float_t sum, v1, v2;
for(size_t i=0; i < m1->n; i++) {
for(size_t j=0; j < m2->m; j++) {
sum = 0.0;
@ -301,14 +301,14 @@ mp_obj_t linalg_det(mp_obj_t oin) {
mp_raise_ValueError("input must be square matrix");
}
float *tmp = m_new(float, in->n*in->n);
mp_float_t *tmp = m_new(mp_float_t, in->n*in->n);
for(size_t i=0; i < in->array->len; i++){
tmp[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);
}
float c;
mp_float_t c;
for(size_t m=0; m < in->m-1; m++){
if(abs(tmp[m*(in->n+1)]) < epsilon) {
m_del(float, tmp, in->n*in->n);
m_del(mp_float_t, tmp, in->n*in->n);
return mp_obj_new_float(0.0);
}
for(size_t n=0; n < in->n; n++){
@ -320,12 +320,12 @@ mp_obj_t linalg_det(mp_obj_t oin) {
}
}
}
float det = 1.0;
mp_float_t det = 1.0;
for(size_t m=0; m < in->m; m++){
det *= tmp[m*(in->n+1)];
}
m_del(float, tmp, in->n*in->n);
m_del(mp_float_t, tmp, in->n*in->n);
return mp_obj_new_float(det);
}
@ -337,7 +337,7 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
if(in->m != in->n) {
mp_raise_ValueError("input must be square matrix");
}
float *array = m_new(float, in->array->len);
mp_float_t *array = m_new(mp_float_t, in->array->len);
for(size_t i=0; i < in->array->len; i++) {
array[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);
}
@ -354,12 +354,12 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
// if we got this far, then the matrix will be symmetric
ndarray_obj_t *eigenvectors = create_new_ndarray(in->m, in->n, NDARRAY_FLOAT);
float *eigvectors = (float *)eigenvectors->array->items;
mp_float_t *eigvectors = (mp_float_t *)eigenvectors->array->items;
// start out with the unit matrix
for(size_t m=0; m < in->m; m++) {
eigvectors[m*(in->n+1)] = 1.0;
}
float largest, w, t, c, s, tau, aMk, aNk, vm, vn;
mp_float_t largest, w, t, c, s, tau, aMk, aNk, vm, vn;
size_t M, N;
size_t iterations = JACOBI_MAX*in->n*in->n;
do {
@ -387,12 +387,12 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
// The following if/else chooses the smaller absolute value for the tangent
// of the rotation angle. Going with the smaller should be numerically stabler.
if(w > 0) {
t = sqrtf(w*w + 1.0) - w;
t = sqrt(w*w + 1.0) - w;
} else {
t = -1.0*(sqrtf(w*w + 1.0) + w);
t = -1.0*(sqrt(w*w + 1.0) + w);
}
s = t / sqrtf(t*t + 1.0); // the sine of the rotation angle
c = 1.0 / sqrtf(t*t + 1.0); // the cosine of the rotation angle
s = t / sqrt(t*t + 1.0); // the sine of the rotation angle
c = 1.0 / sqrt(t*t + 1.0); // the cosine of the rotation angle
tau = (1.0-c)/s; // this is equal to the tangent of the half of the rotation angle
// at this point, we have the rotation angles, so we can transform the matrix
@ -438,15 +438,15 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
if(iterations == 0) {
// the computation did not converge; numpy raises LinAlgError
m_del(float, array, in->array->len);
m_del(mp_float_t, array, in->array->len);
mp_raise_ValueError("iterations did not converge");
}
ndarray_obj_t *eigenvalues = create_new_ndarray(1, in->n, NDARRAY_FLOAT);
float *eigvalues = (float *)eigenvalues->array->items;
mp_float_t *eigvalues = (mp_float_t *)eigenvalues->array->items;
for(size_t i=0; i < in->n; i++) {
eigvalues[i] = array[i*(in->n+1)];
}
m_del(float, array, in->array->len);
m_del(mp_float_t, array, in->array->len);
mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(mp_obj_new_tuple(2, NULL));
tuple->items[0] = MP_OBJ_FROM_PTR(eigenvalues);

10
code/linalg.h
View file

@ -14,13 +14,19 @@
#include "ndarray.h"
#define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }
#define epsilon 1e-6
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
#define epsilon 1.2e-7
#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
#define epsilon 2.3e-16
#endif
#define JACOBI_MAX 20
mp_obj_t linalg_transpose(mp_obj_t );
mp_obj_t linalg_reshape(mp_obj_t , mp_obj_t );
mp_obj_t linalg_size(size_t , const mp_obj_t *, mp_map_t *);
bool linalg_invert_matrix(float *, size_t );
bool linalg_invert_matrix(mp_float_t *, size_t );
mp_obj_t linalg_inv(mp_obj_t );
mp_obj_t linalg_dot(mp_obj_t , mp_obj_t );
mp_obj_t linalg_zeros(size_t , const mp_obj_t *, mp_map_t *);

85
code/ndarray.c
View file

@ -37,26 +37,26 @@ STATIC mp_obj_array_t *array_new(char typecode, size_t n) {
return o;
}
float ndarray_get_float_value(void *data, uint8_t typecode, size_t index) {
mp_float_t ndarray_get_float_value(void *data, uint8_t typecode, size_t index) {
if(typecode == NDARRAY_UINT8) {
return (float)((uint8_t *)data)[index];
return (mp_float_t)((uint8_t *)data)[index];
} else if(typecode == NDARRAY_INT8) {
return (float)((int8_t *)data)[index];
return (mp_float_t)((int8_t *)data)[index];
} else if(typecode == NDARRAY_UINT16) {
return (float)((uint16_t *)data)[index];
return (mp_float_t)((uint16_t *)data)[index];
} else if(typecode == NDARRAY_INT16) {
return (float)((int16_t *)data)[index];
return (mp_float_t)((int16_t *)data)[index];
} else {
return (float)((float_t *)data)[index];
return (mp_float_t)((float_t *)data)[index];
}
}
void fill_array_iterable(float *array, mp_obj_t iterable) {
void fill_array_iterable(mp_float_t *array, mp_obj_t iterable) {
mp_obj_iter_buf_t x_buf;
mp_obj_t x_item, x_iterable = mp_getiter(iterable, &x_buf);
size_t i=0;
while ((x_item = mp_iternext(x_iterable)) != MP_OBJ_STOP_ITERATION) {
array[i] = (float)mp_obj_get_float(x_item);
array[i] = (mp_float_t)mp_obj_get_float(x_item);
i++;
}
}
@ -90,17 +90,22 @@ void ndarray_print(const mp_print_t *print, mp_obj_t self_in, mp_print_kind_t ki
(void)kind;
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
mp_print_str(print, "array(");
if((self->m == 1) || (self->n == 1)) {
ndarray_print_row(print, self->array, 0, self->array->len);
if(self->array->len == 0) {
mp_print_str(print, "[]");
} else {
// TODO: add vertical ellipses for the case, when self->m > PRINT_MAX
mp_print_str(print, "[");
ndarray_print_row(print, self->array, 0, self->n);
for(size_t i=1; i < self->m; i++) {
mp_print_str(print, ",\n\t ");
ndarray_print_row(print, self->array, i*self->n, self->n);
if((self->m == 1) || (self->n == 1)) {
ndarray_print_row(print, self->array, 0, self->array->len);
} else {
// TODO: add vertical ellipses for the case, when self->m > PRINT_MAX
mp_print_str(print, "[");
ndarray_print_row(print, self->array, 0, self->n);
for(size_t i=1; i < self->m; i++) {
mp_print_str(print, ",\n\t ");
ndarray_print_row(print, self->array, i*self->n, self->n);
}
mp_print_str(print, "]");
}
mp_print_str(print, "]");
}
if(self->array->typecode == NDARRAY_UINT8) {
mp_print_str(print, ", dtype=uint8)");
@ -112,7 +117,7 @@ void ndarray_print(const mp_print_t *print, mp_obj_t self_in, mp_print_kind_t ki
mp_print_str(print, ", dtype=int16)");
} else if(self->array->typecode == NDARRAY_FLOAT) {
mp_print_str(print, ", dtype=float)");
}
}
}
void ndarray_assign_elements(mp_obj_array_t *data, mp_obj_t iterable, uint8_t typecode, size_t *idx) {
@ -681,6 +686,9 @@ mp_obj_t ndarray_asbytearray(mp_obj_t self_in) {
// Binary operations
mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) {
// if(op == MP_BINARY_OP_REVERSE_ADD) {
// return ndarray_binary_op(MP_BINARY_OP_ADD, rhs, lhs);
// }
// One of the operands is a scalar
// TODO: conform to numpy with the upcasting
// TODO: implement in-place operators
@ -690,19 +698,18 @@ mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) {
int32_t ivalue = mp_obj_get_int(rhs);
if((ivalue > 0) && (ivalue < 256)) {
CREATE_SINGLE_ITEM(RHS, uint8_t, NDARRAY_UINT8, ivalue);
}
else if((ivalue > 0) && (ivalue < 65535)) {
} else if((ivalue > 255) && (ivalue < 65535)) {
CREATE_SINGLE_ITEM(RHS, uint16_t, NDARRAY_UINT16, ivalue);
}
else if((ivalue < 0) && (ivalue > -128)) {
} else if((ivalue < 0) && (ivalue > -128)) {
CREATE_SINGLE_ITEM(RHS, int8_t, NDARRAY_INT8, ivalue);
}
else if((ivalue < 0) && (ivalue > -32767)) {
} else if((ivalue < -127) && (ivalue > -32767)) {
CREATE_SINGLE_ITEM(RHS, int16_t, NDARRAY_INT16, ivalue);
} else { // the integer value clearly does not fit the ulab types, so move on to float
CREATE_SINGLE_ITEM(RHS, mp_float_t, NDARRAY_FLOAT, ivalue);
}
} else if(mp_obj_is_float(rhs)) {
float fvalue = mp_obj_get_float(rhs);
CREATE_SINGLE_ITEM(RHS, float, NDARRAY_FLOAT, fvalue);
mp_float_t fvalue = mp_obj_get_float(rhs);
CREATE_SINGLE_ITEM(RHS, mp_float_t, NDARRAY_FLOAT, fvalue);
} else {
RHS = rhs;
rhs_is_scalar = false;
@ -765,7 +772,7 @@ mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) {
} else if(or->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int16_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint8_t, float, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ol, or, op);
}
} else if(ol->array->typecode == NDARRAY_INT8) {
if(or->array->typecode == NDARRAY_UINT8) {
@ -777,7 +784,7 @@ mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) {
} else if(or->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, int16_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, int8_t, float, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int8_t, mp_float_t, ol, or, op);
}
} else if(ol->array->typecode == NDARRAY_UINT16) {
if(or->array->typecode == NDARRAY_UINT8) {
@ -787,9 +794,9 @@ mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) {
} else if(or->array->typecode == NDARRAY_UINT16) {
RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint16_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint16_t, int16_t, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, int16_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint8_t, float, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ol, or, op);
}
} else if(ol->array->typecode == NDARRAY_INT16) {
if(or->array->typecode == NDARRAY_UINT8) {
@ -797,23 +804,23 @@ mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) {
} else if(or->array->typecode == NDARRAY_INT8) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int8_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_UINT16) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, int16_t, uint16_t, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int16_t, uint16_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int16_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint16_t, float, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, mp_float_t, ol, or, op);
}
} else if(ol->array->typecode == NDARRAY_FLOAT) {
if(or->array->typecode == NDARRAY_UINT8) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, uint8_t, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint8_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_INT8) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, int8_t, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int8_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_UINT16) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, uint16_t, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint16_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, int16_t, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int16_t, ol, or, op);
} else if(or->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, float, ol, or, op);
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, mp_float_t, ol, or, op);
}
} else { // this should never happen
mp_raise_TypeError("wrong input type");
@ -867,7 +874,7 @@ mp_obj_t ndarray_unary_op(mp_unary_op_t op, mp_obj_t self_in) {
int16_t *array = (int16_t *)ndarray->array->items;
for(size_t i=0; i < self->array->len; i++) array[i] = -array[i];
} else {
float *array = (float *)ndarray->array->items;
mp_float_t *array = (mp_float_t *)ndarray->array->items;
for(size_t i=0; i < self->array->len; i++) array[i] = -array[i];
}
return MP_OBJ_FROM_PTR(ndarray);
@ -892,7 +899,7 @@ mp_obj_t ndarray_unary_op(mp_unary_op_t op, mp_obj_t self_in) {
if(array[i] < 0) array[i] = -array[i];
}
} else {
float *array = (float *)ndarray->array->items;
mp_float_t *array = (mp_float_t *)ndarray->array->items;
for(size_t i=0; i < self->array->len; i++) {
if(array[i] < 0) array[i] = -array[i];
}

15
code/ndarray.h
View file

@ -18,6 +18,11 @@
#define PRINT_MAX 10
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
#define FLOAT_TYPECODE 'f'
#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
#define FLOAT_TYPECODE 'd'
#endif
const mp_obj_type_t ulab_ndarray_type;
@ -26,7 +31,7 @@ enum NDARRAY_TYPE {
NDARRAY_INT8 = 'b',
NDARRAY_UINT16 = 'H',
NDARRAY_INT16 = 'h',
NDARRAY_FLOAT = 'f',
NDARRAY_FLOAT = FLOAT_TYPECODE,
};
typedef struct _ndarray_obj_t {
@ -39,8 +44,8 @@ typedef struct _ndarray_obj_t {
mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t , size_t , mp_obj_iter_buf_t *);
float ndarray_get_float_value(void *, uint8_t , size_t );
void fill_array_iterable(float *, mp_obj_t );
mp_float_t ndarray_get_float_value(void *, uint8_t , size_t );
void fill_array_iterable(mp_float_t *, mp_obj_t );
void ndarray_print_row(const mp_print_t *, mp_obj_array_t *, size_t , size_t );
void ndarray_print(const mp_print_t *, mp_obj_t , mp_print_kind_t );
@ -87,8 +92,8 @@ mp_obj_t ndarray_asbytearray(mp_obj_t );
return MP_OBJ_FROM_PTR(out);\
} else if((op) == MP_BINARY_OP_TRUE_DIVIDE) {\
ndarray_obj_t *out = create_new_ndarray(ol->m, ol->n, NDARRAY_FLOAT);\
float *odata = (float *)out->array->items;\
for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = (float)left[i]/(float)right[j];\
mp_float_t *odata = (mp_float_t *)out->array->items;\
for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = (mp_float_t)left[i]/(mp_float_t)right[j];\
return MP_OBJ_FROM_PTR(out);\
} else if(((op) == MP_BINARY_OP_LESS) || ((op) == MP_BINARY_OP_LESS_EQUAL) || \
((op) == MP_BINARY_OP_MORE) || ((op) == MP_BINARY_OP_MORE_EQUAL)) {\

75
code/numerical.c
View file

@ -44,7 +44,7 @@ mp_obj_t numerical_linspace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *k
if(len < 2) {
mp_raise_ValueError("number of points must be at least 2");
}
float value, step;
mp_float_t value, step;
value = mp_obj_get_float(args[0].u_obj);
uint8_t typecode = args[5].u_int;
if(args[3].u_obj == mp_const_true) step = (mp_obj_get_float(args[1].u_obj)-value)/(len-1);
@ -63,7 +63,7 @@ mp_obj_t numerical_linspace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *k
int16_t *array = (int16_t *)ndarray->array->items;
for(size_t i=0; i < len; i++, value += step) array[i] = (int16_t)value;
} else {
float *array = (float *)ndarray->array->items;
mp_float_t *array = (mp_float_t *)ndarray->array->items;
for(size_t i=0; i < len; i++, value += step) array[i] = value;
}
if(args[4].u_obj == mp_const_false) {
@ -94,7 +94,7 @@ mp_obj_t numerical_sum_mean_std_array(mp_obj_t oin, uint8_t optype) {
return mp_obj_new_float(sum/len);
} else {
sum /= len; // this is now the mean!
return mp_obj_new_float(sqrtf((sq_sum/len-sum*sum)));
return mp_obj_new_float(sqrt((sq_sum/len-sum*sum)));
}
}
@ -119,7 +119,7 @@ STATIC mp_float_t numerical_sum_mean_std_single_line(void *data, size_t start, s
return sum/len;
} else {
sum /= len; // this is now the mean!
return sqrtf((sq_sum/len-sum*sum));
return sqrt((sq_sum/len-sum*sum));
}
}
@ -168,7 +168,7 @@ size_t numerical_argmin_argmax_array(ndarray_obj_t *in, size_t start,
} 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, float, start, stop, stride, op);
ARG_MIN_LOOP(in, mp_float_t, start, stop, stride, op);
}
return best_idx;
}
@ -532,8 +532,71 @@ mp_obj_t numerical_diff(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
} else if(in->array->typecode == NDARRAY_INT16) {
CALCULATE_DIFF(in, out, int16_t, M, N, in->n, increment);
} else {
CALCULATE_DIFF(in, out, float, M, N, in->n, increment);
CALCULATE_DIFF(in, out, mp_float_t, M, N, in->n, increment);
}
m_del(int8_t, stencil, n);
return MP_OBJ_FROM_PTR(out);
}
/*
void heapsort(ndarray_obj_t *ndarray, size_t n, size_t increment) {
uint8_t type;
if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {
type = 0;
}
if((ndarray->array->typecode == NDARRAY_UINT16) || (ndarray->array->typecode == NDARRAY_INT16)) {
type = 1;
} else {
type = 2;
}
size_t i, j, k, l;
l = (n >> 1) + 1;
k = n;
uint8_t ui8_tmp, *ui8_array = NULL;
uint16_t ui16_tmp, *ui16_array = NULL;
mp_float_t float_tmp, *float_array = NULL;
for(;;) {
if (l > 1) {
if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {
ui8_tmp = ui8_array[--l];
}
} else {
tmp = array[k];
array[k] = array[1];
if (--k == 1) {
array[1] = tmp;
break;
}
}
i = l;
j = 2 * l;
while(j <= k) {
if(j < k && array[j] < array[j+1]) {
j++;
}
}
if(tmp < array[j]) {
array[i] = array[j];
i = j;
j <<= 1;
} else break;
}
array[i] = tmp;
}
mp_obj_t numerical_sort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
static const mp_arg_t allowed_args[] = {
{ MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } },
{ MP_QSTR_n, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 1 } },
{ MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = -1 } },
};
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);
if(!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type)) {
mp_raise_TypeError("sort argument must be an ndarray");
}
return mp_const_none;
}*/

1
code/numerical.h
View file

@ -29,6 +29,7 @@ mp_obj_t numerical_maximum(mp_obj_t , mp_obj_t );
mp_obj_t numerical_cumsum(size_t , const mp_obj_t *, mp_map_t *);
mp_obj_t numerical_flip(size_t , const mp_obj_t *, mp_map_t *);
mp_obj_t numerical_diff(size_t , const mp_obj_t *, mp_map_t *);
mp_obj_t numerical_sort(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 {\

42
code/poly.c
View file

@ -59,13 +59,13 @@ mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {
mp_obj_t p_item, p_iterable;
mp_float_t x, y;
float *outf = (float *)out->array->items;
mp_float_t *outf = (mp_float_t *)out->array->items;
uint8_t plen = mp_obj_get_int(mp_obj_len_maybe(o_p));
float *p = m_new(float, plen);
mp_float_t *p = m_new(mp_float_t, plen);
p_iterable = mp_getiter(o_p, &p_buf);
uint16_t i = 0;
while((p_item = mp_iternext(p_iterable)) != MP_OBJ_STOP_ITERATION) {
p[i] = (float)mp_obj_get_float(p_item);
p[i] = mp_obj_get_float(p_item);
i++;
}
i = 0;
@ -78,7 +78,7 @@ mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {
}
outf[i++] = y;
}
m_del(float, p, plen);
m_del(mp_float_t, p, plen);
return MP_OBJ_FROM_PTR(out);
}
@ -91,7 +91,7 @@ mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {
}
uint16_t lenx, leny;
uint8_t deg;
float *x, *XT, *y, *prod;
mp_float_t *x, *XT, *y, *prod;
if(n_args == 2) { // only the y values are supplied
// TODO: this is actually not enough: the first argument can very well be a matrix,
@ -102,11 +102,11 @@ mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {
mp_raise_ValueError("more degrees of freedom than data points");
}
lenx = leny;
x = m_new(float, lenx); // assume uniformly spaced data points
x = m_new(mp_float_t, lenx); // assume uniformly spaced data points
for(size_t i=0; i < lenx; i++) {
x[i] = i;
}
y = m_new(float, leny);
y = m_new(mp_float_t, leny);
fill_array_iterable(y, args[0]);
} else if(n_args == 3) {
lenx = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));
@ -118,15 +118,15 @@ mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {
if(leny < deg) {
mp_raise_ValueError("more degrees of freedom than data points");
}
x = m_new(float, lenx);
x = m_new(mp_float_t, lenx);
fill_array_iterable(x, args[0]);
y = m_new(float, leny);
y = m_new(mp_float_t, leny);
fill_array_iterable(y, args[1]);
}
// one could probably express X as a function of XT,
// and thereby save RAM, because X is used only in the product
XT = m_new(float, (deg+1)*leny); // XT is a matrix of shape (deg+1, len) (rows, columns)
XT = m_new(mp_float_t, (deg+1)*leny); // XT is a matrix of shape (deg+1, len) (rows, columns)
for(uint8_t i=0; i < leny; i++) { // column index
XT[i+0*lenx] = 1.0; // top row
for(uint8_t j=1; j < deg+1; j++) { // row index
@ -134,8 +134,8 @@ mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {
}
}
prod = m_new(float, (deg+1)*(deg+1)); // the product matrix is of shape (deg+1, deg+1)
float sum;
prod = m_new(mp_float_t, (deg+1)*(deg+1)); // the product matrix is of shape (deg+1, deg+1)
mp_float_t sum;
for(uint16_t i=0; i < deg+1; i++) { // column index
for(uint16_t j=0; j < deg+1; j++) { // row index
sum = 0.0;
@ -152,10 +152,10 @@ mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {
// Although X was a Vandermonde matrix, whose inverse is guaranteed to exist,
// we bail out here, if prod couldn't be inverted: if the values in x are not all
// distinct, prod is singular
m_del(float, XT, (deg+1)*lenx);
m_del(float, x, lenx);
m_del(float, y, lenx);
m_del(float, prod, (deg+1)*(deg+1));
m_del(mp_float_t, XT, (deg+1)*lenx);
m_del(mp_float_t, x, lenx);
m_del(mp_float_t, y, lenx);
m_del(mp_float_t, prod, (deg+1)*(deg+1));
mp_raise_ValueError("could not invert Vandermonde matrix");
}
// at this point, we have the inverse of X^T * X
@ -168,10 +168,10 @@ mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {
x[i] = sum;
}
// XT is no longer needed
m_del(float, XT, (deg+1)*leny);
m_del(mp_float_t, XT, (deg+1)*leny);
ndarray_obj_t *beta = create_new_ndarray(deg+1, 1, NDARRAY_FLOAT);
float *betav = (float *)beta->array->items;
mp_float_t *betav = (mp_float_t *)beta->array->items;
// x[0..(deg+1)] contains now the product X^T * y; we can get rid of y
m_del(float, y, leny);
@ -183,11 +183,11 @@ mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {
}
betav[i] = sum;
}
m_del(float, x, lenx);
m_del(float, prod, (deg+1)*(deg+1));
m_del(mp_float_t, x, lenx);
m_del(mp_float_t, prod, (deg+1)*(deg+1));
for(uint8_t i=0; i < (deg+1)/2; i++) {
// We have to reverse the array, for the leading coefficient comes first.
SWAP(float, betav[i], betav[deg-i]);
SWAP(mp_float_t, betav[i], betav[deg-i]);
}
return MP_OBJ_FROM_PTR(beta);
}

6
code/vectorise.c
View file

@ -30,7 +30,7 @@ mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) {
if(MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {
ndarray_obj_t *source = MP_OBJ_TO_PTR(o_in);
ndarray_obj_t *ndarray = create_new_ndarray(source->m, source->n, NDARRAY_FLOAT);
float *dataout = (float *)ndarray->array->items;
mp_float_t *dataout = (mp_float_t *)ndarray->array->items;
if(source->array->typecode == NDARRAY_UINT8) {
ITERATE_VECTOR(uint8_t, source, dataout);
} else if(source->array->typecode == NDARRAY_INT8) {
@ -40,14 +40,14 @@ mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) {
} else if(source->array->typecode == NDARRAY_INT16) {
ITERATE_VECTOR(int16_t, source, dataout);
} else {
ITERATE_VECTOR(float, source, dataout);
ITERATE_VECTOR(mp_float_t, source, dataout);
}
return MP_OBJ_FROM_PTR(ndarray);
} else if(MP_OBJ_IS_TYPE(o_in, &mp_type_tuple) || MP_OBJ_IS_TYPE(o_in, &mp_type_list) ||
MP_OBJ_IS_TYPE(o_in, &mp_type_range)) { // i.e., the input is a generic iterable
mp_obj_array_t *o = MP_OBJ_TO_PTR(o_in);
ndarray_obj_t *out = create_new_ndarray(1, o->len, NDARRAY_FLOAT);
float *dataout = (float *)out->array->items;
mp_float_t *dataout = (mp_float_t *)out->array->items;
mp_obj_iter_buf_t iter_buf;
mp_obj_t item, iterable = mp_getiter(o_in, &iter_buf);
size_t i=0;

57
docs/manual/source/ulab.rst
View file

@ -121,7 +121,26 @@ denotes those cases, where a bit more caution should be exercised,
though this usually means functions that are not supported by ``numpy``.
The detailed discussion of the various functions always contains a link
to the corresponding ``numpy`` documentation.
to the corresponding ``numpy`` documentation. However, before going down
the rabbit hole, the module also defines a constant, the version, which
can always be queried as
.. code::
# code to be run in micropython
import ulab as np
print('you are running ulab version', np.__version__)
.. parsed-literal::
you are running ulab version 0.24
If you find a bug, please, include this number in your report!
Basic ndarray operations
------------------------
@ -231,7 +250,12 @@ type-aware, i.e., one can save RAM by specifying a data type, and using
the smallest reasonable one. Five such types are defined, namely
``uint8``, ``int8``, which occupy a single byte of memory per datum,
``uint16``, and ``int16``, which occupy two bytes per datum, and
``float``, which occupies four bytes per datum.
``float``, which occupies four or eight bytes per datum. The
precision/size of the ``float`` type depends on the definition of
``mp_float_t``. Some platforms, e.g., the PYBD, implement ``double``\ s,
but some, e.g., the pyboard.v.11, don’t. You can find out, what type of
float your particular platform implements by looking at the output of
the `.rawsize <#.rawsize>`__ class method.
On the following pages, we will see how one can work with
``ndarray``\ s. Those familiar with ``numpy`` should find that the
@ -333,6 +357,22 @@ themselves, though it should be pointed out that initialising through
arrays should be faster, because simply a new copy is created, without
inspection, iteration etc.
.. code::
# code to be run in micropython
import ulab as np
a = np.array([1, 2, 3, 4, 5, 6, 7, 8])
print(a+5)
.. parsed-literal::
here is an intarray([6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0], dtype=float)
.. code::
# code to be run in micropython
@ -441,7 +481,8 @@ following data
3. length of the storage (should be equal to the product of 1. and 2.)
4. length of the data storage in bytes
5. datum size in bytes (1 for ``uint8``/``int8``, 2 for
``uint16``/``int16``, and 4 for ``floats``)
``uint16``/``int16``, and 4, or 8 for ``floats``, see `ndarray, the
basic container <#ndarray,-the-basic-container>`__)
**WARNING:** ``rawsize`` is a ``ulab``-only method; it has no equivalent
in ``numpy``.
@ -2370,10 +2411,11 @@ polyval takes two arguments, both arrays or other iterables.
coefficients: [1, 1, 1, 0]
independent values: [0, 1, 2, 3, 4]
values of p(x): array([0.0, 3.0, 14.0, 39.0, 84.00001], dtype=float)
values of p(x): array([0.0, 3.0, 14.0, 39.0, 84.0], dtype=float)
ndarray (a): array([0.0, 1.0, 2.0, 3.0, 4.0], dtype=float)
value of p(a): array([0.0, 3.0, 14.0, 39.0, 84.00001], dtype=float)
value of p(a): array([0.0, 3.0, 14.0, 39.0, 84.0], dtype=float)
@ -2413,10 +2455,11 @@ a ``ValueError``.
independent values: array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0], dtype=float)
dependent values: array([9.0, 4.0, 1.0, 0.0, 1.0, 4.0, 9.0], dtype=float)
fitted values: array([1.0, -6.000003, 9.000003], dtype=float)
fitted values: array([1.0, -6.0, 9.000000000000004], dtype=float)
dependent values: array([9.0, 4.0, 1.0, 0.0, 1.0, 4.0, 9.0], dtype=float)
fitted values: array([1.0, -6.000003, 9.000003], dtype=float)
fitted values: array([1.0, -6.0, 9.000000000000004], dtype=float)

107
docs/ulab-manual.ipynb
View file

@ -31,11 +31,11 @@
},
{
"cell_type": "code",
"execution_count": 133,
"execution_count": 210,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-29T20:09:58.501920Z",
"start_time": "2019-10-29T20:09:48.007101Z"
"end_time": "2019-11-01T13:26:39.807992Z",
"start_time": "2019-11-01T13:26:33.973550Z"
}
},
"outputs": [],
@ -442,8 +442,48 @@
"\n",
"`ulab` supports a number of array operators, which are listed here. I tried to follow the specifications of the `numpy` interface as closely as possible, though, it was not always practical to implement verbatim behaviour. The differences, if any, are in each case small (e.g., a function cannot take all possible keyword arguments), and should not hinder everyday use. In the list below, a single asterisk denotes slight deviations from `numpy`'s nomenclature, and a double asterisk denotes those cases, where a bit more caution should be exercised, though this usually means functions that are not supported by `numpy`.\n",
"\n",
"The detailed discussion of the various functions always contains a link to the corresponding `numpy` documentation. \n",
"The detailed discussion of the various functions always contains a link to the corresponding `numpy` documentation. However, before going down the rabbit hole, the module also defines a constant, the version, which can always be queried as "
]
},
{
"cell_type": "code",
"execution_count": 209,
"metadata": {
"ExecuteTime": {
"end_time": "2019-11-01T13:19:05.158219Z",
"start_time": "2019-11-01T13:19:05.140455Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"you are running ulab version 0.24\n",
"\n",
"\n"
]
}
],
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"\n",
"print('you are running ulab version', np.__version__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you find a bug, please, include this number in your report!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Basic ndarray operations\n",
"\n",
"[Unary operators](#Unary-operators)\n",
@ -539,7 +579,7 @@
"\n",
"The `ndarray` is the underlying container of numerical data. It is derived from micropython's own `array` object, but has a great number of extra features starting with how it can be initialised, how operations can be done on it, and which functions can accept it as an argument.\n",
"\n",
"Since the `ndarray` is a binary container, it is also compact, meaning that it takes only a couple of bytes of extra RAM in addition to what is required for storing the numbers themselves. `ndarray`s are also type-aware, i.e., one can save RAM by specifying a data type, and using the smallest reasonable one. Five such types are defined, namely `uint8`, `int8`, which occupy a single byte of memory per datum, `uint16`, and `int16`, which occupy two bytes per datum, and `float`, which occupies four bytes per datum. \n",
"Since the `ndarray` is a binary container, it is also compact, meaning that it takes only a couple of bytes of extra RAM in addition to what is required for storing the numbers themselves. `ndarray`s are also type-aware, i.e., one can save RAM by specifying a data type, and using the smallest reasonable one. Five such types are defined, namely `uint8`, `int8`, which occupy a single byte of memory per datum, `uint16`, and `int16`, which occupy two bytes per datum, and `float`, which occupies four or eight bytes per datum. The precision/size of the `float` type depends on the definition of `mp_float_t`. Some platforms, e.g., the PYBD, implement `double`s, but some, e.g., the pyboard.v.11, don't. You can find out, what type of float your particular platform implements by looking at the output of the [.rawsize](#.rawsize) class method.\n",
"\n",
"On the following pages, we will see how one can work with `ndarray`s. Those familiar with `numpy` should find that the nomenclature and naming conventions of `numpy` are adhered to as closely as possible. I will point out the few differences, where necessary.\n",
"\n",
@ -661,6 +701,35 @@
"An `ndarray` can be initialised by supplying another array. This statement is almost trivial, since `ndarray`s are iterables themselves, though it should be pointed out that initialising through arrays should be faster, because simply a new copy is created, without inspection, iteration etc."
]
},
{
"cell_type": "code",
"execution_count": 180,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-31T22:21:28.834732Z",
"start_time": "2019-10-31T22:21:28.821423Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"here is an intarray([6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0], dtype=float)\n",
"\n",
"\n"
]
}
],
"source": [
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"\n",
"a = np.array([1, 2, 3, 4, 5, 6, 7, 8])\n",
"print(a+5)"
]
},
{
"cell_type": "code",
"execution_count": 175,
@ -815,7 +884,7 @@
"2. number of columns\n",
"3. length of the storage (should be equal to the product of 1. and 2.)\n",
"4. length of the data storage in bytes \n",
"5. datum size in bytes (1 for `uint8`/`int8`, 2 for `uint16`/`int16`, and 4 for `floats`)\n",
"5. datum size in bytes (1 for `uint8`/`int8`, 2 for `uint16`/`int16`, and 4, or 8 for `floats`, see [ndarray, the basic container](#ndarray,-the-basic-container))\n",
"\n",
"**WARNING:** `rawsize` is a `ulab`-only method; it has no equivalent in `numpy`."
]
@ -3373,11 +3442,11 @@
},
{
"cell_type": "code",
"execution_count": 497,
"execution_count": 187,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-19T13:27:49.950514Z",
"start_time": "2019-10-19T13:27:49.911147Z"
"end_time": "2019-11-01T12:53:22.448303Z",
"start_time": "2019-11-01T12:53:22.435176Z"
}
},
"outputs": [
@ -3388,16 +3457,17 @@
"coefficients: [1, 1, 1, 0]\n",
"independent values: [0, 1, 2, 3, 4]\n",
"\n",
"values of p(x): array([0.0, 3.0, 14.0, 39.0, 84.00001], dtype=float)\n",
"values of p(x): array([0.0, 3.0, 14.0, 39.0, 84.0], dtype=float)\n",
"\n",
"ndarray (a): array([0.0, 1.0, 2.0, 3.0, 4.0], dtype=float)\n",
"value of p(a): array([0.0, 3.0, 14.0, 39.0, 84.00001], dtype=float)\n",
"value of p(a): array([0.0, 3.0, 14.0, 39.0, 84.0], dtype=float)\n",
"\n",
"\n"
]
}
],
"source": [
"%%micropython -pyboard 1\n",
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"\n",
@ -3428,11 +3498,11 @@
},
{
"cell_type": "code",
"execution_count": 499,
"execution_count": 189,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-19T13:28:54.816888Z",
"start_time": "2019-10-19T13:28:54.781067Z"
"end_time": "2019-11-01T12:54:08.326802Z",
"start_time": "2019-11-01T12:54:08.311182Z"
}
},
"outputs": [
@ -3442,16 +3512,17 @@
"text": [
"independent values:\t array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0], dtype=float)\n",
"dependent values:\t array([9.0, 4.0, 1.0, 0.0, 1.0, 4.0, 9.0], dtype=float)\n",
"fitted values:\t\t array([1.0, -6.000003, 9.000003], dtype=float)\n",
"fitted values:\t\t array([1.0, -6.0, 9.000000000000004], dtype=float)\n",
"\n",
"dependent values:\t array([9.0, 4.0, 1.0, 0.0, 1.0, 4.0, 9.0], dtype=float)\n",
"fitted values:\t\t array([1.0, -6.000003, 9.000003], dtype=float)\n",
"fitted values:\t\t array([1.0, -6.0, 9.000000000000004], dtype=float)\n",
"\n",
"\n"
]
}
],
"source": [
"%%micropython -pyboard 1\n",
"%%micropython -unix 1\n",
"\n",
"import ulab as np\n",
"\n",

399
docs/ulab.ipynb
View file

@ -3141,11 +3141,11 @@
},
{
"cell_type": "code",
"execution_count": 2155,
"execution_count": 373,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-25T05:00:47.029659Z",
"start_time": "2019-10-25T05:00:46.987223Z"
"end_time": "2019-11-01T13:05:40.505260Z",
"start_time": "2019-11-01T13:05:40.499409Z"
}
},
"outputs": [
@ -3153,7 +3153,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 4684 bytes to ndarray.h\n"
"written 4891 bytes to ndarray.h\n"
]
}
],
@ -3170,6 +3170,11 @@
"\n",
"#define PRINT_MAX 10\n",
"\n",
"#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT\n",
"#define FLOAT_TYPECODE 'f'\n",
"#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE\n",
"#define FLOAT_TYPECODE 'd'\n",
"#endif\n",
"\n",
"const mp_obj_type_t ulab_ndarray_type;\n",
"\n",
@ -3178,7 +3183,7 @@
" NDARRAY_INT8 = 'b',\n",
" NDARRAY_UINT16 = 'H', \n",
" NDARRAY_INT16 = 'h',\n",
" NDARRAY_FLOAT = 'f',\n",
" NDARRAY_FLOAT = FLOAT_TYPECODE,\n",
"};\n",
"\n",
"typedef struct _ndarray_obj_t {\n",
@ -3191,8 +3196,8 @@
"\n",
"mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t , size_t , mp_obj_iter_buf_t *);\n",
"\n",
"float ndarray_get_float_value(void *, uint8_t , size_t );\n",
"void fill_array_iterable(float *, mp_obj_t );\n",
"mp_float_t ndarray_get_float_value(void *, uint8_t , size_t );\n",
"void fill_array_iterable(mp_float_t *, mp_obj_t );\n",
"\n",
"void ndarray_print_row(const mp_print_t *, mp_obj_array_t *, size_t , size_t );\n",
"void ndarray_print(const mp_print_t *, mp_obj_t , mp_print_kind_t );\n",
@ -3239,8 +3244,8 @@
" return MP_OBJ_FROM_PTR(out);\\\n",
" } else if((op) == MP_BINARY_OP_TRUE_DIVIDE) {\\\n",
" ndarray_obj_t *out = create_new_ndarray(ol->m, ol->n, NDARRAY_FLOAT);\\\n",
" float *odata = (float *)out->array->items;\\\n",
" for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = (float)left[i]/(float)right[j];\\\n",
" mp_float_t *odata = (mp_float_t *)out->array->items;\\\n",
" for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = (mp_float_t)left[i]/(mp_float_t)right[j];\\\n",
" return MP_OBJ_FROM_PTR(out);\\\n",
" } else if(((op) == MP_BINARY_OP_LESS) || ((op) == MP_BINARY_OP_LESS_EQUAL) || \\\n",
" ((op) == MP_BINARY_OP_MORE) || ((op) == MP_BINARY_OP_MORE_EQUAL)) {\\\n",
@ -3268,15 +3273,6 @@
" }\\\n",
"} while(0)\n",
"\n",
"#define ASSIGN_SLICE(ndarray, ndarray_type, value, value_type) do {\\\n",
" ndarray_type *target = (ndarray_type *)(ndarray)->\n",
" ndarray_obj_t *tmp = create_new_ndarray(1, 1, (typecode));\\\n",
" type *tmparr = (type *)tmp->array->items;\\\n",
" tmparr[0] = (type)(value);\\\n",
" (outarray) = MP_OBJ_FROM_PTR(tmp);\\\n",
"} while(0)\n",
"\n",
"\n",
"#endif"
]
},
@ -3289,11 +3285,11 @@
},
{
"cell_type": "code",
"execution_count": 172,
"execution_count": 374,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-29T20:09:06.794242Z",
"start_time": "2019-10-29T20:09:06.778107Z"
"end_time": "2019-11-01T13:08:29.333165Z",
"start_time": "2019-11-01T13:08:29.322821Z"
},
"code_folding": []
},
@ -3302,7 +3298,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 40616 bytes to ndarray.c\n"
"written 41196 bytes to ndarray.c\n"
]
}
],
@ -3338,26 +3334,26 @@
" return o;\n",
"}\n",
"\n",
"float ndarray_get_float_value(void *data, uint8_t typecode, size_t index) {\n",
"mp_float_t ndarray_get_float_value(void *data, uint8_t typecode, size_t index) {\n",
" if(typecode == NDARRAY_UINT8) {\n",
" return (float)((uint8_t *)data)[index];\n",
" return (mp_float_t)((uint8_t *)data)[index];\n",
" } else if(typecode == NDARRAY_INT8) {\n",
" return (float)((int8_t *)data)[index];\n",
" return (mp_float_t)((int8_t *)data)[index];\n",
" } else if(typecode == NDARRAY_UINT16) {\n",
" return (float)((uint16_t *)data)[index];\n",
" return (mp_float_t)((uint16_t *)data)[index];\n",
" } else if(typecode == NDARRAY_INT16) {\n",
" return (float)((int16_t *)data)[index];\n",
" return (mp_float_t)((int16_t *)data)[index];\n",
" } else {\n",
" return (float)((float_t *)data)[index];\n",
" return (mp_float_t)((float_t *)data)[index];\n",
" }\n",
"}\n",
"\n",
"void fill_array_iterable(float *array, mp_obj_t iterable) {\n",
"void fill_array_iterable(mp_float_t *array, mp_obj_t iterable) {\n",
" mp_obj_iter_buf_t x_buf;\n",
" mp_obj_t x_item, x_iterable = mp_getiter(iterable, &x_buf);\n",
" size_t i=0;\n",
" while ((x_item = mp_iternext(x_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" array[i] = (float)mp_obj_get_float(x_item);\n",
" array[i] = (mp_float_t)mp_obj_get_float(x_item);\n",
" i++;\n",
" }\n",
"}\n",
@ -3418,7 +3414,7 @@
" mp_print_str(print, \", dtype=int16)\");\n",
" } else if(self->array->typecode == NDARRAY_FLOAT) {\n",
" mp_print_str(print, \", dtype=float)\");\n",
" } \n",
" }\n",
"}\n",
"\n",
"void ndarray_assign_elements(mp_obj_array_t *data, mp_obj_t iterable, uint8_t typecode, size_t *idx) {\n",
@ -3987,6 +3983,9 @@
"// Binary operations\n",
"\n",
"mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) {\n",
"// if(op == MP_BINARY_OP_REVERSE_ADD) {\n",
" // return ndarray_binary_op(MP_BINARY_OP_ADD, rhs, lhs);\n",
" // } \n",
" // One of the operands is a scalar\n",
" // TODO: conform to numpy with the upcasting\n",
" // TODO: implement in-place operators\n",
@ -3996,19 +3995,18 @@
" int32_t ivalue = mp_obj_get_int(rhs);\n",
" if((ivalue > 0) && (ivalue < 256)) {\n",
" CREATE_SINGLE_ITEM(RHS, uint8_t, NDARRAY_UINT8, ivalue);\n",
" }\n",
" else if((ivalue > 0) && (ivalue < 65535)) {\n",
" } else if((ivalue > 255) && (ivalue < 65535)) {\n",
" CREATE_SINGLE_ITEM(RHS, uint16_t, NDARRAY_UINT16, ivalue);\n",
" }\n",
" else if((ivalue < 0) && (ivalue > -128)) {\n",
" } else if((ivalue < 0) && (ivalue > -128)) {\n",
" CREATE_SINGLE_ITEM(RHS, int8_t, NDARRAY_INT8, ivalue);\n",
" }\n",
" else if((ivalue < 0) && (ivalue > -32767)) {\n",
" } else if((ivalue < -127) && (ivalue > -32767)) {\n",
" CREATE_SINGLE_ITEM(RHS, int16_t, NDARRAY_INT16, ivalue);\n",
" } else { // the integer value clearly does not fit the ulab types, so move on to float\n",
" CREATE_SINGLE_ITEM(RHS, mp_float_t, NDARRAY_FLOAT, ivalue);\n",
" }\n",
" } else if(mp_obj_is_float(rhs)) {\n",
" float fvalue = mp_obj_get_float(rhs); \n",
" CREATE_SINGLE_ITEM(RHS, float, NDARRAY_FLOAT, fvalue);\n",
" mp_float_t fvalue = mp_obj_get_float(rhs); \n",
" CREATE_SINGLE_ITEM(RHS, mp_float_t, NDARRAY_FLOAT, fvalue);\n",
" } else {\n",
" RHS = rhs;\n",
" rhs_is_scalar = false;\n",
@ -4071,7 +4069,7 @@
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint8_t, float, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ol, or, op);\n",
" }\n",
" } else if(ol->array->typecode == NDARRAY_INT8) {\n",
" if(or->array->typecode == NDARRAY_UINT8) {\n",
@ -4083,7 +4081,7 @@
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, int8_t, float, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int8_t, mp_float_t, ol, or, op);\n",
" } \n",
" } else if(ol->array->typecode == NDARRAY_UINT16) {\n",
" if(or->array->typecode == NDARRAY_UINT8) {\n",
@ -4093,9 +4091,9 @@
" } else if(or->array->typecode == NDARRAY_UINT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint16_t, int16_t, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint8_t, float, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ol, or, op);\n",
" }\n",
" } else if(ol->array->typecode == NDARRAY_INT16) {\n",
" if(or->array->typecode == NDARRAY_UINT8) {\n",
@ -4103,23 +4101,23 @@
" } else if(or->array->typecode == NDARRAY_INT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_UINT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, int16_t, uint16_t, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int16_t, uint16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint16_t, float, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, mp_float_t, ol, or, op);\n",
" }\n",
" } else if(ol->array->typecode == NDARRAY_FLOAT) {\n",
" if(or->array->typecode == NDARRAY_UINT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, uint8_t, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, int8_t, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_UINT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, uint16_t, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, int16_t, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, float, ol, or, op);\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, mp_float_t, ol, or, op);\n",
" }\n",
" } else { // this should never happen\n",
" mp_raise_TypeError(\"wrong input type\");\n",
@ -4173,7 +4171,7 @@
" int16_t *array = (int16_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) array[i] = -array[i];\n",
" } else {\n",
" float *array = (float *)ndarray->array->items;\n",
" mp_float_t *array = (mp_float_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) array[i] = -array[i];\n",
" }\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
@ -4198,7 +4196,7 @@
" if(array[i] < 0) array[i] = -array[i];\n",
" }\n",
" } else {\n",
" float *array = (float *)ndarray->array->items;\n",
" mp_float_t *array = (mp_float_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) {\n",
" if(array[i] < 0) array[i] = -array[i];\n",
" } \n",
@ -5059,11 +5057,11 @@
},
{
"cell_type": "code",
"execution_count": 1990,
"execution_count": 362,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-19T09:44:54.199545Z",
"start_time": "2019-10-19T09:44:54.192668Z"
"end_time": "2019-11-01T12:52:24.062644Z",
"start_time": "2019-11-01T12:52:24.057127Z"
}
},
"outputs": [
@ -5071,7 +5069,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 867 bytes to linalg.h\n"
"written 1019 bytes to linalg.h\n"
]
}
],
@ -5084,13 +5082,19 @@
"#include \"ndarray.h\"\n",
"\n",
"#define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }\n",
"#define epsilon 1e-6\n",
"\n",
"#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT\n",
"#define epsilon 1.2e-7\n",
"#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE\n",
"#define epsilon 2.3e-16\n",
"#endif\n",
"\n",
"#define JACOBI_MAX 20\n",
"\n",
"mp_obj_t linalg_transpose(mp_obj_t );\n",
"mp_obj_t linalg_reshape(mp_obj_t , mp_obj_t );\n",
"mp_obj_t linalg_size(size_t , const mp_obj_t *, mp_map_t *);\n",
"bool linalg_invert_matrix(float *, size_t );\n",
"bool linalg_invert_matrix(mp_float_t *, size_t );\n",
"mp_obj_t linalg_inv(mp_obj_t );\n",
"mp_obj_t linalg_dot(mp_obj_t , mp_obj_t );\n",
"mp_obj_t linalg_zeros(size_t , const mp_obj_t *, mp_map_t *);\n",
@ -5112,11 +5116,11 @@
},
{
"cell_type": "code",
"execution_count": 1997,
"execution_count": 358,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-19T10:19:04.522356Z",
"start_time": "2019-10-19T10:19:04.512259Z"
"end_time": "2019-11-01T12:33:00.166042Z",
"start_time": "2019-11-01T12:33:00.158937Z"
},
"code_folding": []
},
@ -5125,7 +5129,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 17287 bytes to linalg.c\n"
"written 17438 bytes to linalg.c\n"
]
}
],
@ -5232,24 +5236,24 @@
" }\n",
"}\n",
"\n",
"bool linalg_invert_matrix(float *data, size_t N) {\n",
"bool linalg_invert_matrix(mp_float_t *data, size_t N) {\n",
" // returns true, of the inversion was successful, \n",
" // false, if the matrix is singular\n",
" \n",
" // initially, this is the unit matrix: the contents of this matrix is what \n",
" // will be returned after all the transformations\n",
" float *unit = m_new(float, N*N);\n",
" mp_float_t *unit = m_new(mp_float_t, N*N);\n",
"\n",
" float elem = 1.0;\n",
" mp_float_t elem = 1.0;\n",
" // initialise the unit matrix\n",
" memset(unit, 0, sizeof(float)*N*N);\n",
" memset(unit, 0, sizeof(mp_float_t)*N*N);\n",
" for(size_t m=0; m < N; m++) {\n",
" memcpy(&unit[m*(N+1)], &elem, sizeof(float));\n",
" memcpy(&unit[m*(N+1)], &elem, sizeof(mp_float_t));\n",
" }\n",
" for(size_t m=0; m < N; m++){\n",
" // this could be faster with ((c < epsilon) && (c > -epsilon))\n",
" if(abs(data[m*(N+1)]) < epsilon) {\n",
" m_del(float, unit, N*N);\n",
" m_del(mp_float_t, unit, N*N);\n",
" return false;\n",
" }\n",
" for(size_t n=0; n < N; n++){\n",
@ -5269,8 +5273,8 @@
" unit[N*m+n] /= elem;\n",
" }\n",
" }\n",
" memcpy(data, unit, sizeof(float)*N*N);\n",
" m_del(float, unit, N*N);\n",
" memcpy(data, unit, sizeof(mp_float_t)*N*N);\n",
" m_del(mp_float_t, unit, N*N);\n",
" return true;\n",
"}\n",
"\n",
@ -5287,21 +5291,21 @@
" mp_raise_ValueError(\"only square matrices can be inverted\");\n",
" }\n",
" ndarray_obj_t *inverted = create_new_ndarray(o->m, o->n, NDARRAY_FLOAT);\n",
" float *data = (float *)inverted->array->items;\n",
" mp_float_t *data = (mp_float_t *)inverted->array->items;\n",
" mp_obj_t elem;\n",
" for(size_t m=0; m < o->m; m++) { // rows first\n",
" for(size_t n=0; n < o->n; n++) { // columns next\n",
" // this could, perhaps, be done in single line... \n",
" // On the other hand, we probably spend little time here\n",
" elem = mp_binary_get_val_array(o->array->typecode, o->array->items, m*o->n+n);\n",
" data[m*o->n+n] = (float)mp_obj_get_float(elem);\n",
" data[m*o->n+n] = (mp_float_t)mp_obj_get_float(elem);\n",
" }\n",
" }\n",
" \n",
" if(!linalg_invert_matrix(data, o->m)) {\n",
" // TODO: I am not sure this is needed here. Otherwise, \n",
" // how should we free up the unused RAM of inverted?\n",
" m_del(float, inverted->array->items, o->n*o->n);\n",
" m_del(mp_float_t, inverted->array->items, o->n*o->n);\n",
" mp_raise_ValueError(\"input matrix is singular\");\n",
" }\n",
" return MP_OBJ_FROM_PTR(inverted);\n",
@ -5316,8 +5320,8 @@
" }\n",
" // TODO: numpy uses upcasting here\n",
" ndarray_obj_t *out = create_new_ndarray(m1->m, m2->n, NDARRAY_FLOAT);\n",
" float *outdata = (float *)out->array->items;\n",
" float sum, v1, v2;\n",
" mp_float_t *outdata = (mp_float_t *)out->array->items;\n",
" mp_float_t sum, v1, v2;\n",
" for(size_t i=0; i < m1->n; i++) {\n",
" for(size_t j=0; j < m2->m; j++) {\n",
" sum = 0.0;\n",
@ -5425,14 +5429,14 @@
" mp_raise_ValueError(\"input must be square matrix\");\n",
" }\n",
" \n",
" float *tmp = m_new(float, in->n*in->n);\n",
" mp_float_t *tmp = m_new(mp_float_t, in->n*in->n);\n",
" for(size_t i=0; i < in->array->len; i++){\n",
" tmp[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);\n",
" }\n",
" float c;\n",
" mp_float_t c;\n",
" for(size_t m=0; m < in->m-1; m++){\n",
" if(abs(tmp[m*(in->n+1)]) < epsilon) {\n",
" m_del(float, tmp, in->n*in->n);\n",
" m_del(mp_float_t, tmp, in->n*in->n);\n",
" return mp_obj_new_float(0.0);\n",
" }\n",
" for(size_t n=0; n < in->n; n++){\n",
@ -5444,12 +5448,12 @@
" }\n",
" }\n",
" }\n",
" float det = 1.0;\n",
" mp_float_t det = 1.0;\n",
" \n",
" for(size_t m=0; m < in->m; m++){ \n",
" det *= tmp[m*(in->n+1)];\n",
" }\n",
" m_del(float, tmp, in->n*in->n);\n",
" m_del(mp_float_t, tmp, in->n*in->n);\n",
" return mp_obj_new_float(det);\n",
"}\n",
"\n",
@ -5461,7 +5465,7 @@
" if(in->m != in->n) {\n",
" mp_raise_ValueError(\"input must be square matrix\");\n",
" }\n",
" float *array = m_new(float, in->array->len);\n",
" mp_float_t *array = m_new(mp_float_t, in->array->len);\n",
" for(size_t i=0; i < in->array->len; i++) {\n",
" array[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);\n",
" }\n",
@ -5478,12 +5482,12 @@
" // if we got this far, then the matrix will be symmetric\n",
" \n",
" ndarray_obj_t *eigenvectors = create_new_ndarray(in->m, in->n, NDARRAY_FLOAT);\n",
" float *eigvectors = (float *)eigenvectors->array->items;\n",
" mp_float_t *eigvectors = (mp_float_t *)eigenvectors->array->items;\n",
" // start out with the unit matrix\n",
" for(size_t m=0; m < in->m; m++) {\n",
" eigvectors[m*(in->n+1)] = 1.0;\n",
" }\n",
" float largest, w, t, c, s, tau, aMk, aNk, vm, vn;\n",
" mp_float_t largest, w, t, c, s, tau, aMk, aNk, vm, vn;\n",
" size_t M, N;\n",
" size_t iterations = JACOBI_MAX*in->n*in->n;\n",
" do {\n",
@ -5511,12 +5515,12 @@
" // The following if/else chooses the smaller absolute value for the tangent \n",
" // of the rotation angle. Going with the smaller should be numerically stabler.\n",
" if(w > 0) {\n",
" t = sqrtf(w*w + 1.0) - w;\n",
" t = sqrt(w*w + 1.0) - w;\n",
" } else {\n",
" t = -1.0*(sqrtf(w*w + 1.0) + w);\n",
" t = -1.0*(sqrt(w*w + 1.0) + w);\n",
" }\n",
" s = t / sqrtf(t*t + 1.0); // the sine of the rotation angle\n",
" c = 1.0 / sqrtf(t*t + 1.0); // the cosine of the rotation angle\n",
" s = t / sqrt(t*t + 1.0); // the sine of the rotation angle\n",
" c = 1.0 / sqrt(t*t + 1.0); // the cosine of the rotation angle\n",
" tau = (1.0-c)/s; // this is equal to the tangent of the half of the rotation angle\n",
" \n",
" // at this point, we have the rotation angles, so we can transform the matrix\n",
@ -5562,15 +5566,15 @@
" \n",
" if(iterations == 0) { \n",
" // the computation did not converge; numpy raises LinAlgError\n",
" m_del(float, array, in->array->len);\n",
" m_del(mp_float_t, array, in->array->len);\n",
" mp_raise_ValueError(\"iterations did not converge\");\n",
" }\n",
" ndarray_obj_t *eigenvalues = create_new_ndarray(1, in->n, NDARRAY_FLOAT);\n",
" float *eigvalues = (float *)eigenvalues->array->items;\n",
" mp_float_t *eigvalues = (mp_float_t *)eigenvalues->array->items;\n",
" for(size_t i=0; i < in->n; i++) {\n",
" eigvalues[i] = array[i*(in->n+1)];\n",
" }\n",
" m_del(float, array, in->array->len);\n",
" m_del(mp_float_t, array, in->array->len);\n",
" \n",
" mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(mp_obj_new_tuple(2, NULL));\n",
" tuple->items[0] = MP_OBJ_FROM_PTR(eigenvalues);\n",
@ -5767,11 +5771,11 @@
},
{
"cell_type": "code",
"execution_count": 1666,
"execution_count": 370,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-17T05:28:21.378139Z",
"start_time": "2019-10-17T05:28:21.371481Z"
"end_time": "2019-11-01T13:03:15.833279Z",
"start_time": "2019-11-01T13:03:15.826010Z"
}
},
"outputs": [
@ -5779,7 +5783,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 3233 bytes to vectorise.c\n"
"written 3258 bytes to vectorise.c\n"
]
}
],
@ -5808,7 +5812,7 @@
" if(MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {\n",
" ndarray_obj_t *source = MP_OBJ_TO_PTR(o_in);\n",
" ndarray_obj_t *ndarray = create_new_ndarray(source->m, source->n, NDARRAY_FLOAT);\n",
" float *dataout = (float *)ndarray->array->items;\n",
" mp_float_t *dataout = (mp_float_t *)ndarray->array->items;\n",
" if(source->array->typecode == NDARRAY_UINT8) {\n",
" ITERATE_VECTOR(uint8_t, source, dataout);\n",
" } else if(source->array->typecode == NDARRAY_INT8) {\n",
@ -5818,14 +5822,14 @@
" } else if(source->array->typecode == NDARRAY_INT16) {\n",
" ITERATE_VECTOR(int16_t, source, dataout);\n",
" } else {\n",
" ITERATE_VECTOR(float, source, dataout);\n",
" ITERATE_VECTOR(mp_float_t, source, dataout);\n",
" }\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
" } else if(MP_OBJ_IS_TYPE(o_in, &mp_type_tuple) || MP_OBJ_IS_TYPE(o_in, &mp_type_list) || \n",
" MP_OBJ_IS_TYPE(o_in, &mp_type_range)) { // i.e., the input is a generic iterable\n",
" mp_obj_array_t *o = MP_OBJ_TO_PTR(o_in);\n",
" ndarray_obj_t *out = create_new_ndarray(1, o->len, NDARRAY_FLOAT);\n",
" float *dataout = (float *)out->array->items;\n",
" mp_float_t *dataout = (mp_float_t *)out->array->items;\n",
" mp_obj_iter_buf_t iter_buf;\n",
" mp_obj_t item, iterable = mp_getiter(o_in, &iter_buf);\n",
" size_t i=0;\n",
@ -6117,9 +6121,7 @@
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"metadata": {},
"source": [
"## poly.h"
]
@ -6131,8 +6133,7 @@
"ExecuteTime": {
"end_time": "2019-09-21T05:29:23.021150Z",
"start_time": "2019-09-21T05:29:23.014023Z"
},
"hidden": true
}
},
"outputs": [
{
@ -6157,29 +6158,26 @@
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"metadata": {},
"source": [
"## poly.c"
]
},
{
"cell_type": "code",
"execution_count": 1328,
"execution_count": 346,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-12T20:24:49.994676Z",
"start_time": "2019-10-12T20:24:47.825965Z"
},
"hidden": true
"end_time": "2019-11-01T12:09:15.148115Z",
"start_time": "2019-11-01T12:09:15.073086Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 6751 bytes to poly.c\n"
"written 6859 bytes to poly.c\n"
]
}
],
@ -6237,13 +6235,13 @@
" mp_obj_t p_item, p_iterable;\n",
"\n",
" mp_float_t x, y;\n",
" float *outf = (float *)out->array->items;\n",
" mp_float_t *outf = (mp_float_t *)out->array->items;\n",
" uint8_t plen = mp_obj_get_int(mp_obj_len_maybe(o_p));\n",
" float *p = m_new(float, plen);\n",
" mp_float_t *p = m_new(mp_float_t, plen);\n",
" p_iterable = mp_getiter(o_p, &p_buf);\n",
" uint16_t i = 0; \n",
" while((p_item = mp_iternext(p_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" p[i] = (float)mp_obj_get_float(p_item);\n",
" p[i] = mp_obj_get_float(p_item);\n",
" i++;\n",
" }\n",
" i = 0;\n",
@ -6256,7 +6254,7 @@
" }\n",
" outf[i++] = y;\n",
" }\n",
" m_del(float, p, plen);\n",
" m_del(mp_float_t, p, plen);\n",
" return MP_OBJ_FROM_PTR(out);\n",
"}\n",
"\n",
@ -6269,7 +6267,7 @@
" }\n",
" uint16_t lenx, leny;\n",
" uint8_t deg;\n",
" float *x, *XT, *y, *prod;\n",
" mp_float_t *x, *XT, *y, *prod;\n",
"\n",
" if(n_args == 2) { // only the y values are supplied\n",
" // TODO: this is actually not enough: the first argument can very well be a matrix, \n",
@ -6280,11 +6278,11 @@
" mp_raise_ValueError(\"more degrees of freedom than data points\");\n",
" }\n",
" lenx = leny;\n",
" x = m_new(float, lenx); // assume uniformly spaced data points\n",
" x = m_new(mp_float_t, lenx); // assume uniformly spaced data points\n",
" for(size_t i=0; i < lenx; i++) {\n",
" x[i] = i;\n",
" }\n",
" y = m_new(float, leny);\n",
" y = m_new(mp_float_t, leny);\n",
" fill_array_iterable(y, args[0]);\n",
" } else if(n_args == 3) {\n",
" lenx = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));\n",
@ -6296,15 +6294,15 @@
" if(leny < deg) {\n",
" mp_raise_ValueError(\"more degrees of freedom than data points\");\n",
" }\n",
" x = m_new(float, lenx);\n",
" x = m_new(mp_float_t, lenx);\n",
" fill_array_iterable(x, args[0]);\n",
" y = m_new(float, leny);\n",
" y = m_new(mp_float_t, leny);\n",
" fill_array_iterable(y, args[1]);\n",
" }\n",
" \n",
" // one could probably express X as a function of XT, \n",
" // and thereby save RAM, because X is used only in the product\n",
" XT = m_new(float, (deg+1)*leny); // XT is a matrix of shape (deg+1, len) (rows, columns)\n",
" XT = m_new(mp_float_t, (deg+1)*leny); // XT is a matrix of shape (deg+1, len) (rows, columns)\n",
" for(uint8_t i=0; i < leny; i++) { // column index\n",
" XT[i+0*lenx] = 1.0; // top row\n",
" for(uint8_t j=1; j < deg+1; j++) { // row index\n",
@ -6312,8 +6310,8 @@
" }\n",
" }\n",
" \n",
" prod = m_new(float, (deg+1)*(deg+1)); // the product matrix is of shape (deg+1, deg+1)\n",
" float sum;\n",
" prod = m_new(mp_float_t, (deg+1)*(deg+1)); // the product matrix is of shape (deg+1, deg+1)\n",
" mp_float_t sum;\n",
" for(uint16_t i=0; i < deg+1; i++) { // column index\n",
" for(uint16_t j=0; j < deg+1; j++) { // row index\n",
" sum = 0.0;\n",
@ -6330,10 +6328,10 @@
" // Although X was a Vandermonde matrix, whose inverse is guaranteed to exist, \n",
" // we bail out here, if prod couldn't be inverted: if the values in x are not all \n",
" // distinct, prod is singular\n",
" m_del(float, XT, (deg+1)*lenx);\n",
" m_del(float, x, lenx);\n",
" m_del(float, y, lenx);\n",
" m_del(float, prod, (deg+1)*(deg+1));\n",
" m_del(mp_float_t, XT, (deg+1)*lenx);\n",
" m_del(mp_float_t, x, lenx);\n",
" m_del(mp_float_t, y, lenx);\n",
" m_del(mp_float_t, prod, (deg+1)*(deg+1));\n",
" mp_raise_ValueError(\"could not invert Vandermonde matrix\");\n",
" } \n",
" // at this point, we have the inverse of X^T * X\n",
@ -6346,10 +6344,10 @@
" x[i] = sum;\n",
" }\n",
" // XT is no longer needed\n",
" m_del(float, XT, (deg+1)*leny);\n",
" m_del(mp_float_t, XT, (deg+1)*leny);\n",
" \n",
" ndarray_obj_t *beta = create_new_ndarray(deg+1, 1, NDARRAY_FLOAT);\n",
" float *betav = (float *)beta->array->items;\n",
" mp_float_t *betav = (mp_float_t *)beta->array->items;\n",
" // x[0..(deg+1)] contains now the product X^T * y; we can get rid of y\n",
" m_del(float, y, leny);\n",
" \n",
@ -6361,11 +6359,11 @@
" }\n",
" betav[i] = sum;\n",
" }\n",
" m_del(float, x, lenx);\n",
" m_del(float, prod, (deg+1)*(deg+1));\n",
" m_del(mp_float_t, x, lenx);\n",
" m_del(mp_float_t, prod, (deg+1)*(deg+1));\n",
" for(uint8_t i=0; i < (deg+1)/2; i++) {\n",
" // We have to reverse the array, for the leading coefficient comes first. \n",
" SWAP(float, betav[i], betav[deg-i]);\n",
" SWAP(mp_float_t, betav[i], betav[deg-i]);\n",
" }\n",
" return MP_OBJ_FROM_PTR(beta);\n",
"}"
@ -6610,11 +6608,11 @@
},
{
"cell_type": "code",
"execution_count": 1994,
"execution_count": 380,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-19T09:56:10.777011Z",
"start_time": "2019-10-19T09:56:10.770009Z"
"end_time": "2019-11-01T13:12:49.126204Z",
"start_time": "2019-11-01T13:12:49.121112Z"
}
},
"outputs": [
@ -6622,7 +6620,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 5295 bytes to fft.c\n"
"written 5352 bytes to fft.c\n"
]
}
],
@ -6647,13 +6645,13 @@
" FFT_SPECTRUM,\n",
"};\n",
"\n",
"void fft_kernel(float *real, float *imag, int n, int isign) {\n",
"void fft_kernel(mp_float_t *real, mp_float_t *imag, int n, int isign) {\n",
" // This is basically a modification of four1 from Numerical Recipes\n",
" // The main difference is that this function takes two arrays, one \n",
" // for the real, and one for the imaginary parts. \n",
" int j, m, mmax, istep;\n",
" float tempr, tempi;\n",
" float wtemp, wr, wpr, wpi, wi, theta;\n",
" mp_float_t tempr, tempi;\n",
" mp_float_t wtemp, wr, wpr, wpi, wi, theta;\n",
"\n",
" j = 0;\n",
" for(int i = 0; i < n; i++) {\n",
@ -6673,9 +6671,9 @@
" while (n > mmax) {\n",
" istep = mmax << 1;\n",
" theta = -1.0*isign*6.28318530717959/istep;\n",
" wtemp = sinf(0.5 * theta);\n",
" wtemp = sin(0.5 * theta);\n",
" wpr = -2.0 * wtemp * wtemp;\n",
" wpi = sinf(theta);\n",
" wpi = sin(theta);\n",
" wr = 1.0;\n",
" wi = 0.0;\n",
" for(m = 0; m < mmax; m++) {\n",
@ -6713,19 +6711,19 @@
" }\n",
" \n",
" ndarray_obj_t *out_re = create_new_ndarray(1, len, NDARRAY_FLOAT);\n",
" float *data_re = (float *)out_re->array->items;\n",
" mp_float_t *data_re = (mp_float_t *)out_re->array->items;\n",
" \n",
" if(re->array->typecode == NDARRAY_FLOAT) { \n",
" // By treating this case separately, we can save a bit of time.\n",
" // I don't know if it is worthwhile, though...\n",
" memcpy((float *)out_re->array->items, (float *)re->array->items, re->bytes);\n",
" memcpy((mp_float_t *)out_re->array->items, (mp_float_t *)re->array->items, re->bytes);\n",
" } else {\n",
" for(size_t i=0; i < len; i++) {\n",
" data_re[i] = ndarray_get_float_value(re->array->items, re->array->typecode, i);\n",
" }\n",
" }\n",
" ndarray_obj_t *out_im = create_new_ndarray(1, len, NDARRAY_FLOAT);\n",
" float *data_im = (float *)out_im->array->items;\n",
" mp_float_t *data_im = (mp_float_t *)out_im->array->items;\n",
"\n",
" if(n_args == 2) {\n",
" ndarray_obj_t *im = MP_OBJ_TO_PTR(arg_im);\n",
@ -6733,7 +6731,7 @@
" mp_raise_ValueError(\"real and imaginary parts must be of equal length\");\n",
" }\n",
" if(im->array->typecode == NDARRAY_FLOAT) {\n",
" memcpy((float *)out_im->array->items, (float *)im->array->items, im->bytes);\n",
" memcpy((mp_float_t *)out_im->array->items, (mp_float_t *)im->array->items, im->bytes);\n",
" } else {\n",
" for(size_t i=0; i < len; i++) {\n",
" data_im[i] = ndarray_get_float_value(im->array->items, im->array->typecode, i);\n",
@ -6744,7 +6742,7 @@
" fft_kernel(data_re, data_im, len, 1);\n",
" if(type == FFT_SPECTRUM) {\n",
" for(size_t i=0; i < len; i++) {\n",
" data_re[i] = sqrtf(data_re[i]*data_re[i] + data_im[i]*data_im[i]);\n",
" data_re[i] = sqrt(data_re[i]*data_re[i] + data_im[i]*data_im[i]);\n",
" }\n",
" }\n",
" } else { // inverse transform\n",
@ -6986,11 +6984,11 @@
},
{
"cell_type": "code",
"execution_count": 290,
"execution_count": 305,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-31T11:48:57.625930Z",
"start_time": "2019-10-31T11:48:57.619375Z"
"end_time": "2019-10-31T13:13:29.838649Z",
"start_time": "2019-10-31T13:13:29.833361Z"
}
},
"outputs": [
@ -6998,7 +6996,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 2215 bytes to numerical.h\n"
"written 2279 bytes to numerical.h\n"
]
}
],
@ -7026,6 +7024,7 @@
"mp_obj_t numerical_cumsum(size_t , const mp_obj_t *, mp_map_t *);\n",
"mp_obj_t numerical_flip(size_t , const mp_obj_t *, mp_map_t *);\n",
"mp_obj_t numerical_diff(size_t , const mp_obj_t *, mp_map_t *);\n",
"mp_obj_t numerical_sort(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",
@ -7072,11 +7071,11 @@
},
{
"cell_type": "code",
"execution_count": 291,
"execution_count": 378,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-31T11:49:41.515945Z",
"start_time": "2019-10-31T11:49:41.507030Z"
"end_time": "2019-11-01T13:12:18.077600Z",
"start_time": "2019-11-01T13:12:18.066200Z"
}
},
"outputs": [
@ -7084,7 +7083,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"written 22886 bytes to numerical.c\n"
"written 24860 bytes to numerical.c\n"
]
}
],
@ -7127,7 +7126,7 @@
" if(len < 2) {\n",
" mp_raise_ValueError(\"number of points must be at least 2\");\n",
" }\n",
" float value, step;\n",
" mp_float_t value, step;\n",
" value = mp_obj_get_float(args[0].u_obj);\n",
" uint8_t typecode = args[5].u_int;\n",
" if(args[3].u_obj == mp_const_true) step = (mp_obj_get_float(args[1].u_obj)-value)/(len-1);\n",
@ -7146,7 +7145,7 @@
" int16_t *array = (int16_t *)ndarray->array->items;\n",
" for(size_t i=0; i < len; i++, value += step) array[i] = (int16_t)value;\n",
" } else {\n",
" float *array = (float *)ndarray->array->items;\n",
" mp_float_t *array = (mp_float_t *)ndarray->array->items;\n",
" for(size_t i=0; i < len; i++, value += step) array[i] = value;\n",
" }\n",
" if(args[4].u_obj == mp_const_false) {\n",
@ -7177,7 +7176,7 @@
" return mp_obj_new_float(sum/len);\n",
" } else {\n",
" sum /= len; // this is now the mean!\n",
" return mp_obj_new_float(sqrtf((sq_sum/len-sum*sum)));\n",
" return mp_obj_new_float(sqrt((sq_sum/len-sum*sum)));\n",
" }\n",
"}\n",
"\n",
@ -7202,7 +7201,7 @@
" return sum/len;\n",
" } else {\n",
" sum /= len; // this is now the mean!\n",
" return sqrtf((sq_sum/len-sum*sum));\n",
" return sqrt((sq_sum/len-sum*sum));\n",
" }\n",
"}\n",
"\n",
@ -7251,7 +7250,7 @@
" } 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, float, start, stop, stride, op);\n",
" ARG_MIN_LOOP(in, mp_float_t, start, stop, stride, op);\n",
" }\n",
" return best_idx;\n",
"}\n",
@ -7615,11 +7614,74 @@
" } else if(in->array->typecode == NDARRAY_INT16) {\n",
" CALCULATE_DIFF(in, out, int16_t, M, N, in->n, increment);\n",
" } else {\n",
" CALCULATE_DIFF(in, out, float, M, N, in->n, increment);\n",
" CALCULATE_DIFF(in, out, mp_float_t, M, N, in->n, increment);\n",
" }\n",
" m_del(int8_t, stencil, n);\n",
" return MP_OBJ_FROM_PTR(out);\n",
"}"
"}\n",
"/*\n",
"void heapsort(ndarray_obj_t *ndarray, size_t n, size_t increment) {\n",
" uint8_t type;\n",
" if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {\n",
" type = 0;\n",
" }\n",
" if((ndarray->array->typecode == NDARRAY_UINT16) || (ndarray->array->typecode == NDARRAY_INT16)) {\n",
" type = 1;\n",
" } else {\n",
" type = 2;\n",
" }\n",
"\n",
" size_t i, j, k, l;\n",
" l = (n >> 1) + 1;\n",
" k = n;\n",
" uint8_t ui8_tmp, *ui8_array = NULL;\n",
" uint16_t ui16_tmp, *ui16_array = NULL;\n",
" mp_float_t float_tmp, *float_array = NULL;\n",
" \n",
" for(;;) {\n",
" if (l > 1) {\n",
" if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {\n",
" ui8_tmp = ui8_array[--l];\n",
" }\n",
" } else {\n",
" tmp = array[k];\n",
" array[k] = array[1];\n",
" if (--k == 1) {\n",
" array[1] = tmp;\n",
" break;\n",
" }\n",
" }\n",
" i = l;\n",
" j = 2 * l;\n",
" while(j <= k) {\n",
" if(j < k && array[j] < array[j+1]) {\n",
" j++;\n",
" }\n",
" } \n",
" if(tmp < array[j]) {\n",
" array[i] = array[j];\n",
" i = j;\n",
" j <<= 1;\n",
" } else break;\n",
" }\n",
" array[i] = tmp;\n",
"}\n",
"\n",
"mp_obj_t numerical_sort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\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_n, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 1 } },\n",
" { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = -1 } },\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",
" \n",
" if(!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(\"sort argument must be an ndarray\");\n",
" }\n",
" return mp_const_none;\n",
"}*/"
]
},
{
@ -7898,11 +7960,11 @@
},
{
"cell_type": "code",
"execution_count": 270,
"execution_count": 307,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-31T11:13:22.819965Z",
"start_time": "2019-10-31T11:13:22.813107Z"
"end_time": "2019-10-31T22:09:52.722470Z",
"start_time": "2019-10-31T22:09:52.717256Z"
}
},
"outputs": [
@ -7920,11 +7982,11 @@
},
{
"cell_type": "code",
"execution_count": 292,
"execution_count": 381,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-31T11:49:48.565012Z",
"start_time": "2019-10-31T11:49:47.525535Z"
"end_time": "2019-11-01T13:12:54.652197Z",
"start_time": "2019-11-01T13:12:53.380723Z"
},
"scrolled": false
},
@ -7941,13 +8003,14 @@
"GEN build/genhdr/qstrdefs.collected.h\n",
"QSTR not updated\n",
"CC ../../py/objmodule.c\n",
"CC ../../../ulab/code/fft.c\n",
"CC ../../../ulab/code/numerical.c\n",
"CC ../../../ulab/code/ulab.c\n",
"LINK micropython\n",
" text\t data\t bss\t dec\t hex\tfilename\n",
" 2117\t 6862\t 0\t 8979\t 2313\tbuild/build/frozen_mpy.o\n",
" 34\t 0\t 0\t 34\t 22\tbuild/build/frozen.o\n",
" 495084\t 58344\t 2104\t 555532\t 87a0c\tmicropython\n"
" 494949\t 58328\t 2104\t 555381\t 87975\tmicropython\n"
]
}
],