jepler/circuitpython-ulab
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
circuitpython-ulab/code/ndarray.c
Zoltán Vörös 05b6e0a464 ndim code backup without any particular milestone
2019-12-24 08:17:28 +01:00

1015 lines
44 KiB
C
Raw Permalink Blame History

/*
* This file is part of the micropython-ulab project,
*
* https://github.com/v923z/micropython-ulab
*
* The MIT License (MIT)
*
* Copyright (c) 2019 Zoltán Vörös
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "py/runtime.h"
#include "py/binary.h"
#include "py/obj.h"
#include "py/objtuple.h"
#include "ndarray.h"
// This function is copied verbatim from objarray.c
STATIC mp_obj_array_t *array_new(char typecode, size_t n) {
int typecode_size = mp_binary_get_size('@', typecode, NULL);
mp_obj_array_t *o = m_new_obj(mp_obj_array_t);
// this step could probably be skipped: we are never going to store a bytearray per se
#if MICROPY_PY_BUILTINS_BYTEARRAY && MICROPY_PY_ARRAY
o->base.type = (typecode == BYTEARRAY_TYPECODE) ? &mp_type_bytearray : &mp_type_array;
#elif MICROPY_PY_BUILTINS_BYTEARRAY
o->base.type = &mp_type_bytearray;
#else
o->base.type = &mp_type_array;
#endif
o->typecode = typecode;
o->free = 0;
o->len = n;
o->items = m_new(byte, typecode_size * o->len);
return o;
}
// helper functions
mp_float_t ndarray_get_float_value(void *data, uint8_t typecode, size_t index) {
if(typecode == NDARRAY_UINT8) {
return (mp_float_t)((uint8_t *)data)[index];
} else if(typecode == NDARRAY_INT8) {
return (mp_float_t)((int8_t *)data)[index];
} else if(typecode == NDARRAY_UINT16) {
return (mp_float_t)((uint16_t *)data)[index];
} else if(typecode == NDARRAY_INT16) {
return (mp_float_t)((int16_t *)data)[index];
} else {
return (mp_float_t)((mp_float_t *)data)[index];
}
}
size_t ndarray_index_from_contracted(size_t index, ndarray_obj_t *ndarray, int32_t *strides, uint8_t ndim, uint8_t axis) {
// calculates the index in the original (linear) array from the index in the contracted (linear) array
size_t q, new_index = 0;
for(size_t i=0; i <= ndim-1; i++) {
q = index / strides[i];
if(i < axis) {
new_index += q * ndarray->strides[i];
} else {
new_index += q * ndarray->strides[i+1];
}
index -= q * strides[i];
}
return new_index + ndarray->offset;
}
size_t *ndarray_new_coords(uint8_t ndim) {
size_t *coords = m_new(size_t, ndim);
memset(coords, 0, ndim*sizeof(size_t));
return coords;
}
void ndarray_print(const mp_print_t *print, mp_obj_t self_in, mp_print_kind_t kind) {
(void)kind;
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
size_t offset = self->offset;
uint8_t print_extra = self->ndim;
size_t *coords = ndarray_new_coords(self->ndim);
mp_print_str(print, "array(");
for(size_t i=0; i < self->len; i++) {
for(uint8_t j=0; j < print_extra; j++) {
mp_print_str(print, "[");
}
print_extra = 0;
if(!self->boolean) {
mp_obj_print_helper(print, mp_binary_get_val_array(self->array->typecode, self->array->items, offset), PRINT_REPR);
} else {
if(((uint8_t *)self->array->items)[offset]) {
mp_print_str(print, "True");
} else {
mp_print_str(print, "False");
}
}
offset += self->strides[self->ndim-1];
coords[self->ndim-1] += 1;
if(coords[self->ndim-1] != self->shape[self->ndim-1]) {
mp_print_str(print, ", ");
}
for(uint8_t j=self->ndim-1; j > 0; j--) {
if(coords[j] == self->shape[j]) {
offset -= self->shape[j] * self->strides[j];
offset += self->strides[j-1];
print_extra += 1;
coords[j] = 0;
coords[j-1] += 1;
mp_print_str(print, "]");
} else { // coordinates can change only, if the last coordinate changes
break;
}
}
if(print_extra && (i != self->len-1)) {
mp_print_str(print, ",\n");
if(print_extra > 1) {
mp_print_str(print, "\n");
}
}
}
mp_print_str(print, "]");
m_del(size_t, coords, self->ndim);
if(self->boolean) {
mp_print_str(print, ", dtype=bool)");
} else if(self->array->typecode == NDARRAY_UINT8) {
mp_print_str(print, ", dtype=uint8)");
} else if(self->array->typecode == NDARRAY_INT8) {
mp_print_str(print, ", dtype=int8)");
} else if(self->array->typecode == NDARRAY_UINT16) {
mp_print_str(print, ", dtype=uint16)");
} else if(self->array->typecode == NDARRAY_INT16) {
mp_print_str(print, ", dtype=int16)");
} else if(self->array->typecode == NDARRAY_FLOAT) {
mp_print_str(print, ", dtype=float)");
}
}
ndarray_obj_t *ndarray_new_ndarray(uint8_t ndim, size_t *shape, int32_t *strides, uint8_t typecode) {
// Creates the base ndarray with shape, and initialises the values to straight 0s
ndarray_obj_t *ndarray = m_new_obj(ndarray_obj_t);
ndarray->base.type = &ulab_ndarray_type;
ndarray->ndim = ndim;
ndarray->shape = shape;
ndarray->strides = strides;
ndarray->offset = 0;
if(typecode == NDARRAY_BOOL) {
typecode = NDARRAY_UINT8;
ndarray->boolean = NDARRAY_BOOLEAN;
} else {
ndarray->boolean = NDARRAY_NUMERIC;
}
ndarray->len = 1;
for(uint8_t i=0; i < ndim; i++) {
ndarray->len *= shape[i];
}
mp_obj_array_t *array = array_new(typecode, ndarray->len);
// this should set all elements to 0, irrespective of the of the typecode (all bits are zero)
// we could, perhaps, leave this step out, and initialise the array only, when needed
memset(array->items, 0, array->len);
ndarray->array = array;
return ndarray;
}
ndarray_obj_t *ndarray_new_dense_ndarray(uint8_t ndim, size_t *shape, uint8_t typecode) {
// creates a dense array, i.e., one, where the strides are derived directly from the shapes
int32_t *strides = m_new(int32_t, ndim);
strides[ndim-1] = 1;
for(size_t i=ndim-1; i > 0; i--) {
strides[i-1] = strides[i] * shape[i];
}
return ndarray_new_ndarray(ndim, shape, strides, typecode);
}
ndarray_obj_t *ndarray_new_view(mp_obj_array_t *array, uint8_t ndim, size_t *shape, int32_t *strides, size_t offset, uint8_t boolean) {
ndarray_obj_t *ndarray = m_new_obj(ndarray_obj_t);
ndarray->base.type = &ulab_ndarray_type;
ndarray->boolean = boolean;
ndarray->ndim = ndim;
ndarray->shape = shape;
ndarray->strides = strides;
ndarray->len = 1;
for(uint8_t i=0; i < ndim; i++) {
ndarray->len *= shape[i];
}
ndarray->offset = offset;
ndarray->array = array;
return ndarray;
}
ndarray_obj_t *ndarray_copy_view(ndarray_obj_t *input, uint8_t typecode) {
// Creates a new ndarray from the input
// If the input was a sliced view, the output will inherit the shape, but not the strides
int32_t *strides = m_new(int32_t, input->ndim);
strides[input->ndim-1] = 1;
for(uint8_t i=input->ndim-1; i > 0; i--) {
strides[i-1] = strides[i] * input->shape[i];
}
ndarray_obj_t *ndarray = ndarray_new_ndarray(input->ndim, input->shape, strides, typecode);
ndarray->boolean = input->boolean;
mp_obj_t item;
size_t offset = input->offset;
size_t *coords = ndarray_new_coords(input->ndim);
for(size_t i=0; i < ndarray->len; i++) {
item = mp_binary_get_val_array(input->array->typecode, input->array->items, offset);
mp_binary_set_val_array(typecode, ndarray->array->items, i, item);
offset += input->strides[input->ndim-1];
coords[input->ndim-1] += 1;
for(uint8_t j=ndarray->ndim-1; j > 0; j--) {
if(coords[j] == input->shape[j]) {
offset -= input->shape[j] * input->strides[j];
offset += input->strides[j-1];
coords[j] = 0;
coords[j-1] += 1;
} else { // coordinates can change only, if the last coordinate changes
break;
}
}
}
m_del(size_t, coords, input->ndim);
return ndarray;
}
ndarray_obj_t *ndarray_new_linear_array(size_t len, uint8_t dtype) {
size_t *shape = m_new(size_t, 1);
int32_t *strides = m_new(int32_t, 1);
shape[0] = len;
strides[0] = 1;
return ndarray_new_ndarray(1, shape, strides, dtype);
}
STATIC uint8_t ndarray_init_helper(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_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT } },
};
mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
uint8_t dtype = args[1].u_int;
// at this point, dtype can still be `?` for Boolean arrays
return dtype;
}
mp_obj_t ndarray_make_new(const mp_obj_type_t *type, size_t n_args, size_t n_kw, const mp_obj_t *args) {
// TODO: implement dtype, and copy keywords
mp_arg_check_num(n_args, n_kw, 1, 2, true);
mp_map_t kw_args;
mp_map_init_fixed_table(&kw_args, n_kw, args + n_args);
uint8_t dtype = ndarray_init_helper(n_args, args, &kw_args);
mp_obj_t len_in = mp_obj_len_maybe(args[0]);
size_t len = MP_OBJ_SMALL_INT_VALUE(len_in);
ndarray_obj_t *self, *ndarray;
if(MP_OBJ_IS_TYPE(args[0], &ulab_ndarray_type)) {
ndarray = MP_OBJ_TO_PTR(args[0]);
self = ndarray_copy_view(ndarray, dtype);
return MP_OBJ_FROM_PTR(self);
}
// work with a single dimension for now
self = ndarray_new_linear_array(len, dtype);
size_t i=0;
mp_obj_iter_buf_t iter_buf;
mp_obj_t item, iterable = mp_getiter(args[0], &iter_buf);
if(self->boolean) {
uint8_t *array = (uint8_t *)self->array->items;
while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
if(mp_obj_get_float(item)) {
*array = 1;
}
array++;
}
} else {
while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
mp_binary_set_val_array(dtype, self->array->items, i++, item);
}
}
return MP_OBJ_FROM_PTR(self);
}
mp_bound_slice_t generate_slice(mp_uint_t n, mp_obj_t index) {
// micropython seems to have difficulties with negative steps
mp_bound_slice_t slice;
if(MP_OBJ_IS_TYPE(index, &mp_type_slice)) {
mp_seq_get_fast_slice_indexes(n, index, &slice);
} else if(mp_obj_is_int(index)) {
int32_t _index = mp_obj_get_int(index);
if(_index < 0) {
_index += n;
}
if((_index >= n) || (_index < 0)) {
mp_raise_msg(&mp_type_IndexError, "index is out of bounds");
}
slice.start = _index;
slice.stop = _index + 1;
slice.step = 1;
} else {
mp_raise_msg(&mp_type_IndexError, "indices must be integers, slices, or Boolean lists");
}
return slice;
}
size_t slice_length(mp_bound_slice_t slice) {
int32_t len, correction = 1;
if(slice.step > 0) correction = -1;
len = (slice.stop - slice.start + (slice.step + correction)) / slice.step;
if(len < 0) return 0;
return (size_t)len;
}
ndarray_obj_t *ndarray_new_view_from_tuple(ndarray_obj_t *ndarray, mp_obj_tuple_t *slices) {
mp_bound_slice_t slice;
size_t *shape_array = m_new(size_t, ndarray->ndim);
int32_t *strides_array = m_new(int32_t, ndarray->ndim);
size_t offset = ndarray->offset;
for(uint8_t i=0; i < ndarray->ndim; i++) {
if(i < slices->len) {
slice = generate_slice(ndarray->shape[i], slices->items[i]);
offset += ndarray->offset + slice.start * ndarray->strides[i];
shape_array[i] = slice_length(slice);
strides_array[i] = ndarray->strides[i] * slice.step;
} else {
shape_array[i] = ndarray->shape[i];
strides_array[i] = ndarray->strides[i];
}
}
return ndarray_new_view(ndarray->array, ndarray->ndim, shape_array, strides_array, offset, ndarray->boolean);
}
bool ndarray_check_compatibility(ndarray_obj_t *lhs, ndarray_obj_t *rhs) {
if(rhs->ndim > lhs->ndim) {
return false;
}
for(uint8_t i=0; i < rhs->ndim; i++) {
if((rhs->shape[rhs->ndim-1-i] != 1) && (rhs->shape[rhs->ndim-1-i] != lhs->shape[lhs->ndim-1-i])) {
return false;
}
}
return true;
}
mp_obj_t ndarray_assign_view_from_tuple(ndarray_obj_t *ndarray, mp_obj_tuple_t *slices, mp_obj_t value) {
ndarray_obj_t *lhs = ndarray_new_view_from_tuple(ndarray, slices);
ndarray_obj_t *rhs;
if(MP_OBJ_IS_TYPE(value, &ulab_ndarray_type)) {
rhs = MP_OBJ_TO_PTR(value);
// since this is an assignment, the left hand side should definitely be able to contain the right hand side
if(!ndarray_check_compatibility(lhs, rhs)) {
mp_raise_ValueError("could not broadcast input array into output array");
}
} else { // we have a scalar, so create an ndarray for it
size_t *shape = m_new(size_t, lhs->ndim*sizeof(size_t));
for(uint8_t i=0; i < lhs->ndim; i++) {
shape[i] = 1;
}
rhs = ndarray_new_dense_ndarray(lhs->ndim, shape, lhs->array->typecode);
mp_binary_set_val_array(rhs->array->typecode, rhs->array->items, 0, value);
}
size_t roffset = rhs->offset;
size_t loffset = lhs->offset;
size_t *lcoords = ndarray_new_coords(lhs->ndim);
mp_obj_t item;
uint8_t diff_ndim = lhs->ndim - rhs->ndim;
for(size_t i=0; i < lhs->len; i++) {
item = mp_binary_get_val_array(rhs->array->typecode, rhs->array->items, roffset);
mp_binary_set_val_array(lhs->array->typecode, lhs->array->items, loffset, item);
for(uint8_t j=lhs->ndim-1; j > 0; j--) {
loffset += lhs->strides[j];
lcoords[j] += 1;
if(j >= diff_ndim) {
if(rhs->shape[j-diff_ndim] != 1) {
roffset += rhs->strides[j-diff_ndim];
}
}
if(lcoords[j] == lhs->shape[j]) { // we are at a dimension boundary
if(j > diff_ndim) { // this means right-hand-side coordinates that haven't been prepended
if(rhs->shape[j-diff_ndim] != 1) {
// if rhs->shape[j-diff_ndim] != 1, then rhs->shape[j-diff_ndim] == lhs->shape[j],
// so we can advance the offset counter
roffset -= rhs->shape[j-diff_ndim] * rhs->strides[j-diff_ndim];
roffset += rhs->strides[j-diff_ndim-1];
} else { // rhs->shape[j-diff_ndim] == 1
roffset -= rhs->strides[j-diff_ndim];
}
} else {
roffset = rhs->offset;
}
loffset -= lhs->shape[j] * lhs->strides[j];
loffset += lhs->strides[j-1];
lcoords[j] = 0;
lcoords[j-1] += 1;
} else { // coordinates can change only, if the last coordinate changes
break;
}
}
}
m_del(size_t, lcoords, lhs->ndim);
return mp_const_none;
}
mp_obj_t ndarray_subscr(mp_obj_t self, mp_obj_t index, mp_obj_t value) {
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(self);
if (value == MP_OBJ_SENTINEL) { // return value(s)
if(mp_obj_is_int(index) || MP_OBJ_IS_TYPE(index, &mp_type_slice)) {
mp_obj_t *items = m_new(mp_obj_t, 1);
items[0] = index;
mp_obj_t tuple = mp_obj_new_tuple(1, items);
return ndarray_new_view_from_tuple(ndarray, tuple);
}
// first, check, whether all members of the tuple are integer scalars, or slices
if(MP_OBJ_IS_TYPE(index, &mp_type_tuple)) {
mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(index);
for(uint8_t i=0; i < tuple->len; i++) {
if(!MP_OBJ_IS_TYPE(tuple->items[i], &mp_type_slice) && !mp_obj_is_int(tuple->items[i])) {
// TODO: we have to return a copy here
mp_raise_msg(&mp_type_IndexError, "wrong index type");
}
}
// now we know that we can return a view
ndarray_obj_t *result = ndarray_new_view_from_tuple(ndarray, tuple);
return MP_OBJ_FROM_PTR(result);
}
} else { // assignment; the value must be an ndarray, or a scalar
if(!MP_OBJ_IS_TYPE(value, &ulab_ndarray_type) && !mp_obj_is_int(value) && !mp_obj_is_float(value)) {
mp_raise_ValueError("right hand side must be an ndarray, or a scalar");
} else {
if(mp_obj_is_int(index) || MP_OBJ_IS_TYPE(index, &mp_type_slice)) {
mp_obj_t *items = m_new(mp_obj_t, 1);
items[0] = index;
mp_obj_t tuple = mp_obj_new_tuple(1, items);
return ndarray_assign_view_from_tuple(ndarray, tuple, value);
}
if(MP_OBJ_IS_TYPE(index, &mp_type_tuple)) {
mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(index);
for(uint8_t i=0; i < tuple->len; i++) {
if(!MP_OBJ_IS_TYPE(tuple->items[i], &mp_type_slice) && !mp_obj_is_int(tuple->items[i])) {
// TODO: we have to return a copy here
mp_raise_msg(&mp_type_IndexError, "wrong index type");
}
}
return ndarray_assign_view_from_tuple(ndarray, tuple, value);
} else {
mp_raise_NotImplementedError("wrong index type");
}
}
}
return mp_const_none;
}
// itarray iterator
mp_obj_t ndarray_getiter(mp_obj_t o_in, mp_obj_iter_buf_t *iter_buf) {
return mp_obj_new_ndarray_iterator(o_in, 0, iter_buf);
}
typedef struct _mp_obj_ndarray_it_t {
mp_obj_base_t base;
mp_fun_1_t iternext;
mp_obj_t ndarray;
size_t cur;
} mp_obj_ndarray_it_t;
mp_obj_t ndarray_iternext(mp_obj_t self_in) {
mp_obj_ndarray_it_t *self = MP_OBJ_TO_PTR(self_in);
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(self->ndarray);
size_t iter_end = ndarray->shape[0];
if(self->cur < iter_end) {
if(ndarray->ndim == 1) { // we have a linear array
// read the current value
self->cur++;
return mp_binary_get_val_array(ndarray->array->typecode, ndarray->array->items, self->cur-1);
} else { // we have a tensor, return the reduced view
size_t offset = ndarray->offset + self->cur * ndarray->strides[0];
ndarray_obj_t *value = ndarray_new_view(ndarray->array, ndarray->ndim-1, ndarray->shape+1, ndarray->strides+1, offset, ndarray->boolean);
self->cur++;
return MP_OBJ_FROM_PTR(value);
}
} else {
return MP_OBJ_STOP_ITERATION;
}
}
mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t ndarray, size_t cur, mp_obj_iter_buf_t *iter_buf) {
assert(sizeof(mp_obj_ndarray_it_t) <= sizeof(mp_obj_iter_buf_t));
mp_obj_ndarray_it_t *o = (mp_obj_ndarray_it_t*)iter_buf;
o->base.type = &mp_type_polymorph_iter;
o->iternext = ndarray_iternext;
o->ndarray = ndarray;
o->cur = cur;
return MP_OBJ_FROM_PTR(o);
}
mp_obj_t ndarray_shape(mp_obj_t self_in) {
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
mp_obj_t *items = m_new(mp_obj_t, self->ndim);
size_t *shape = (size_t *)self->shape;
for(uint8_t i=0; i < self->ndim; i++) {
items[i] = mp_obj_new_int(shape[i]);
}
mp_obj_t tuple = mp_obj_new_tuple(self->ndim, items);
m_del(mp_obj_t, items, self->ndim);
return tuple;
}
mp_obj_t ndarray_strides(mp_obj_t self_in) {
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
mp_obj_t *items = m_new(mp_obj_t, self->ndim);
int32_t *strides = (int32_t *)self->strides;
for(int8_t i=0; i < self->ndim; i++) {
items[i] = mp_obj_new_int(strides[i]);
}
mp_obj_t tuple = mp_obj_new_tuple(self->ndim, items);
m_del(mp_obj_t, items, self->ndim);
return tuple;
}
mp_obj_t ndarray_info(mp_obj_t self_in) {
// TODO: the proper way of handling this would be to use mp_print_str()
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
printf("class: ndarray\n");
printf("shape: (%ld", (size_t)self->shape[0]);
for(uint8_t i=1; i < self->ndim; i++) {
printf(", %ld", (size_t)self->shape[i]);
}
printf(")");
printf("\nstrides: (%ld", (size_t)self->strides[0]);
for(uint8_t i=1; i < self->ndim; i++) {
printf(", %ld", (size_t)self->strides[i]);
}
printf(")");
printf("\nitemsize: %ld\n", mp_binary_get_size('@', self->array->typecode, NULL));
printf("data pointer: %p\n", self->array->items);
if(self->array->typecode == NDARRAY_BOOL) {
printf("type: bool\n");
} else if(self->array->typecode == NDARRAY_UINT8) {
printf("type: uint8\n");
} else if(self->array->typecode == NDARRAY_INT8) {
printf("type: int8\n");
} else if(self->array->typecode == NDARRAY_UINT16) {
printf("type: uint16\n");
} else if(self->array->typecode == NDARRAY_INT16) {
printf("type: int16\n");
} else {
printf("type: float\n");
}
return mp_const_none;
}
mp_obj_t ndarray_flatten(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
static const mp_arg_t allowed_args[] = {
{ MP_QSTR_order, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_QSTR(MP_QSTR_C)} },
};
mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
mp_arg_parse_all(n_args - 1, pos_args + 1, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(pos_args[0]);
GET_STR_DATA_LEN(args[0].u_obj, order, clen);
if((clen != 1) || ((memcmp(order, "C", 1) != 0) && (memcmp(order, "F", 1) != 0))) {
mp_raise_ValueError("flattening order must be either 'C', or 'F'");
}
ndarray_obj_t *result = ndarray_new_linear_array(ndarray->len, ndarray->array->typecode);
uint8_t itemsize = mp_binary_get_size('@', ndarray->array->typecode, NULL);
if(ndarray->len == ndarray->array->len) { // this is a dense array, we can simply copy everything
if(memcmp(order, "C", 1) == 0) { // C order; this should be fast, because no re-ordering is required
memcpy(result->array->items, ndarray->array->items, itemsize*ndarray->len);
} else { // Fortran order
mp_raise_NotImplementedError("flatten is implemented in C order only");
}
return MP_OBJ_FROM_PTR(result);
} else {
uint8_t *rarray = (uint8_t *)result->array->items;
uint8_t *narray = (uint8_t *)ndarray->array->items;
size_t *coords = ndarray_new_coords(ndarray->ndim);
size_t offset = ndarray->offset;
if(memcmp(order, "C", 1) == 0) { // C order; this is a view, so we have to collect the items
for(size_t i=0; i < result->len; i++) {
memcpy(rarray, &narray[offset*itemsize], itemsize);
rarray += itemsize;
offset += ndarray->strides[ndarray->ndim-1];
coords[ndarray->ndim-1] += 1;
for(uint8_t j=ndarray->ndim-1; j > 0; j--) {
if(coords[j] == ndarray->shape[j]) {
offset -= ndarray->shape[j] * ndarray->strides[j];
offset += ndarray->strides[j-1];
coords[j] = 0;
coords[j-1] += 1;
} else { // coordinates can change only, if the last coordinate changes
break;
}
}
}
m_del(size_t, coords, ndarray->ndim);
} else { // Fortran order
mp_raise_NotImplementedError("flatten is implemented for C order only");
}
//m_del(int32_t, shape_strides, 1);
return MP_OBJ_FROM_PTR(result);
}
}
mp_obj_t ndarray_reshape(mp_obj_t oin, mp_obj_t new_shape) {
// TODO: not all reshaping operations can be realised via views!
// returns a new view with the specified shape
ndarray_obj_t *ndarray_in = MP_OBJ_TO_PTR(oin);
mp_obj_tuple_t *shape = MP_OBJ_TO_PTR(new_shape);
if(ndarray_in->len != ndarray_in->array->len) {
mp_raise_ValueError("input and output shapes are not compatible");
}
size_t *shape_array = m_new(size_t, shape->len);
int32_t *strides_array = m_new(int32_t, shape->len);
size_t new_offset = ndarray_in->offset; // this has to be re-calculated
strides_array[shape->len-1] = 1;
shape_array[shape->len-1] = mp_obj_get_int(shape->items[shape->len-1]);
for(uint8_t i=shape->len-1; i > 0; i--) {
shape_array[i-1] = mp_obj_get_int(shape->items[i-1]);
strides_array[i-1] = strides_array[i] * shape_array[i];
}
ndarray_obj_t *ndarray = ndarray_new_view(ndarray_in->array, shape->len, shape_array, strides_array, new_offset, ndarray_in->boolean);
return MP_OBJ_FROM_PTR(ndarray);
}
mp_obj_t ndarray_transpose(mp_obj_t self_in) {
// TODO: check, what happens to the offset here!
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
size_t *shape = m_new(size_t, self->ndim);
int32_t *strides = m_new(int32_t, self->ndim);
for(uint8_t i=0; i < self->ndim; i++) {
shape[i] = self->shape[self->ndim-1-i];
strides[i] = self->strides[self->ndim-1-i];
}
ndarray_obj_t *ndarray = ndarray_new_view(self->array, self->ndim, shape, strides, self->offset, self->boolean);
return MP_OBJ_FROM_PTR(ndarray);
}
mp_obj_t ndarray_itemsize(mp_obj_t self_in) {
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
return mp_obj_new_int(mp_binary_get_size('@', self->array->typecode, NULL));
}
uint8_t upcasting(ndarray_obj_t *in1, ndarray_obj_t *in2) {
if((in1->array->typecode = NDARRAY_FLOAT) || (in2->array->typecode = NDARRAY_FLOAT)) { // 9 cases
return NDARRAY_FLOAT;
} else if(in1->array->typecode == in2->array->typecode) { // 4 cases
return in1->array->typecode;
} else if(in1->array->typecode == NDARRAY_UINT8) { // 3 cases
if(in2->array->typecode == NDARRAY_UINT16) return NDARRAY_UINT16;
else return NDARRAY_INT16; // in2->array->typecode == NDARRAY_INT8, NDARRAY_INT16,
} else if(in1->array->typecode == NDARRAY_INT8) { // 3 cases
return NDARRAY_INT16; // in2->array->typecode == NDARRAY_UINT8, NDARRAY_UINT16, NDARRAY_INT16
} else if(in1->array->typecode == NDARRAY_UINT16) { // 3 cases
if((in2->array->typecode == NDARRAY_UINT8) || (in2->array->typecode == NDARRAY_INT8)) return NDARRAY_UINT16;
return NDARRAY_FLOAT; // in2->array->typecode == NDARRAY_INT16
} else if(in1->array->typecode == NDARRAY_INT16) { // 3 cases
if((in2->array->typecode == NDARRAY_UINT8) || (in2->array->typecode == NDARRAY_INT8)) return NDARRAY_INT16;
return NDARRAY_FLOAT; // in2->array->typecode == NDARRAY_UINT16
}
}
size_t *broadcasting(ndarray_obj_t *in1, ndarray_obj_t *in2) {
// creates a new dense array with the dimensions
// and shapes dictated by the broadcasting rules
uint8_t ndim = in2->ndim;
if(in1->ndim > in2->ndim) {
// overwrite ndim, if the first array is larger than the second
ndim = in1->ndim;
}
size_t *shape = m_new(size_t, ndim);
size_t *shape1 = m_new(size_t, ndim);
size_t *shape2 = m_new(size_t, ndim);
for(uint8_t i=0; i < ndim; i++) {
// initialise the shapes with straight ones
shape1[i] = shape2[i] = 1;
shape[i] = 0;
}
// overwrite the first in1->ndim (in2->ndim) shape values from the right
for(uint8_t i=ndim; i > ndim-in1->ndim; i--) {
shape1[i-1] = in1->shape[i-1];
}
for(uint8_t i=ndim; i > ndim-in2->ndim; i--) {
shape2[i-1] = in2->shape[i-1];
}
// check, whether the new shapes conform to the broadcasting rules
for(uint8_t i=0; i < ndim; i++) {
if((shape1[i] == shape2[i]) || (shape1[i] == 1) || (shape2[i] == 1)) {
shape[i] = shape1[i];
if(shape1[i] < shape2[i]) {
shape[i] = shape2[i];
}
} else {
m_del(size_t, shape1, ndim);
m_del(size_t, shape2, ndim);
break;
}
}
m_del(size_t, shape1, ndim);
m_del(size_t, shape2, ndim);
return shape;
}
ndarray_obj_t *ndarray_do_binary_op(ndarray_obj_t *lhs, ndarray_obj_t *rhs, uint8_t op) {
uint8_t typecode = upcasting(lhs, rhs);
size_t *shape = broadcasting(lhs, rhs);
uint8_t ndim = lhs->ndim;
if(rhs->ndim > lhs->ndim) ndim = rhs->ndim;
// this is going to be the result array
ndarray_obj_t *ndarray = ndarray_new_dense_ndarray(ndim, shape, typecode);
size_t offset = ndarray->offset;
size_t roffset = rhs->offset;
size_t loffset = lhs->offset;
size_t *coords = ndarray_new_coords(ndim); // coordinates in ndarray
size_t *lcoords = ndarray_new_coords(ndim); // coordinates in the left-hand-side ndarray
size_t *rcoords = ndarray_new_coords(ndim); // coordinates in the right-hand-side ndarray
int32_t *lstride = m_new(int32_t, ndim);
int32_t *rstride = m_new(int32_t, ndim);
for(uint8_t i=0; i < lhs->ndim; i++) {
lstride[ndim-1-i] = lhs->strides[lhs->ndim-1-i];
}
for(uint8_t i=lhs->ndim; i < ndim; i++) {
lstride[ndim-1-i] = lhs->strides[0];
}
for(uint8_t i=0; i < rhs->ndim; i++) {
rstride[ndim-1-i] = rhs->strides[rhs->ndim-1-i];
}
for(uint8_t i=rhs->ndim; i < ndim; i++) {
rstride[ndim-1-i] = rhs->strides[0];
}
mp_float_t lvalue, rvalue, result = 0.0;
int32_t ilvalue, irvalue, iresult = 0;
uint8_t float_size = sizeof(mp_float_t);
void *narray = ndarray->array->items;
uint8_t *larray = (uint8_t *)lhs->array->items;
uint8_t *rarray = (uint8_t *)rhs->array->items;
for(size_t i=0; i < ndarray->len; i++) {
// we could make this prettier by moving everything into a function,
// but the function call might be too expensive, especially that it is in the loop...
// left hand side
if(lhs->array->typecode == NDARRAY_UINT8) {
ilvalue = (int32_t)((uint8_t *)larray)[loffset];
} else if(lhs->array->typecode == NDARRAY_INT8) {
ilvalue = (int32_t)((int8_t *)larray)[loffset];
} else if(lhs->array->typecode == NDARRAY_UINT16) {
ilvalue = (int32_t)((uint16_t *)larray)[loffset*2];
} else if(lhs->array->typecode == NDARRAY_INT16) {
ilvalue = (int32_t)((int16_t *)larray)[loffset*2];
} else {
lvalue = (mp_float_t)((mp_float_t *)larray)[loffset*float_size];
}
// right hand side
if(rhs->array->typecode == NDARRAY_UINT8) {
irvalue = (int32_t)((uint8_t *)rarray)[roffset];
} else if(rhs->array->typecode == NDARRAY_INT8) {
irvalue = (int32_t)((int8_t *)rarray)[roffset];
} else if(rhs->array->typecode == NDARRAY_UINT16) {
irvalue = (int32_t)((uint16_t *)rarray)[roffset*2];
} else if(lhs->array->typecode == NDARRAY_INT16) {
irvalue = (int32_t)((int16_t *)rarray)[roffset*2];
} else {
rvalue = (mp_float_t)((mp_float_t *)rarray)[roffset*float_size];
}
// this is the place, where the actual operations take place
if(op == MP_BINARY_OP_ADD) {
result = lvalue + rvalue;
} else if(op == MP_BINARY_OP_SUBTRACT) {
result = lvalue - rvalue;
} else if(op == MP_BINARY_OP_MULTIPLY) {
result = lvalue * rvalue;
}
// cast the result to the proper output type
if(typecode == NDARRAY_UINT8) {
((uint8_t *)narray)[offset] = (uint8_t)result;
} else if(typecode == NDARRAY_INT8) {
((int8_t *)narray)[offset] = (int8_t)result;
} else if(typecode == NDARRAY_UINT16) {
((uint16_t *)narray)[offset] = (uint16_t)result;
} else if(typecode == NDARRAY_INT16) {
((int16_t *)narray)[offset] = (int16_t)result;
} else {
((mp_float_t *)narray)[offset] = result;
}
offset += ndarray->strides[ndim-1];
coords[ndim-1] += 1;
for(uint8_t j=ndim-1; j > 0; j--) {
if(coords[j] == ndarray->shape[j]) {
offset -= ndarray->shape[j] * ndarray->strides[j];
offset += ndarray->strides[j-1];
coords[j] = 0;
coords[j-1] += 1;
} else { // coordinates can change only, if the last coordinate changes
break;
}
}
}
m_del(size_t, coords, ndim);
m_del(size_t, lcoords, ndim);
m_del(size_t, rcoords, ndim);
m_del(int32_t, lstride, ndim);
m_del(int32_t, rstride, ndim);
return ndarray;
}
// Binary operations
mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t LHS, mp_obj_t RHS) {
// First, handle the case, when the operand on the right hand side is a scalar
ndarray_obj_t *lhs = MP_OBJ_TO_PTR(LHS);
ndarray_obj_t *rhs;
if(mp_obj_is_int(RHS) || mp_obj_is_float(RHS)) {
size_t *shape = m_new(size_t, lhs->ndim*sizeof(size_t));
for(uint8_t i=0; i < lhs->ndim; i++) {
shape[i] = 1;
}
if(mp_obj_is_int(RHS)) {
int32_t ivalue = mp_obj_get_int(RHS);
if((ivalue > 0) && (ivalue < 256)) {
rhs = ndarray_new_dense_ndarray(lhs->ndim, shape, NDARRAY_UINT8);
} else if((ivalue > 255) && (ivalue < 65535)) {
rhs = ndarray_new_dense_ndarray(lhs->ndim, shape, NDARRAY_UINT16);
} else if((ivalue < 0) && (ivalue > -128)) {
rhs = ndarray_new_dense_ndarray(lhs->ndim, shape, NDARRAY_INT8);
} else if((ivalue < -127) && (ivalue > -32767)) {
rhs = ndarray_new_dense_ndarray(lhs->ndim, shape, NDARRAY_INT16);
} else { // the integer value clearly does not fit the ulab types, so move on to float
rhs = ndarray_new_dense_ndarray(lhs->ndim, shape, NDARRAY_FLOAT);
}
mp_binary_set_val_array(rhs->array->typecode, rhs->array->items, 0, RHS);
} else { // we have a float
rhs = ndarray_new_dense_ndarray(lhs->ndim, shape, NDARRAY_FLOAT);
mp_binary_set_val_array(rhs->array->typecode, rhs->array->items, 0, RHS);
}
} else {
// the right hand side is an ndarray
rhs = MP_OBJ_TO_PTR(RHS);
}
ndarray_obj_t *ndarray = NULL;
size_t *new_shape = broadcasting(lhs, rhs);
switch(op) {
case MP_BINARY_OP_MAT_MULTIPLY:
break;
case MP_BINARY_OP_ADD:
if(new_shape[0] == 0) {
mp_raise_ValueError("operands could not be cast together");
}
ndarray = ndarray_do_binary_op(lhs, rhs, op);
// The parameters of RUN_BINARY_LOOP are
// typecode of type_out, type_left, type_right, out_array, lhs_array, rhs_array, shape, ndim, operator
//RUN_BINARY_LOOP(ndarray, NDARRAY_FLOAT, mp_float_t, mp_float_t, mp_float_t, lhs, rhs, new_shape, ndim, operator);
/* if(lhs->array->typecode == NDARRAY_UINT8) {
if(rhs->array->typecode == NDARRAY_UINT8) {
RUN_BINARY_LOOP(NDARRAY_UINT8, uint8_t, uint8_t, uint8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT8) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_UINT16) {
RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint8_t, uint16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, lhs, rhs, new_shape, ndim, operator);
}
} else if(lhs->array->typecode == NDARRAY_INT8) {
if(rhs->array->typecode == NDARRAY_UINT8) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT8) {
RUN_BINARY_LOOP(NDARRAY_INT8, int8_t, int8_t, int8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_UINT16) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, int16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int8_t, mp_float_t, lhs, rhs, new_shape, ndim, operator);
}
} else if(lhs->array->typecode == NDARRAY_UINT16) {
if(rhs->array->typecode == NDARRAY_UINT8) {
RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT8) {
RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, int8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_UINT16) {
RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, int16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, lhs, rhs, new_shape, ndim, operator);
}
} else if(lhs->array->typecode == NDARRAY_INT16) {
if(rhs->array->typecode == NDARRAY_UINT8) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, uint8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT8) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_UINT16) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int16_t, uint16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, mp_float_t, lhs, rhs, new_shape, ndim, operator);
}
} else if(lhs->array->typecode == NDARRAY_FLOAT) {
if(rhs->array->typecode == NDARRAY_UINT8) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT8) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int8_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_UINT16) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_INT16) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int16_t, lhs, rhs, new_shape, ndim, operator);
} else if(rhs->array->typecode == NDARRAY_FLOAT) {
RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, mp_float_t, lhs, rhs, new_shape, ndim, operator);
}
} else { // this should never happen
mp_raise_TypeError("wrong input type");
} */
break;
default:
// m_del(size_t, shape, lhs->ndim*sizeof(size_t));
return mp_const_none;
}
return MP_OBJ_FROM_PTR(ndarray);
// return mp_const_none;
}
mp_obj_t ndarray_unary_op(mp_unary_op_t op, mp_obj_t self_in) {
ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
ndarray_obj_t *ndarray;
switch(op) {
case MP_UNARY_OP_LEN:
return mp_obj_new_int(self->shape[0]);
break;
case MP_UNARY_OP_INVERT:
if(self->array->typecode == NDARRAY_FLOAT) {
mp_raise_ValueError("operation is not supported for given type");
}
// we can invert the content byte by byte, there is no need to distinguish between different typecodes
ndarray = ndarray_copy_view(self, self->array->typecode);
uint8_t *array = (uint8_t *)ndarray->array->items;
if(self->boolean == NDARRAY_BOOLEAN) {
for(size_t i=0; i < self->len; i++) array[i] = 1 - array[i];
}
else {
for(size_t i=0; i < self->len; i++) array[i] = ~array[i];
}
return MP_OBJ_FROM_PTR(ndarray);
break;
case MP_UNARY_OP_NEGATIVE:
if(self->boolean == NDARRAY_BOOLEAN) {
mp_raise_TypeError("boolean negative '-' is not supported, use the '~' operator instead");
}
ndarray = ndarray_copy_view(self, self->array->typecode);
if(self->array->typecode == NDARRAY_UINT8) {
uint8_t *array = (uint8_t *)ndarray->array->items;
for(size_t i=0; i < self->len; i++) array[i] = -array[i];
} else if(self->array->typecode == NDARRAY_INT8) {
int8_t *array = (int8_t *)ndarray->array->items;
for(size_t i=0; i < self->len; i++) array[i] = -array[i];
} else if(self->array->typecode == NDARRAY_UINT16) {
uint16_t *array = (uint16_t *)ndarray->array->items;
for(size_t i=0; i < self->len; i++) array[i] = -array[i];
} else if(self->array->typecode == NDARRAY_INT16) {
int16_t *array = (int16_t *)ndarray->array->items;
for(size_t i=0; i < self->len; i++) array[i] = -array[i];
} else {
mp_float_t *array = (mp_float_t *)ndarray->array->items;
for(size_t i=0; i < self->len; i++) array[i] = -array[i];
}
return MP_OBJ_FROM_PTR(ndarray);
break;
case MP_UNARY_OP_POSITIVE:
return MP_OBJ_FROM_PTR(ndarray_copy_view(self, self->array->typecode));
case MP_UNARY_OP_ABS:
if((self->array->typecode == NDARRAY_UINT8) || (self->array->typecode == NDARRAY_UINT16)) {
return MP_OBJ_FROM_PTR(ndarray_copy_view(self, self->array->typecode));
}
ndarray = ndarray_copy_view(self, self->array->typecode);
if(self->array->typecode == NDARRAY_INT8) {
int8_t *array = (int8_t *)ndarray->array->items;
for(size_t i=0; i < self->len; i++) {
if(array[i] < 0) array[i] = -array[i];
}
} else if(self->array->typecode == NDARRAY_INT16) {
int16_t *array = (int16_t *)ndarray->array->items;
for(size_t i=0; i < self->len; i++) {
if(array[i] < 0) array[i] = -array[i];
}
} else {
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];
}
}
return MP_OBJ_FROM_PTR(ndarray);
break;
default: return MP_OBJ_NULL; // operator not supported
}
}
Reference in a new issue View git blame Copy permalink