ulab_ndarray_type is extern now

code/fft.c

@@ -87,35 +87,26 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
     ndarray_obj_t *re = MP_OBJ_TO_PTR(arg_re);
     uint16_t len = re->array->len;
     if((len & (len-1)) != 0) {
+        // TODO: pad the input vector, if the length is not a power of 2
         mp_raise_ValueError("input array length must be power of 2");
     }
     
-    ndarray_obj_t *out_re = create_new_ndarray(1, len, NDARRAY_FLOAT);
-    mp_float_t *data_re = (mp_float_t *)out_re->array->items;
+    ndarray_obj_t *ndarray_re = ndarray_new_linear_array(len, NDARRAY_FLOAT);
+    mp_float_t *data_re = (mp_float_t *)ndarray_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((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);
+        data_re[i] = ndarray_get_float_value(re->array->items, re->array->typecode, re->offset+i*re->strides[0]);
     }
-    }
-    ndarray_obj_t *out_im = create_new_ndarray(1, len, NDARRAY_FLOAT);
-    mp_float_t *data_im = (mp_float_t *)out_im->array->items;
+    ndarray_obj_t *ndarray_im = ndarray_new_linear_array(len, NDARRAY_FLOAT);
+    mp_float_t *data_im = (mp_float_t *)ndarray_im->array->items;
 
     if(n_args == 2) {
         ndarray_obj_t *im = MP_OBJ_TO_PTR(arg_im);
-        if (re->array->len != im->array->len) {
+        if (re->len != im->len) {
             mp_raise_ValueError("real and imaginary parts must be of equal length");
         }
-        if(im->array->typecode == NDARRAY_FLOAT) {
-            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);
-            }
+            data_im[i] = ndarray_get_float_value(im->array->items, im->array->typecode, im->offset+i*im->strides[0]);
         }
     }
     if((type == FFT_FFT) || (type == FFT_SPECTRUM)) {
@@ -134,11 +125,11 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
         }
     }
     if(type == FFT_SPECTRUM) {
-        return MP_OBJ_TO_PTR(out_re);
+        return MP_OBJ_TO_PTR(ndarray_re);
     } else {
         mp_obj_t tuple[2];
-        tuple[0] = out_re;
-        tuple[1] = out_im;
+        tuple[0] = ndarray_re;
+        tuple[1] = ndarray_im;
         return mp_obj_new_tuple(2, tuple);
     }
 }

code/fft.h

@@ -20,4 +20,5 @@
 mp_obj_t fft_fft(size_t , const mp_obj_t *);
 mp_obj_t fft_ifft(size_t , const mp_obj_t *);
 mp_obj_t fft_spectrum(size_t , const mp_obj_t *);
+
 #endif

code/linalg.c

@@ -16,98 +16,6 @@
 #include "py/misc.h"
 #include "linalg.h"
 
-mp_obj_t linalg_transpose(mp_obj_t self_in) {
-    ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
-    // the size of a single item in the array
-    uint8_t _sizeof = mp_binary_get_size('@', self->array->typecode, NULL);
-    
-    // NOTE: 
-    // if the matrices are square, we can simply swap items, but 
-    // generic matrices can't be transposed in place, so we have to 
-    // declare a temporary variable
-    
-    // NOTE: 
-    //  In the old matrix, the coordinate (m, n) is m*self->n + n
-    //  We have to assign this to the coordinate (n, m) in the new 
-    //  matrix, i.e., to n*self->m + m (since the new matrix has self->m columns)
-    
-    // one-dimensional arrays can be transposed by simply swapping the dimensions
-    if((self->m != 1) && (self->n != 1)) {
-        uint8_t *c = (uint8_t *)self->array->items;
-        // self->bytes is the size of the bytearray, irrespective of the typecode
-        uint8_t *tmp = m_new(uint8_t, self->bytes);
-        for(size_t m=0; m < self->m; m++) {
-            for(size_t n=0; n < self->n; n++) {
-                memcpy(tmp+_sizeof*(n*self->m + m), c+_sizeof*(m*self->n + n), _sizeof);
-            }
-        }
-        memcpy(self->array->items, tmp, self->bytes);
-        m_del(uint8_t, tmp, self->bytes);
-    } 
-    SWAP(size_t, self->m, self->n);
-    return mp_const_none;
-}
-
-mp_obj_t linalg_reshape(mp_obj_t self_in, mp_obj_t shape) {
-    ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
-    if(!MP_OBJ_IS_TYPE(shape, &mp_type_tuple) || (MP_OBJ_SMALL_INT_VALUE(mp_obj_len_maybe(shape)) != 2)) {
-        mp_raise_ValueError("shape must be a 2-tuple");
-    }
-
-    mp_obj_iter_buf_t iter_buf;
-    mp_obj_t item, iterable = mp_getiter(shape, &iter_buf);
-    uint16_t m, n;
-    item = mp_iternext(iterable);
-    m = mp_obj_get_int(item);
-    item = mp_iternext(iterable);
-    n = mp_obj_get_int(item);
-    if(m*n != self->m*self->n) {
-        // TODO: the proper error message would be "cannot reshape array of size %d into shape (%d, %d)"
-        mp_raise_ValueError("cannot reshape array (incompatible input/output shape)");
-    }
-    self->m = m;
-    self->n = n;
-    return MP_OBJ_FROM_PTR(self);
-}
-
-mp_obj_t linalg_size(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
-    static const mp_arg_t allowed_args[] = {
-        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } },
-        { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },
-    };
-
-    mp_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("size is defined for ndarrays only");
-    } else {
-        ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
-        if(args[1].u_obj == mp_const_none) {
-            return mp_obj_new_int(ndarray->array->len);
-        } else if(mp_obj_is_int(args[1].u_obj)) {
-            uint8_t ax = mp_obj_get_int(args[1].u_obj);
-            if(ax == 0) {
-                if(ndarray->m == 1) {
-                    return mp_obj_new_int(ndarray->n);
-                } else {
-                    return mp_obj_new_int(ndarray->m);                    
-                }
-            } else if(ax == 1) {
-                if(ndarray->m == 1) {
-                    mp_raise_ValueError("tuple index out of range");
-                } else {
-                    return mp_obj_new_int(ndarray->n);
-                }
-            } else {
-                    mp_raise_ValueError("tuple index out of range");                
-            }
-        } else {
-            mp_raise_TypeError("wrong argument type");
-        }
-    }
-}
-
 bool linalg_invert_matrix(mp_float_t *data, size_t N) {
     // returns true, of the inversion was successful, 
     // false, if the matrix is singular
@@ -151,39 +59,41 @@ bool linalg_invert_matrix(mp_float_t *data, size_t N) {
 }
 
 mp_obj_t linalg_inv(mp_obj_t o_in) {
-    // since inv is not a class method, we have to inspect the input argument first
-    if(!MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {
-        mp_raise_TypeError("only ndarrays can be inverted");
-    }
-    ndarray_obj_t *o = MP_OBJ_TO_PTR(o_in);
     if(!MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {
         mp_raise_TypeError("only ndarray objects can be inverted");
     }
-    if(o->m != o->n) {
+    ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(o_in);
+    if(ndarray->ndim != 2) {
+        mp_raise_ValueError("only two-dimensional tensors can be inverted");
+    }
+    if(ndarray->shape[0] != ndarray->shape[1]) {
         mp_raise_ValueError("only square matrices can be inverted");
     }
-    ndarray_obj_t *inverted = create_new_ndarray(o->m, o->n, NDARRAY_FLOAT);
+    size_t *shape = m_new(size_t, 2);
+    shape[0] = shape[1] = ndarray->shape[0];
+    ndarray_obj_t *inverted = ndarray_new_dense_ndarray(2, shape, NDARRAY_FLOAT);
     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
+    for(size_t m=0; m < ndarray->shape[0]; m++) { // rows first
+        for(size_t n=0; n < ndarray->shape[1]; 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] = (mp_float_t)mp_obj_get_float(elem);
+            elem = mp_binary_get_val_array(ndarray->array->typecode, ndarray->array->items, m*ndarray->shape[1]+n);
+            data[m*ndarray->shape[1]+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(mp_float_t, inverted->array->items, o->n*o->n);
+    if(!linalg_invert_matrix(data, ndarray->shape[0])) {
+        // TODO: I am not sure this is needed here. Otherwise, how should we free up the unused RAM of inverted?
+        m_del(mp_float_t, inverted->array->items, ndarray->shape[0]*ndarray->shape[1]);
         mp_raise_ValueError("input matrix is singular");
     }
     return MP_OBJ_FROM_PTR(inverted);
 }
 
 mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) {
+    return mp_const_none;
+    /*
     // TODO: should the results be upcast?
     ndarray_obj_t *m1 = MP_OBJ_TO_PTR(_m1);
     ndarray_obj_t *m2 = MP_OBJ_TO_PTR(_m2);    
@@ -206,7 +116,7 @@ mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) {
             outdata[i*m1->m+j] = sum;
         }
     }
-    return MP_OBJ_FROM_PTR(out);
+    return MP_OBJ_FROM_PTR(out); */
 } 
 
 mp_obj_t linalg_zeros_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t kind) {
@@ -220,19 +130,21 @@ mp_obj_t linalg_zeros_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
     
     uint8_t dtype = args[1].u_int;
     if(!mp_obj_is_int(args[0].u_obj) && !mp_obj_is_type(args[0].u_obj, &mp_type_tuple)) {
-        mp_raise_TypeError("input argument must be an integer or a 2-tuple");
+        mp_raise_TypeError("input argument must be an integer or a tuple");
     }
     ndarray_obj_t *ndarray = NULL;
     if(mp_obj_is_int(args[0].u_obj)) {
         size_t n = mp_obj_get_int(args[0].u_obj);
-        ndarray = create_new_ndarray(1, n, dtype);
+        size_t *shape = m_new(size_t, 1);
+        shape[0] = n;
+        ndarray = ndarray_new_dense_ndarray(1, shape, dtype);
     } else if(mp_obj_is_type(args[0].u_obj, &mp_type_tuple)) {
         mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(args[0].u_obj);
-        if(tuple->len != 2) {
-            mp_raise_TypeError("input argument must be an integer or a 2-tuple");            
+        size_t *shape = m_new(size_t, tuple->len);
+        for(uint8_t i=0; i < tuple->len; i++) {
+            shape[i] = mp_obj_get_int(tuple->items[i]);
         }
-        ndarray = create_new_ndarray(mp_obj_get_int(tuple->items[0]), 
-                                                  mp_obj_get_int(tuple->items[1]), dtype);
+        ndarray = ndarray_new_dense_ndarray(tuple->len, shape, dtype);
     }
     if(kind == 1) {
         mp_obj_t one = mp_obj_new_int(1);
@@ -252,6 +164,9 @@ mp_obj_t linalg_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args)
 }
 
 mp_obj_t linalg_eye(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+    // TODO: this is a bit more problematic in higher dimensions
+    return mp_const_none;
+    /*
     static const mp_arg_t allowed_args[] = {
         { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_INT, {.u_int = 0} },
         { MP_QSTR_M, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } },
@@ -289,43 +204,47 @@ mp_obj_t linalg_eye(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args)
             i++;
         }
     }
-    return MP_OBJ_FROM_PTR(ndarray);
+    return MP_OBJ_FROM_PTR(ndarray); */
 }
 
 mp_obj_t linalg_det(mp_obj_t oin) {
     if(!mp_obj_is_type(oin, &ulab_ndarray_type)) {
         mp_raise_TypeError("function defined for ndarrays only");
     }
-    ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
-    if(in->m != in->n) {
-        mp_raise_ValueError("input must be square matrix");
+    ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);
+    if(ndarray->ndim != 2) {
+        mp_raise_ValueError("only two-dimensional tensors can be inverted");
+    }
+    if(ndarray->shape[0] != ndarray->shape[1]) {
+        mp_raise_ValueError("only square matrices can be inverted");
     }
     
-    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);
+    mp_float_t *tmp = m_new(mp_float_t, ndarray->shape[0]*ndarray->shape[1]);
+    // TODO: this won't work for sliced arrays
+    for(size_t i=0; i < ndarray->len; i++){
+        tmp[i] = ndarray_get_float_value(ndarray->array->items, ndarray->array->typecode, i);
     }
     mp_float_t c;
-    for(size_t m=0; m < in->m-1; m++){
-        if(MICROPY_FLOAT_C_FUN(fabs)(tmp[m*(in->n+1)]) < epsilon) {
-            m_del(mp_float_t, tmp, in->n*in->n);
+    for(size_t m=0; m < ndarray->shape[0]-1; m++){
+        if(MICROPY_FLOAT_C_FUN(fabs)(tmp[m*(ndarray->shape[1]+1)]) < epsilon) {
+            m_del(mp_float_t, tmp, ndarray->shape[0]*ndarray->shape[1]);
             return mp_obj_new_float(0.0);
         }
-        for(size_t n=0; n < in->n; n++){
+        for(size_t n=0; n < ndarray->shape[1]; n++){
             if(m != n) {
-                c = tmp[in->n*n+m] / tmp[m*(in->n+1)];
-                for(size_t k=0; k < in->n; k++){
-                    tmp[in->n*n+k] -= c * tmp[in->n*m+k];
+                c = tmp[ndarray->shape[0]*n+m] / tmp[m*(ndarray->shape[1]+1)];
+                for(size_t k=0; k < ndarray->shape[1]; k++){
+                    tmp[ndarray->shape[1]*n+k] -= c * tmp[ndarray->shape[1]*m+k];
                 }
             }
         }
     }
     mp_float_t det = 1.0;
                             
-    for(size_t m=0; m < in->m; m++){ 
-        det *= tmp[m*(in->n+1)];
+    for(size_t m=0; m < ndarray->shape[0]; m++){ 
+        det *= tmp[m*(ndarray->shape[1]+1)];
     }
-    m_del(mp_float_t, tmp, in->n*in->n);
+    m_del(mp_float_t, tmp, ndarray->shape[0]*ndarray->shape[1]);
     return mp_obj_new_float(det);
 }
 
@@ -333,45 +252,50 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
     if(!mp_obj_is_type(oin, &ulab_ndarray_type)) {
         mp_raise_TypeError("function defined for ndarrays only");
     }
-    ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
-    if(in->m != in->n) {
-        mp_raise_ValueError("input must be square matrix");
+    ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);
+    if(ndarray->ndim != 2) {
+        mp_raise_ValueError("only two-dimensional tensors can be inverted");
     }
-    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);
+    if(ndarray->shape[0] != ndarray->shape[1]) {
+        mp_raise_ValueError("only square matrices can be inverted");
+    }
+    mp_float_t *array = m_new(mp_float_t, ndarray->len);
+    // TODO: this won't work for sliced arrays
+    for(size_t i=0; i < ndarray->len; i++) {
+        array[i] = ndarray_get_float_value(ndarray->array->items, ndarray->array->typecode, i);
     }
     // make sure the matrix is symmetric
-    for(size_t m=0; m < in->m; m++) {
-        for(size_t n=m+1; n < in->n; n++) {
+    for(size_t m=0; m < ndarray->shape[0]; m++) {
+        for(size_t n=m+1; n < ndarray->shape[1]; n++) {
             // compare entry (m, n) to (n, m)
             // TODO: this must probably be scaled!
-            if(epsilon < MICROPY_FLOAT_C_FUN(fabs)(array[m*in->n + n] - array[n*in->n + m])) {
+            if(epsilon < MICROPY_FLOAT_C_FUN(fabs)(array[m*ndarray->shape[0] + n] - array[n*ndarray->shape[0] + m])) {
                 mp_raise_ValueError("input matrix is asymmetric");
             }
         }
     }
     
     // if we got this far, then the matrix will be symmetric
-    
-    ndarray_obj_t *eigenvectors = create_new_ndarray(in->m, in->n, NDARRAY_FLOAT);
+    size_t *shape = m_new(size_t, 2);
+    shape[0] = shape[1] = ndarray->shape[0];
+    ndarray_obj_t *eigenvectors = ndarray_new_dense_ndarray(2, shape, NDARRAY_FLOAT);
     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;
+    for(size_t m=0; m < ndarray->shape[0]; m++) {
+        eigvectors[m*(ndarray->shape[1]+1)] = 1.0;
     }
     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;
+    size_t iterations = JACOBI_MAX*ndarray->shape[0]*ndarray->shape[0];
     do {
         iterations--;
         // find the pivot here
         M = 0;
         N = 0;
         largest = 0.0;
-        for(size_t m=0; m < in->m-1; m++) { // -1: no need to inspect last row
-            for(size_t n=m+1; n < in->n; n++) {
-                w = MICROPY_FLOAT_C_FUN(fabs)(array[m*in->n + n]);
+        for(size_t m=0; m < ndarray->shape[0]-1; m++) { // -1: no need to inspect last row
+            for(size_t n=m+1; n < ndarray->shape[0]; n++) {
+                w = MICROPY_FLOAT_C_FUN(fabs)(array[m*ndarray->shape[0] + n]);
                 if((largest < w) && (epsilon < w)) {
                     M = m;
                     N = n;
@@ -384,7 +308,7 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
         }
         // at this point, we have the pivot, and it is the entry (M, N)
         // now we have to find the rotation angle
-        w = (array[N*in->n + N] - array[M*in->n + M]) / (2.0*array[M*in->n + N]);
+        w = (array[N*ndarray->shape[0] + N] - array[M*ndarray->shape[0] + M]) / (2.0*array[M*ndarray->shape[0] + N]);
         // 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) {
@@ -399,26 +323,26 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
         // at this point, we have the rotation angles, so we can transform the matrix
         // first the two diagonal elements
         // a(M, M) = a(M, M) - t*a(M, N)
-        array[M*in->n + M] = array[M*in->n + M] - t * array[M*in->n + N];
+        array[M*ndarray->shape[0] + M] = array[M*ndarray->shape[0] + M] - t * array[M*ndarray->shape[0] + N];
         // a(N, N) = a(N, N) + t*a(M, N)
-        array[N*in->n + N] = array[N*in->n + N] + t * array[M*in->n + N];
+        array[N*ndarray->shape[0] + N] = array[N*ndarray->shape[0] + N] + t * array[M*ndarray->shape[0] + N];
         // after the rotation, the a(M, N), and a(N, M) entries should become zero
-        array[M*in->n + N] = array[N*in->n + M] = 0.0;
+        array[M*ndarray->shape[0] + N] = array[N*ndarray->shape[0] + M] = 0.0;
         // then all other elements in the column
-        for(size_t k=0; k < in->m; k++) {
+        for(size_t k=0; k < ndarray->shape[0]; k++) {
             if((k == M) || (k == N)) {
                 continue;
             }
-            aMk = array[M*in->n + k];
-            aNk = array[N*in->n + k];
+            aMk = array[M*ndarray->shape[0] + k];
+            aNk = array[N*ndarray->shape[0] + k];
             // a(M, k) = a(M, k) - s*(a(N, k) + tau*a(M, k))
-            array[M*in->n + k] -= s*(aNk + tau*aMk);
+            array[M*ndarray->shape[0] + k] -= s*(aNk + tau*aMk);
             // a(N, k) = a(N, k) + s*(a(M, k) - tau*a(N, k))
-            array[N*in->n + k] += s*(aMk - tau*aNk);
+            array[N*ndarray->shape[0] + k] += s*(aMk - tau*aNk);
             // a(k, M) = a(M, k)
-            array[k*in->n + M] = array[M*in->n + k];
+            array[k*ndarray->shape[0] + M] = array[M*ndarray->shape[0] + k];
             // a(k, N) = a(N, k)
-            array[k*in->n + N] = array[N*in->n + k];
+            array[k*ndarray->shape[0] + N] = array[N*ndarray->shape[0] + k];
         }
         // now we have to update the eigenvectors
         // the rotation matrix, R, multiplies from the right
@@ -427,27 +351,29 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
         // R(N, M) = s
         // (M, N) = -s
         // entries. This means that only the Mth, and Nth columns will change
-        for(size_t m=0; m < in->m; m++) {
-            vm = eigvectors[m*in->n+M];
-            vn = eigvectors[m*in->n+N];
+        for(size_t m=0; m < ndarray->shape[0]; m++) {
+            vm = eigvectors[m*ndarray->shape[0]+M];
+            vn = eigvectors[m*ndarray->shape[0]+N];
             // the new value of eigvectors(m, M)
-            eigvectors[m*in->n+M] = c * vm - s * vn;
+            eigvectors[m*ndarray->shape[0]+M] = c * vm - s * vn;
             // the new value of eigvectors(m, N)
-            eigvectors[m*in->n+N] = s * vm + c * vn;
+            eigvectors[m*ndarray->shape[0]+N] = s * vm + c * vn;
         }
     } while(iterations > 0);
     
     if(iterations == 0) { 
         // the computation did not converge; numpy raises LinAlgError
-        m_del(mp_float_t, array, in->array->len);
+        m_del(mp_float_t, array, ndarray->len);
         mp_raise_ValueError("iterations did not converge");
     }
-    ndarray_obj_t *eigenvalues = create_new_ndarray(1, in->n, NDARRAY_FLOAT);
+    size_t *eigen_shape = m_new(size_t, 1);
+    eigen_shape[0] = ndarray->shape[0];
+    ndarray_obj_t *eigenvalues = ndarray_new_dense_ndarray(1, eigen_shape, NDARRAY_FLOAT);
     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)];
+    for(size_t i=0; i < ndarray->shape[0]; i++) {
+        eigvalues[i] = array[i*(ndarray->shape[0]+1)];
     }
-    m_del(mp_float_t, array, in->array->len);
+    m_del(mp_float_t, array, ndarray->len);
     
     mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(mp_obj_new_tuple(2, NULL));
     tuple->items[0] = MP_OBJ_FROM_PTR(eigenvalues);

code/linalg.h

@@ -13,8 +13,6 @@
 
 #include "ndarray.h"
 
-#define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }
-
 #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
 #define epsilon        1.2e-7
 #elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
@@ -23,9 +21,6 @@
 
 #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(mp_float_t *, size_t );
 mp_obj_t linalg_inv(mp_obj_t );
 mp_obj_t linalg_dot(mp_obj_t , mp_obj_t );

code/micropython.mk

@@ -3,13 +3,11 @@ USERMODULES_DIR := $(USERMOD_DIR)
 
 # Add all C files to SRC_USERMOD.
 SRC_USERMOD += $(USERMODULES_DIR)/ndarray.c
-SRC_USERMOD += $(USERMODULES_DIR)/linalg.c
 SRC_USERMOD += $(USERMODULES_DIR)/vectorise.c
+SRC_USERMOD += $(USERMODULES_DIR)/linalg.c
 SRC_USERMOD += $(USERMODULES_DIR)/poly.c
-SRC_USERMOD += $(USERMODULES_DIR)/fft.c
 SRC_USERMOD += $(USERMODULES_DIR)/numerical.c
+SRC_USERMOD += $(USERMODULES_DIR)/fft.c
 SRC_USERMOD += $(USERMODULES_DIR)/ulab.c
 
-# We can add our module folder to include paths if needed
-# This is not actually needed in this example.
 CFLAGS_USERMOD += -I$(USERMODULES_DIR)

code/ndarray.h

@@ -16,7 +16,10 @@
 #include "py/objstr.h"
 #include "py/objlist.h"
 
-#define PRINT_MAX  10
+#define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }
+
+#define NDARRAY_NUMERIC   0
+#define NDARRAY_BOOLEAN   1
 
 #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
 #define FLOAT_TYPECODE 'f'
@@ -24,9 +27,10 @@
 #define FLOAT_TYPECODE 'd'
 #endif
 
-const mp_obj_type_t ulab_ndarray_type;
+extern const mp_obj_type_t ulab_ndarray_type;
 
 enum NDARRAY_TYPE {
+    NDARRAY_BOOL = '?', // this must never be assigned to the typecode!
     NDARRAY_UINT8 = 'B',
     NDARRAY_INT8 = 'b',
     NDARRAY_UINT16 = 'H', 
@@ -36,21 +40,77 @@ enum NDARRAY_TYPE {
 
 typedef struct _ndarray_obj_t {
     mp_obj_base_t base;
-    size_t m, n;
+    uint8_t boolean;
+    uint8_t ndim;
+    size_t *shape;
+    int32_t *strides;
     size_t len;
+    size_t offset;
     mp_obj_array_t *array;
-    size_t bytes;
 } ndarray_obj_t;
 
-mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t , size_t , mp_obj_iter_buf_t *);
+// this is a helper structure, so that we can return shape AND strides from a function
+typedef struct _ndarray_header_obj_t {
+    size_t *shape;
+    int32_t *strides;
+    int8_t axis;
+    size_t offset;
+} ndarray_header_obj_t;
 
+// various helper functions
+size_t ndarray_index_from_flat(size_t , ndarray_obj_t *, int32_t *);
+size_t ndarray_index_from_contracted(size_t  , ndarray_obj_t * , int32_t * , uint8_t  , uint8_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 );
+// calculates the index (in the original linear array) of an item, if the index in the flat array is given
+// this is the macro equivalent of ndarray_index_from_flat()
+// TODO: This fails, when the last stride is not 1!!!
+#define NDARRAY_INDEX_FROM_FLAT(ndarray, shape_strides, index, _tindex, _nindex) do {\
+    size_t Q;\
+    (_tindex) = (index);\
+    (_nindex) = (ndarray)->offset;\
+    for(size_t _x=0; _x < (ndarray)->ndim; _x++) {\
+        Q = (_tindex) / (shape_strides)[_x];\
+        (_tindex) -= Q * (shape_strides)[_x];\
+        (_nindex) += Q * (ndarray)->strides[_x];\
+    }\
+} while(0)
+
+#define NDARRAY_INDEX_FROM_FLAT2(ndarray, stride_array, shape_strides, index, _tindex, _nindex) do {\
+    size_t Q;\
+    (_tindex) = (index);\
+    (_nindex) = (ndarray)->offset;\
+    for(size_t _x=0; _x < (ndarray)->ndim; _x++) {\
+        Q = (_tindex) / (shape_strides)[_x];\
+        (_tindex) -= Q * (shape_strides)[_x];\
+        (_nindex) += Q * (stride_array)[_x];\
+    }\
+} while(0)
+    
+#define CREATE_SINGLE_ITEM(ndarray, type, typecode, value) do {\
+    (ndarray) = ndarray_new_linear_array(1, (typecode));\
+    type *tmparr = (type *)(ndarray)->array->items;\
+    tmparr[0] = (type)(value);\
+} while(0)
+
+#define RUN_BINARY_LOOP(ndarray, typecode, type_out, type_left, type_right, lhs, rhs, shape, ndim, operator) do {\
+    uint8_t *left = (uint8_t *)(lhs)->array->items;\
+    uint8_t *right = (uint8_t *)(rhs)->array->items;\
+    (ndarray) = ndarray_new_dense_ndarray((ndim), (shape), (typecode));\
+    uint8_t size_left = sizeof(type_left), size_right = sizeof(type_right);\
+    type_out *out = (type_out *)ndarray->array->items;\
+    for(size_t i=0; i < (lhs)->len; i++) {\
+        out[i] = *(type_left *)left + *(type_right *)right;\
+        left += size_left; right += size_right;\
+    }\
+} while(0)
+
+mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t , size_t , mp_obj_iter_buf_t *);
 void ndarray_print(const mp_print_t *, mp_obj_t , mp_print_kind_t );
-void ndarray_assign_elements(mp_obj_array_t *, mp_obj_t , uint8_t , size_t *);
-ndarray_obj_t *create_new_ndarray(size_t , size_t , uint8_t );
+ndarray_obj_t *ndarray_new_ndarray(uint8_t , size_t *, int32_t *, uint8_t );
+ndarray_obj_t *ndarray_new_dense_ndarray(uint8_t , size_t *, uint8_t );
+ndarray_obj_t *ndarray_new_linear_array(size_t , uint8_t );
+ndarray_obj_t *ndarray_copy_view(ndarray_obj_t *, uint8_t );
 
 mp_obj_t ndarray_copy(mp_obj_t );
 mp_obj_t ndarray_make_new(const mp_obj_type_t *, size_t , size_t , const mp_obj_t *);
@@ -60,65 +120,11 @@ mp_obj_t ndarray_binary_op(mp_binary_op_t , mp_obj_t , mp_obj_t );
 mp_obj_t ndarray_unary_op(mp_unary_op_t , mp_obj_t );
 
 mp_obj_t ndarray_shape(mp_obj_t );
-mp_obj_t ndarray_rawsize(mp_obj_t );
+mp_obj_t ndarray_reshape(mp_obj_t , mp_obj_t );
+mp_obj_t ndarray_transpose(mp_obj_t );
 mp_obj_t ndarray_flatten(size_t , const mp_obj_t *, mp_map_t *);
-mp_obj_t ndarray_asbytearray(mp_obj_t );
-
-#define CREATE_SINGLE_ITEM(outarray, type, typecode, value) do {\
-    ndarray_obj_t *tmp = create_new_ndarray(1, 1, (typecode));\
-    type *tmparr = (type *)tmp->array->items;\
-    tmparr[0] = (type)(value);\
-    (outarray) = MP_OBJ_FROM_PTR(tmp);\
-} while(0)
-
-/*  
-    mp_obj_t row = mp_obj_new_list(n, NULL);
-    mp_obj_list_t *row_ptr = MP_OBJ_TO_PTR(row);
-    
-    should work outside the loop, but it doesn't. Go figure! 
-*/
-
-#define RUN_BINARY_LOOP(typecode, type_out, type_left, type_right, ol, or, op) do {\
-    type_left *left = (type_left *)(ol)->array->items;\
-    type_right *right = (type_right *)(or)->array->items;\
-    uint8_t inc = 0;\
-    if((or)->array->len > 1) inc = 1;\
-    if(((op) == MP_BINARY_OP_ADD) || ((op) == MP_BINARY_OP_SUBTRACT) || ((op) == MP_BINARY_OP_MULTIPLY)) {\
-        ndarray_obj_t *out = create_new_ndarray(ol->m, ol->n, typecode);\
-        type_out *(odata) = (type_out *)out->array->items;\
-        if((op) == MP_BINARY_OP_ADD) { for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = left[i] + right[j];}\
-        if((op) == MP_BINARY_OP_SUBTRACT) { for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = left[i] - right[j];}\
-        if((op) == MP_BINARY_OP_MULTIPLY) { for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = left[i] * right[j];}\
-        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);\
-        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)) {\
-        mp_obj_t out_list = mp_obj_new_list(0, NULL);\
-        size_t m = (ol)->m, n = (ol)->n;\
-        for(size_t i=0, r=0; i < m; i++, r+=inc) {\
-            mp_obj_t row = mp_obj_new_list(n, NULL);\
-            mp_obj_list_t *row_ptr = MP_OBJ_TO_PTR(row);\
-            for(size_t j=0, s=0; j < n; j++, s+=inc) {\
-                row_ptr->items[j] = mp_const_false;\
-                if((op) == MP_BINARY_OP_LESS) {\
-                    if(left[i*n+j] < right[r*n+s]) row_ptr->items[j] = mp_const_true;\
-                } else if((op) == MP_BINARY_OP_LESS_EQUAL) {\
-                    if(left[i*n+j] <= right[r*n+s]) row_ptr->items[j] = mp_const_true;\
-                } else if((op) == MP_BINARY_OP_MORE) {\
-                    if(left[i*n+j] > right[r*n+s]) row_ptr->items[j] = mp_const_true;\
-                } else if((op) == MP_BINARY_OP_MORE_EQUAL) {\
-                    if(left[i*n+j] >= right[r*n+s]) row_ptr->items[j] = mp_const_true;\
-                }\
-            }\
-            if(m == 1) return row;\
-            mp_obj_list_append(out_list, row);\
-        }\
-        return out_list;\
-    }\
-} while(0)
+mp_obj_t ndarray_itemsize(mp_obj_t );
+mp_obj_t ndarray_strides(mp_obj_t );
+mp_obj_t ndarray_info(mp_obj_t );
 
 #endif

code/numerical.c

@@ -28,6 +28,32 @@ enum NUMERICAL_FUNCTION_TYPE {
     NUMERICAL_STD,
 };
 
+// creates the shape and strides arrays of the contracted ndarray
+ndarray_header_obj_t contracted_shape_strides(ndarray_obj_t *ndarray, int8_t axis) {
+    if(axis < 0) axis += ndarray->ndim;
+    if((axis > ndarray->ndim-1) || (axis < 0)) {
+        mp_raise_ValueError("tuple index out of range");
+    }
+    size_t *shape = m_new(size_t, ndarray->ndim-1);
+    int32_t *strides = m_new(int32_t, ndarray->ndim-1);
+    for(size_t i=0, j=0; i < ndarray->ndim; i++) {
+        if(axis != i) {
+            shape[j] = ndarray->shape[j];
+            j++;
+        }
+    }
+    int32_t stride = 1;
+    for(size_t i=0; i < ndarray->ndim-1; i++) {
+        strides[ndarray->ndim-2-i] = stride;
+        stride *= shape[ndarray->ndim-2-i];
+    }
+    ndarray_header_obj_t header;
+    header.shape = shape;
+    header.strides = strides;
+    header.axis = axis;
+    return header;
+}
+
 mp_obj_t numerical_linspace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
     static const mp_arg_t allowed_args[] = {
         { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } },
@@ -50,7 +76,7 @@ mp_obj_t numerical_linspace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *k
     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);
     else step = (mp_obj_get_float(args[1].u_obj)-value)/len;
-    ndarray_obj_t *ndarray = create_new_ndarray(1, len, typecode);
+    ndarray_obj_t *ndarray = ndarray_new_linear_array(len, typecode);
     if(typecode == NDARRAY_UINT8) {
         uint8_t *array = (uint8_t *)ndarray->array->items;
         for(size_t i=0; i < len; i++, value += step) array[i] = (uint8_t)value;
@@ -77,120 +103,71 @@ mp_obj_t numerical_linspace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *k
     }
 }
 
-mp_obj_t numerical_sum_mean_std_array(mp_obj_t oin, uint8_t optype) {
-    mp_float_t value, sum = 0.0, sq_sum = 0.0;
-    mp_obj_iter_buf_t iter_buf;
-    mp_obj_t item, iterable = mp_getiter(oin, &iter_buf);
-    mp_int_t len = mp_obj_get_int(mp_obj_len(oin));
-    while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
-        value = mp_obj_get_float(item);
-        sum += value;
-        if(optype == NUMERICAL_STD) {
-            sq_sum += value*value;
+mp_obj_t numerical_flat_sum_mean_std(ndarray_obj_t *ndarray, uint8_t optype, size_t ddof) {
+    mp_float_t value;
+    int32_t *shape_strides = m_new(int32_t, ndarray->ndim);
+    shape_strides[ndarray->ndim-1] = ndarray->strides[ndarray->ndim-1];
+    for(uint8_t i=ndarray->ndim-1; i > 0; i--) {
+        shape_strides[i-1] = shape_strides[i] * ndarray->shape[i-1];
     }
+    if(ndarray->array->typecode == NDARRAY_UINT8) {
+        CALCULATE_FLAT_SUM_STD(ndarray, uint8_t, value, shape_strides, optype);
+    } else if(ndarray->array->typecode == NDARRAY_INT8) {
+        CALCULATE_FLAT_SUM_STD(ndarray, int8_t, value, shape_strides, optype);
+    } if(ndarray->array->typecode == NDARRAY_UINT16) {
+        CALCULATE_FLAT_SUM_STD(ndarray, uint16_t, value, shape_strides, optype);
+    } else if(ndarray->array->typecode == NDARRAY_INT16) {
+        CALCULATE_FLAT_SUM_STD(ndarray, int16_t, value, shape_strides, optype);
+    } else {
+        CALCULATE_FLAT_SUM_STD(ndarray, mp_float_t, value, shape_strides, optype);
     }
+    m_del(int32_t, shape_strides, ndarray->ndim);
     if(optype == NUMERICAL_SUM) {
-        return mp_obj_new_float(sum);
+        return mp_obj_new_float(value);
     } else if(optype == NUMERICAL_MEAN) {
-        return mp_obj_new_float(sum/len);
+        return mp_obj_new_float(value/ndarray->len);
     } else {
-        sum /= len; // this is now the mean!
-        return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/len-sum*sum));
+        return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(value/(ndarray->len-ddof)));
     }
 }   
 
-STATIC mp_float_t numerical_sum_mean_std_single_line(void *data, size_t start, size_t stop, 
-                                                  size_t stride, uint8_t typecode, uint8_t optype) {
-    
-    mp_float_t sum = 0.0, sq_sum = 0.0, value;
-    size_t len = 0;
-    for(size_t i=start; i < stop; i+=stride, len++) {
-        value = ndarray_get_float_value(data, typecode, i);        
-        sum += value;
-        if(optype == NUMERICAL_STD) {
-            sq_sum += value*value;
-        }
-    }
-    if(len == 0) {
-        mp_raise_ValueError("data length is 0!");
-    }
-    if(optype ==  NUMERICAL_SUM) {
-        return sum;
-    } else if(optype == NUMERICAL_MEAN) {
-        return sum/len;
+// numerical functions for ndarrays
+mp_obj_t numerical_sum_mean_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {
+    if(axis == mp_const_none) {
+        return numerical_flat_sum_mean_std(ndarray, optype, 0);
     } else {
-        sum /= len; // this is now the mean!
-        return MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/len-sum*sum);
-    }
-}
-
-STATIC mp_obj_t numerical_sum_mean_std_matrix(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {
-    ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
-    if((axis == mp_const_none) || (in->m == 1) || (in->n == 1)) { 
-        // return the value for the flattened array
-        return mp_obj_new_float(numerical_sum_mean_std_single_line(in->array->items, 0, 
-                                                      in->array->len, 1, in->array->typecode, optype));
+        int8_t ax = mp_obj_get_int(axis);
+        ndarray_header_obj_t header = contracted_shape_strides(ndarray, ax);
+        ndarray_obj_t *result = ndarray_new_ndarray(ndarray->ndim-1, header.shape, header.strides, NDARRAY_FLOAT);
+        mp_float_t *farray = (mp_float_t *)result->array->items;
+        size_t offset;
+        // iterate along the length of the output array, so as to avoid recursion
+        for(size_t i=0; i < result->array->len; i++) {
+            offset = ndarray_index_from_contracted(i, ndarray, result->strides, result->ndim, header.axis) + ndarray->offset;
+            if(ndarray->array->typecode == NDARRAY_UINT8) {
+                CALCULATE_SUM(ndarray, uint8_t, farray, ndarray->shape[header.axis], i, ndarray->strides[header.axis], offset, optype);
+            } else if(ndarray->array->typecode == NDARRAY_INT8) {
+                CALCULATE_SUM(ndarray, int8_t, farray, ndarray->shape[header.axis], i, ndarray->strides[header.axis], offset, optype); 
+            } else if(ndarray->array->typecode == NDARRAY_UINT16) {
+                CALCULATE_SUM(ndarray, uint16_t, farray, ndarray->shape[header.axis], i, ndarray->strides[header.axis], offset, optype);     
+            } else if(ndarray->array->typecode == NDARRAY_INT16) {
+                CALCULATE_SUM(ndarray, int16_t, farray, ndarray->shape[header.axis], i, ndarray->strides[header.axis], offset, optype);
             } else {
-        uint8_t _axis = mp_obj_get_int(axis);
-        size_t m = (_axis == 0) ? 1 : in->m;
-        size_t n = (_axis == 0) ? in->n : 1;
-        size_t len = in->array->len;
-        mp_float_t sms;
-        // TODO: pass in->array->typcode to create_new_ndarray
-        ndarray_obj_t *out = create_new_ndarray(m, n, NDARRAY_FLOAT);
-
-        // TODO: these two cases could probably be combined in a more elegant fashion...
-        if(_axis == 0) { // vertical
-            for(size_t i=0; i < n; i++) {
-                sms = numerical_sum_mean_std_single_line(in->array->items, i, len, 
-                                                               n, in->array->typecode, optype);
-                ((float_t *)out->array->items)[i] = sms;
+                CALCULATE_SUM(ndarray, mp_float_t, farray, ndarray->shape[header.axis], i, ndarray->strides[header.axis], offset, optype);       
             }
-        } else { // horizontal
-            for(size_t i=0; i < m; i++) {
-                sms = numerical_sum_mean_std_single_line(in->array->items, i*in->n, 
-                                                               (i+1)*in->n, 1, in->array->typecode, optype);
-                ((float_t *)out->array->items)[i] = sms;
+            if(optype == NUMERICAL_MEAN) farray[i] /= ndarray->shape[header.axis];
         }
+        return MP_OBJ_FROM_PTR(result);
     }
-    return MP_OBJ_FROM_PTR(out);
-    }
+    return mp_const_none;
 }
 
-size_t numerical_argmin_argmax_array(ndarray_obj_t *in, size_t start, 
-                                       size_t stop, size_t stride, uint8_t op) {
-    size_t best_idx = start;
-    if(in->array->typecode == NDARRAY_UINT8) {
-        ARG_MIN_LOOP(in, uint8_t, start, stop, stride, op);
-    } else if(in->array->typecode == NDARRAY_INT8) {
-        ARG_MIN_LOOP(in, int8_t, start, stop, stride, op);
-    } else if(in->array->typecode == NDARRAY_UINT16) {
-        ARG_MIN_LOOP(in, uint16_t, start, stop, stride, op);
-    } else if(in->array->typecode == NDARRAY_INT16) {
-        ARG_MIN_LOOP(in, uint16_t, start, stop, stride, op);
-    } else if(in->array->typecode == NDARRAY_FLOAT) {
-        ARG_MIN_LOOP(in, mp_float_t, start, stop, stride, op);
-    }
-    return best_idx;
+mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {
+    return mp_const_none;
 }
 
-void copy_value_into_ndarray(ndarray_obj_t *target, ndarray_obj_t *source, size_t target_idx, size_t source_idx) {
-    // since we are simply copying, it doesn't matter, whether the arrays are signed or unsigned, 
-    // we can cast them in any way we like
-    // This could also be done with byte copies. I don't know, whether that would have any benefits
-    if((target->array->typecode == NDARRAY_UINT8) || (target->array->typecode == NDARRAY_INT8)) {
-        ((uint8_t *)target->array->items)[target_idx] = ((uint8_t *)source->array->items)[source_idx];
-    } else if((target->array->typecode == NDARRAY_UINT16) || (target->array->typecode == NDARRAY_INT16)) {
-        ((uint16_t *)target->array->items)[target_idx] = ((uint16_t *)source->array->items)[source_idx];
-    } else { 
-        ((float *)target->array->items)[target_idx] = ((float *)source->array->items)[source_idx];
-    }
-}
- 
-STATIC mp_obj_t numerical_argmin_argmax(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {
-    if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) || 
-        MP_OBJ_IS_TYPE(oin, &mp_type_range)) {
-        // This case will work for single iterables only 
+// numerical function for interables (single axis)
+mp_obj_t numerical_argmin_argmax_iterable(mp_obj_t oin, uint8_t optype) {
     size_t idx = 0, best_idx = 0;
     mp_obj_iter_buf_t iter_buf;
     mp_obj_t iterable = mp_getiter(oin, &iter_buf);
@@ -210,105 +187,71 @@ STATIC mp_obj_t numerical_argmin_argmax(mp_obj_t oin, mp_obj_t axis, uint8_t opt
     } else {
         return best_obj;
     }    
-    } else if(mp_obj_is_type(oin, &ulab_ndarray_type)) {
-            ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
-            size_t best_idx;
-            if((axis == mp_const_none) || (in->m == 1) || (in->n == 1)) {
-                // return the value for the flattened array                
-                best_idx = numerical_argmin_argmax_array(in, 0, in->array->len, 1, optype);
-                if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {
-                    return MP_OBJ_NEW_SMALL_INT(best_idx);
-                } else {
-                    if(in->array->typecode == NDARRAY_FLOAT) {
-                        return mp_obj_new_float(ndarray_get_float_value(in->array->items, in->array->typecode, best_idx));
-                    } else {
-                        return mp_binary_get_val_array(in->array->typecode, in->array->items, best_idx);
-                    }
-                }
-            } else { // we have to work with a full matrix here
-                uint8_t _axis = mp_obj_get_int(axis);
-                size_t m = (_axis == 0) ? 1 : in->m;
-                size_t n = (_axis == 0) ? in->n : 1;
-                size_t len = in->array->len;
-                ndarray_obj_t *ndarray = NULL;
-                if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {
-                    ndarray = create_new_ndarray(m, n, in->array->typecode);
-                } else { // argmin/argmax
-                    // TODO: one might get away with uint8_t, if both m, and n < 255
-                    ndarray = create_new_ndarray(m, n, NDARRAY_UINT16);
 }
 
-                // TODO: these two cases could probably be combined in a more elegant fashion...
-                if(_axis == 0) { // vertical
-                    for(size_t i=0; i < n; i++) {
-                        best_idx = numerical_argmin_argmax_array(in, i, len, n, optype);
-                        if((optype == NUMERICAL_MIN) || (optype == NUMERICAL_MAX)) {
-                            copy_value_into_ndarray(ndarray, in, i, best_idx);
-                        } else {
-                            ((uint16_t *)ndarray->array->items)[i] = (uint16_t)(best_idx / n);
+mp_obj_t numerical_sum_mean_std_iterable(mp_obj_t oin, uint8_t optype, size_t ddof) {
+    mp_float_t value, sum = 0.0, sq_sum = 0.0;
+    mp_obj_iter_buf_t iter_buf;
+    mp_obj_t item, iterable = mp_getiter(oin, &iter_buf);
+    mp_int_t len = mp_obj_get_int(mp_obj_len(oin));
+    while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
+        value = mp_obj_get_float(item);
+        sum += value;
     }
+    if(optype ==  NUMERICAL_SUM) {
+        return mp_obj_new_float(sum);
+    } else if(optype == NUMERICAL_MEAN) {
+        return mp_obj_new_float(sum/len);
+    } else { // this should be the case of the standard deviation
+        sum /= len; // this is the mean now
+        iterable = mp_getiter(oin, &iter_buf);
+        while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
+            value = mp_obj_get_float(item) - sum;
+            sq_sum += value * value;
         }
-                } else { // horizontal
-                    for(size_t i=0; i < m; i++) {
-                        best_idx = numerical_argmin_argmax_array(in, i*in->n, (i+1)*in->n, 1, optype);
-                        if((optype == NUMERICAL_MIN) || (optype == NUMERICAL_MAX)) {
-                             copy_value_into_ndarray(ndarray, in, i, best_idx);
-                        } else {
-                            ((uint16_t *)ndarray->array->items)[i] = (uint16_t)(best_idx - i*in->n);
+        return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/(len-ddof)));
     }
 }
-                }
-                return MP_OBJ_FROM_PTR(ndarray);
-            }
-            return mp_const_none;
-        }
-    mp_raise_TypeError("input type is not supported");
-}
 
-STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t type) {
+STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t optype) {
     static const mp_arg_t allowed_args[] = {
         { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} } ,
-        { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },
+        { MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },
     };
 
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
     
     mp_obj_t oin = args[0].u_obj;
     mp_obj_t axis = args[1].u_obj;
-    if((axis != mp_const_none) && (mp_obj_get_int(axis) != 0) && (mp_obj_get_int(axis) != 1)) {
-        // this seems to pass with False, and True...
-        mp_raise_ValueError("axis must be None, 0, or 1");
-    }
     
     if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) || 
         MP_OBJ_IS_TYPE(oin, &mp_type_range)) {
-        switch(type) {
+        switch(optype) {
             case NUMERICAL_MIN:
             case NUMERICAL_ARGMIN:
             case NUMERICAL_MAX:
             case NUMERICAL_ARGMAX:
-                return numerical_argmin_argmax(oin, axis, type);
+                return numerical_argmin_argmax_iterable(oin, optype);
             case NUMERICAL_SUM:
             case NUMERICAL_MEAN:
-            case NUMERICAL_STD:
-                return numerical_sum_mean_std_array(oin, type);
-            default: // we should never reach this point, but whatever
+                return numerical_sum_mean_std_iterable(oin, optype, 0);
+            default: // we should never end up here
                 return mp_const_none;
         }
     } else if(MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {
-        switch(type) {
+        ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);
+        switch(optype) {
             case NUMERICAL_MIN:
-            case NUMERICAL_MAX:
             case NUMERICAL_ARGMIN:
+            case NUMERICAL_MAX:
             case NUMERICAL_ARGMAX:
-                return numerical_argmin_argmax(oin, axis, type);
+                return numerical_argmin_argmax_ndarray(ndarray, axis, optype);
             case NUMERICAL_SUM:
             case NUMERICAL_MEAN:
-            case NUMERICAL_STD:
-                return numerical_sum_mean_std_matrix(oin, axis, type);            
+                return numerical_sum_mean_ndarray(ndarray, args[1].u_obj, optype);
             default:
-                mp_raise_NotImplementedError("operation is not implemented on ndarrays");
+                return mp_const_none;
         }
     } else {
         mp_raise_TypeError("input must be tuple, list, range, or ndarray");
@@ -316,22 +259,6 @@ STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_m
     return mp_const_none;
 }
 
-mp_obj_t numerical_min(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
-    return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MIN);
-}
-
-mp_obj_t numerical_max(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
-    return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MAX);
-}
-
-mp_obj_t numerical_argmin(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
-    return numerical_function(n_args, pos_args, kw_args, NUMERICAL_ARGMIN);
-}
-
-mp_obj_t numerical_argmax(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
-    return numerical_function(n_args, pos_args, kw_args, NUMERICAL_ARGMAX);
-}
-
 mp_obj_t numerical_sum(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
     return numerical_function(n_args, pos_args, kw_args, NUMERICAL_SUM);
 }
@@ -341,351 +268,49 @@ mp_obj_t numerical_mean(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
 }
 
 mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
-    return numerical_function(n_args, pos_args, kw_args, NUMERICAL_STD);
-}
-
-mp_obj_t numerical_roll(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_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } },
-        { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },
+        { MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },
+        { MP_QSTR_ddof, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} },
     };
 
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(2, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
     
     mp_obj_t oin = args[0].u_obj;
-    int16_t shift = mp_obj_get_int(args[1].u_obj);
-    if((args[2].u_obj != mp_const_none) && 
-           (mp_obj_get_int(args[2].u_obj) != 0) && 
-           (mp_obj_get_int(args[2].u_obj) != 1)) {
-        mp_raise_ValueError("axis must be None, 0, or 1");
-    }
-
-    ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
-    uint8_t _sizeof = mp_binary_get_size('@', in->array->typecode, NULL);
-    size_t len;
-    int16_t _shift;
-    uint8_t *array = (uint8_t *)in->array->items;
-    // TODO: transpose the matrix, if axis == 0. Though, that is hard on the RAM...
-    if(shift < 0) {
-        _shift = -shift;
+    mp_obj_t axis = args[1].u_obj;
+    size_t ddof = args[2].u_int;
+    if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) || MP_OBJ_IS_TYPE(oin, &mp_type_range)) {
+        return numerical_sum_mean_std_iterable(oin, NUMERICAL_STD, ddof);
+    } else if(MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {
+        ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);
+        if(axis == mp_const_none) { // calculate for the flat array
+            return numerical_flat_sum_mean_std(ndarray, NUMERICAL_STD, ddof);
         } else {
-        _shift = shift;
-    }
-    if((args[2].u_obj == mp_const_none) || (mp_obj_get_int(args[2].u_obj) == 1)) { // shift horizontally
-        uint16_t M;
-        if(args[2].u_obj == mp_const_none) {
-            len = in->array->len;
-            M = 1;
+            int8_t ax = mp_obj_get_int(axis);
+            ndarray_header_obj_t header = contracted_shape_strides(ndarray, ax);
+            ndarray_obj_t *result = ndarray_new_ndarray(ndarray->ndim-1, header.shape, header.strides, NDARRAY_FLOAT);
+            mp_float_t *farray = (mp_float_t *)result->array->items, sum_sq;
+            size_t offset;
+            // iterate along the length of the output array, so as to avoid recursion
+            for(size_t i=0; i < result->array->len; i++) {
+                offset = ndarray_index_from_contracted(i, ndarray, result->strides, result->ndim, header.axis) + ndarray->offset;
+                if(ndarray->array->typecode == NDARRAY_UINT8) {
+                    CALCULATE_STD(ndarray, uint8_t, sum_sq, ndarray->shape[header.axis], ndarray->strides[header.axis], offset);
+                } else if(ndarray->array->typecode == NDARRAY_INT8) {
+                    CALCULATE_STD(ndarray, int8_t, sum_sq, ndarray->shape[header.axis], ndarray->strides[header.axis], offset); 
+                } else if(ndarray->array->typecode == NDARRAY_UINT16) {
+                    CALCULATE_STD(ndarray, uint16_t, sum_sq, ndarray->shape[header.axis], ndarray->strides[header.axis], offset);     
+                } else if(ndarray->array->typecode == NDARRAY_INT16) {
+                    CALCULATE_STD(ndarray, int16_t, sum_sq, ndarray->shape[header.axis], ndarray->strides[header.axis], offset);
                 } else {
-            len = in->n;
-            M = in->m;
+                    CALCULATE_STD(ndarray, mp_float_t, sum_sq, ndarray->shape[header.axis], ndarray->strides[header.axis], offset);
                 }
-        _shift = _shift % len;
-        if(shift < 0) _shift = len - _shift;
-        // TODO: if(shift > len/2), we should move in the opposite direction. That would save RAM
-        _shift *= _sizeof;
-        uint8_t *tmp = m_new(uint8_t, _shift);
-        for(size_t m=0; m < M; m++) {
-            memmove(tmp, &array[m*len*_sizeof], _shift);
-            memmove(&array[m*len*_sizeof], &array[m*len*_sizeof+_shift], len*_sizeof-_shift);
-            memmove(&array[(m+1)*len*_sizeof-_shift], tmp, _shift);
+                farray[i] = MICROPY_FLOAT_C_FUN(sqrt)(sum_sq / (ndarray->shape[header.axis] - ddof));
             }
-        m_del(uint8_t, tmp, _shift);
-        return mp_const_none;
-    } else {
-        len = in->m;
-        // temporary buffer
-        uint8_t *_data = m_new(uint8_t, _sizeof*len);
-        
-        _shift = _shift % len;
-        if(shift < 0) _shift = len - _shift;
-        _shift *= _sizeof;
-        uint8_t *tmp = m_new(uint8_t, _shift);
-
-        for(size_t n=0; n < in->n; n++) {
-            for(size_t m=0; m < len; m++) {
-                // this loop should fill up the temporary buffer
-                memmove(&_data[m*_sizeof], &array[(m*in->n+n)*_sizeof], _sizeof);
+            return MP_OBJ_FROM_PTR(result);
         }
-            // now, the actual shift
-            memmove(tmp, _data, _shift);
-            memmove(_data, &_data[_shift], len*_sizeof-_shift);
-            memmove(&_data[len*_sizeof-_shift], tmp, _shift);
-            for(size_t m=0; m < len; m++) {
-                // this loop should dump the content of the temporary buffer into data
-                memmove(&array[(m*in->n+n)*_sizeof], &_data[m*_sizeof], _sizeof);
+        mp_raise_TypeError("input must be tuple, list, range, or ndarray");
     }
-        }
-        m_del(uint8_t, tmp, _shift);
-        m_del(uint8_t, _data, _sizeof*len);
     return mp_const_none;
 }
-}
-
-mp_obj_t numerical_flip(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_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} },
-    };
-
-    mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
-    
-    if(!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type)) {
-        mp_raise_TypeError("flip argument must be an ndarray");
-    }
-    if((args[1].u_obj != mp_const_none) && 
-           (mp_obj_get_int(args[1].u_obj) != 0) && 
-           (mp_obj_get_int(args[1].u_obj) != 1)) {
-        mp_raise_ValueError("axis must be None, 0, or 1");
-    }
-
-    ndarray_obj_t *in = MP_OBJ_TO_PTR(args[0].u_obj);
-    mp_obj_t oout = ndarray_copy(args[0].u_obj);
-    ndarray_obj_t *out = MP_OBJ_TO_PTR(oout);
-    uint8_t _sizeof = mp_binary_get_size('@', in->array->typecode, NULL);
-    uint8_t *array_in = (uint8_t *)in->array->items;
-    uint8_t *array_out = (uint8_t *)out->array->items;    
-    size_t len;
-    if((args[1].u_obj == mp_const_none) || (mp_obj_get_int(args[1].u_obj) == 1)) { // flip horizontally
-        uint16_t M = in->m;
-        len = in->n;
-        if(args[1].u_obj == mp_const_none) { // flip flattened array
-            len = in->array->len;
-            M = 1;
-        }
-        for(size_t m=0; m < M; m++) {
-            for(size_t n=0; n < len; n++) {
-                memcpy(array_out+_sizeof*(m*len+n), array_in+_sizeof*((m+1)*len-n-1), _sizeof);
-            }
-        }
-    } else { // flip vertically
-        for(size_t m=0; m < in->m; m++) {
-            for(size_t n=0; n < in->n; n++) {
-                memcpy(array_out+_sizeof*(m*in->n+n), array_in+_sizeof*((in->m-m-1)*in->n+n), _sizeof);
-            }
-        }
-    }
-    return out;
-}
-
-mp_obj_t numerical_diff(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("diff argument must be an ndarray");
-    }
-    
-    ndarray_obj_t *in = MP_OBJ_TO_PTR(args[0].u_obj);
-    size_t increment, N, M;
-    if((args[2].u_int == -1) || (args[2].u_int == 1)) { // differentiate along the horizontal axis
-        increment = 1;
-    } else if(args[2].u_int == 0) { // differtiate along vertical axis
-        increment = in->n;
-    } else {
-        mp_raise_ValueError("axis must be -1, 0, or 1");        
-    }
-    if((args[1].u_int < 0) || (args[1].u_int > 9)) {
-        mp_raise_ValueError("n must be between 0, and 9");
-    }
-    uint8_t n = args[1].u_int;
-    int8_t *stencil = m_new(int8_t, n+1);
-    stencil[0] = 1;
-    for(uint8_t i=1; i < n+1; i++) {
-        stencil[i] = -stencil[i-1]*(n-i+1)/i;
-    }
-
-    ndarray_obj_t *out;
-    
-    if(increment == 1) { // differentiate along the horizontal axis 
-        if(n >= in->n) {
-            out = create_new_ndarray(in->m, 0, in->array->typecode);
-            m_del(uint8_t, stencil, n);
-            return MP_OBJ_FROM_PTR(out);
-        }
-        N = in->n - n;
-        M = in->m;
-    } else { // differentiate along vertical axis
-        if(n >= in->m) {
-            out = create_new_ndarray(0, in->n, in->array->typecode);
-            m_del(uint8_t, stencil, n);
-            return MP_OBJ_FROM_PTR(out);
-        }
-        M = in->m - n;
-        N = in->n;
-    }
-    out = create_new_ndarray(M, N, in->array->typecode);
-    if(in->array->typecode == NDARRAY_UINT8) {
-        CALCULATE_DIFF(in, out, uint8_t, M, N, in->n, increment);
-    } else if(in->array->typecode == NDARRAY_INT8) {
-        CALCULATE_DIFF(in, out, int8_t, M, N, in->n, increment);
-    }  else if(in->array->typecode == NDARRAY_UINT16) {
-        CALCULATE_DIFF(in, out, uint16_t, M, N, in->n, increment);
-    } else if(in->array->typecode == NDARRAY_INT16) {
-        CALCULATE_DIFF(in, out, int16_t, M, N, in->n, increment);
-    } else {
-        CALCULATE_DIFF(in, out, mp_float_t, M, N, in->n, increment);
-    }
-    m_del(int8_t, stencil, n);
-    return MP_OBJ_FROM_PTR(out);
-}
-
-mp_obj_t numerical_sort_helper(mp_obj_t oin, mp_obj_t axis, uint8_t inplace) {
-    if(!mp_obj_is_type(oin, &ulab_ndarray_type)) {
-        mp_raise_TypeError("sort argument must be an ndarray");
-    }
-
-    ndarray_obj_t *ndarray;
-    mp_obj_t out;
-    if(inplace == 1) {
-        ndarray = MP_OBJ_TO_PTR(oin);
-    } else {
-        out = ndarray_copy(oin);
-        ndarray = MP_OBJ_TO_PTR(out);
-    }
-    size_t increment, start_inc, end, N;
-    if(axis == mp_const_none) { // flatten the array
-        ndarray->m = 1;
-        ndarray->n = ndarray->array->len;
-        increment = 1;
-        start_inc = ndarray->n;
-        end = ndarray->n;
-        N = ndarray->n;
-    } else if((mp_obj_get_int(axis) == -1) || 
-              (mp_obj_get_int(axis) == 1)) { // sort along the horizontal axis
-        increment = 1;
-        start_inc = ndarray->n;
-        end = ndarray->array->len;
-        N = ndarray->n;
-    } else if(mp_obj_get_int(axis) == 0) { // sort along vertical axis
-        increment = ndarray->n;
-        start_inc = 1;
-        end = ndarray->m;
-        N = ndarray->m;
-    } else {
-        mp_raise_ValueError("axis must be -1, 0, None, or 1");        
-    }
-    
-    size_t q, k, p, c;
-
-    for(size_t start=0; start < end; start+=start_inc) {
-        q = N; 
-        k = (q >> 1);
-        if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {
-            HEAPSORT(uint8_t, ndarray);
-        } else if((ndarray->array->typecode == NDARRAY_INT16) || (ndarray->array->typecode == NDARRAY_INT16)) {
-            HEAPSORT(uint16_t, ndarray);
-        } else {
-            HEAPSORT(mp_float_t, ndarray);
-        }
-    }
-    if(inplace == 1) {
-        return mp_const_none;
-    } else {
-        return out;
-    }
-}
-
-// numpy function
-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_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.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);
-
-    return numerical_sort_helper(args[0].u_obj, args[1].u_obj, 0);
-}
-// method of an ndarray
-mp_obj_t numerical_sort_inplace(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_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.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);
-
-    return numerical_sort_helper(args[0].u_obj, args[1].u_obj, 1);
-}
-
-mp_obj_t numerical_argsort(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_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.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("argsort argument must be an ndarray");
-    }
-
-    ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
-    size_t increment, start_inc, end, N, m, n;
-    if(args[1].u_obj == mp_const_none) { // flatten the array
-        m = 1;
-        n = ndarray->array->len;
-        ndarray->m = m;
-        ndarray->n = n;
-        increment = 1;
-        start_inc = ndarray->n;
-        end = ndarray->n;
-        N = n;
-    } else if((mp_obj_get_int(args[1].u_obj) == -1) || 
-              (mp_obj_get_int(args[1].u_obj) == 1)) { // sort along the horizontal axis
-        m = ndarray->m;
-        n = ndarray->n;
-        increment = 1;
-        start_inc = n;
-        end = ndarray->array->len;
-        N = n;
-    } else if(mp_obj_get_int(args[1].u_obj) == 0) { // sort along vertical axis
-        m = ndarray->m;
-        n = ndarray->n;
-        increment = n;
-        start_inc = 1;
-        end = m;
-        N = m;
-    } else {
-        mp_raise_ValueError("axis must be -1, 0, None, or 1");
-    }
-
-    // at the expense of flash, we could save RAM by creating 
-    // an NDARRAY_UINT16 ndarray only, if needed, otherwise, NDARRAY_UINT8
-    ndarray_obj_t *indices = create_new_ndarray(m, n, NDARRAY_UINT16);
-    uint16_t *index_array = (uint16_t *)indices->array->items;
-    // initialise the index array
-    // if array is flat: 0 to indices->n
-    // if sorting vertically, identical indices are arranged row-wise
-    // if sorting horizontally, identical indices are arranged colunn-wise
-    for(uint16_t start=0; start < end; start+=start_inc) {
-        for(uint16_t s=0; s < N; s++) {
-            index_array[start+s*increment] = s;
-        }
-    }
-
-    size_t q, k, p, c;
-    for(size_t start=0; start < end; start+=start_inc) {
-        q = N; 
-        k = (q >> 1);
-        if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {
-            HEAP_ARGSORT(uint8_t, ndarray, index_array);
-        } else if((ndarray->array->typecode == NDARRAY_INT16) || (ndarray->array->typecode == NDARRAY_INT16)) {
-            HEAP_ARGSORT(uint16_t, ndarray, index_array);
-        } else {
-            HEAP_ARGSORT(mp_float_t, ndarray, index_array);
-        }
-    }
-    return MP_OBJ_FROM_PTR(indices);
-}

code/numerical.h

@@ -14,121 +14,52 @@
 #include "ndarray.h"
 
 mp_obj_t numerical_linspace(size_t , const mp_obj_t *, mp_map_t *);
+
 mp_obj_t numerical_sum(size_t , const mp_obj_t *, mp_map_t *);
 mp_obj_t numerical_mean(size_t , const mp_obj_t *, mp_map_t *);
 mp_obj_t numerical_std(size_t , const mp_obj_t *, mp_map_t *);
-mp_obj_t numerical_min(size_t , const mp_obj_t *, mp_map_t *);
-mp_obj_t numerical_max(size_t , const mp_obj_t *, mp_map_t *);
-mp_obj_t numerical_argmin(size_t , const mp_obj_t *, mp_map_t *);
-mp_obj_t numerical_argmax(size_t , const mp_obj_t *, mp_map_t *);
-mp_obj_t numerical_roll(size_t , const mp_obj_t *, mp_map_t *);
 
-// TODO: implement minimum/maximum, and cumsum
-mp_obj_t numerical_minimum(mp_obj_t , mp_obj_t );
-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 *);
-mp_obj_t numerical_sort_inplace(size_t , const mp_obj_t *, mp_map_t *);
-mp_obj_t numerical_argsort(size_t , const mp_obj_t *, mp_map_t *);
-
-// this macro could be tighter, if we moved the ifs to the argmin function, assigned <, as well as >
-#define ARG_MIN_LOOP(in, type, start, stop, stride, op) do {\
-    type *array = (type *)(in)->array->items;\
-    if(((op) == NUMERICAL_MAX) || ((op) == NUMERICAL_ARGMAX)) {\
-        for(size_t i=(start)+(stride); i < (stop); i+=(stride)) {\
-            if((array)[i] > (array)[best_idx]) {\
-                best_idx = i;\
-            }\
-        }\
-    } else{\
-        for(size_t i=(start)+(stride); i < (stop); i+=(stride)) {\
-            if((array)[i] < (array)[best_idx]) best_idx = i;\
-        }\
-    }\
-} while(0)
-
-#define CALCULATE_DIFF(in, out, type, M, N, inn, increment) do {\
-    type *source = (type *)(in)->array->items;\
-    type *target = (type *)(out)->array->items;\
-    for(size_t i=0; i < (M); i++) {\
-        for(size_t j=0; j < (N); j++) {\
-            for(uint8_t k=0; k < n+1; k++) {\
-                target[i*(N)+j] -= stencil[k]*source[i*(inn)+j+k*(increment)];\
-            }\
-        }\
-    }\
-} while(0)
-
-#define HEAPSORT(type, ndarray) do {\
+#define CALCULATE_SUM(ndarray, type, farray, shape_ax, index, stride, offset, optype) do {\
     type *array = (type *)(ndarray)->array->items;\
-    type tmp;\
-    for (;;) {\
-        if (k > 0) {\
-            tmp = array[start+(--k)*increment];\
-        } else {\
-            q--;\
-            if(q == 0) {\
-                break;\
-            }\
-            tmp = array[start+q*increment];\
-            array[start+q*increment] = array[start];\
-        }\
-        p = k;\
-        c = k + k + 1;\
-        while (c < q) {\
-            if((c + 1 < q)  &&  (array[start+(c+1)*increment] > array[start+c*increment])) {\
-                c++;\
-            }\
-            if(array[start+c*increment] > tmp) {\
-                array[start+p*increment] = array[start+c*increment];\
-                p = c;\
-                c = p + p + 1;\
-            } else {\
-                break;\
-            }\
-        }\
-        array[start+p*increment] = tmp;\
+    (farray)[(index)] = 0.0;\
+    for(size_t j=0; j < (shape_ax); j++, (offset) += (stride)) {\
+        (farray)[(index)] += array[(offset)];\
     }\
 } while(0)
 
-// This is pretty similar to HEAPSORT above; perhaps, the two could be combined somehow
-// On the other hand, since this is a macro, it doesn't really matter
-// Keep in mind that initially, index_array[start+s*increment] = s
-#define HEAP_ARGSORT(type, ndarray, index_array) do {\
+// TODO: this can be done without the NDARRAY_INDEX_FROM_FLAT macro
+// Welford algorithm for the standard deviation
+#define CALCULATE_FLAT_SUM_STD(ndarray, type, value, shape_strides, optype) do {\
     type *array = (type *)(ndarray)->array->items;\
-    type tmp;\
-    uint16_t itmp;\
-    for (;;) {\
-        if (k > 0) {\
-            k--;\
-            tmp = array[start+index_array[start+k*increment]*increment];\
-            itmp = index_array[start+k*increment];\
+    (value) = 0.0;\
+    mp_float_t m = 0.0, mtmp;\
+    size_t index, nindex;\
+    for(size_t j=0; j < (ndarray)->len; j++) {\
+        NDARRAY_INDEX_FROM_FLAT((ndarray), (shape_strides), j, index, nindex);\
+        if((optype) == NUMERICAL_STD) {\
+            mtmp = m;\
+            m = mtmp + (array[nindex] - mtmp) / (j+1);\
+            (value) += (array[nindex] - mtmp) * (array[nindex] - m);\
         } else {\
-            q--;\
-            if(q == 0) {\
-                break;\
-            }\
-            tmp = array[start+index_array[start+q*increment]*increment];\
-            itmp = index_array[start+q*increment];\
-            index_array[start+q*increment] = index_array[start];\
-        }\
-        p = k;\
-        c = k + k + 1;\
-        while (c < q) {\
-            if((c + 1 < q)  &&  (array[start+index_array[start+(c+1)*increment]*increment] > array[start+index_array[start+c*increment]*increment])) {\
-                c++;\
-            }\
-            if(array[start+index_array[start+c*increment]*increment] > tmp) {\
-                index_array[start+p*increment] = index_array[start+c*increment];\
-                p = c;\
-                c = p + p + 1;\
-            } else {\
-                break;\
+            (value) += array[nindex];\
         }\
     }\
-        index_array[start+p*increment] = itmp;\
+} while(0)
+
+// we calculate the standard deviation in two passes, in order to avoid negative values through truncation errors
+// We could do in a single pass, if we resorted to the Welford algorithm above
+#define CALCULATE_STD(ndarray, type, sq_sum, shape_ax, stride, offset) do {\
+    type *array = (type *)(ndarray)->array->items;\
+    mp_float_t x, ave = 0.0;\
+    (sq_sum) = 0.0;\
+    size_t j, _offset = (offset);\
+    for(j=0; j < (shape_ax); j++, _offset += (stride)) {\
+        ave += array[_offset];\
+    }\
+    ave /= j;\
+    for(j=0; j < (shape_ax); j++, (offset) += (stride)) {\
+        x = array[(offset)] - ave;\
+        (sq_sum) += x * x;\
     }\
 } while(0)
 

code/poly.c

@@ -15,102 +15,62 @@
 #include "linalg.h"
 #include "poly.h"
 
-
-bool object_is_nditerable(mp_obj_t o_in) {
-    if(mp_obj_is_type(o_in, &ulab_ndarray_type) || 
-      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)) {
-        return true;
-    }
-    return false;
-}
-
-size_t get_nditerable_len(mp_obj_t o_in) {
-    if(mp_obj_is_type(o_in, &ulab_ndarray_type)) {
-        ndarray_obj_t *in = MP_OBJ_TO_PTR(o_in);
-        return in->array->len;
-    } else {
-        return (size_t)mp_obj_get_int(mp_obj_len_maybe(o_in));
+void fill_array_iterable(mp_float_t *array, mp_obj_t oin) {
+    mp_obj_iter_buf_t buf;
+    mp_obj_t item, iterable = mp_getiter(oin, &buf);
+    while((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
+        *array++ = mp_obj_get_float(item);
     }
 }
 
 mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {
-    // TODO: return immediately, if o_p is not an iterable
-    // TODO: there is a bug here: matrices won't work, 
-    // because there is a single iteration loop
-    size_t m, n;
-    if(MP_OBJ_IS_TYPE(o_x, &ulab_ndarray_type)) {
-        ndarray_obj_t *ndx = MP_OBJ_TO_PTR(o_x);
-        m = ndx->m;
-        n = ndx->n;
-    } else {
-        mp_obj_array_t *ix = MP_OBJ_TO_PTR(o_x);
-        m = 1;
-        n = ix->len;
-    }
-    // polynomials are going to be of type float, except, when both 
-    // the coefficients and the independent variable are integers
-    ndarray_obj_t *out = create_new_ndarray(m, n, NDARRAY_FLOAT);
-    mp_obj_iter_buf_t x_buf;
-    mp_obj_t x_item, x_iterable = mp_getiter(o_x, &x_buf);
-
-    mp_obj_iter_buf_t p_buf;
-    mp_obj_t p_item, p_iterable;
-
-    mp_float_t x, y;
-    mp_float_t *outf = (mp_float_t *)out->array->items;
+    // we always return floats: polynomials are going to be of type float, except, 
+    // when both the coefficients and the independent variable are integers; 
     uint8_t plen = mp_obj_get_int(mp_obj_len_maybe(o_p));
     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] = mp_obj_get_float(p_item);
-        i++;
+    fill_array_iterable(p, o_p);
+    ndarray_obj_t *ndarray;
+    mp_float_t *array;
+    if(MP_OBJ_IS_TYPE(o_x, &ulab_ndarray_type)) {
+        ndarray_obj_t *input = MP_OBJ_TO_PTR(o_x);
+        ndarray = ndarray_copy_view(input, NDARRAY_FLOAT);
+        array = (mp_float_t *)ndarray->array->items;
+    } else { // at this point, we should have a 1-D iterable
+        size_t len = mp_obj_get_int(mp_obj_len_maybe(o_x));
+        ndarray = ndarray_new_linear_array(len, NDARRAY_FLOAT);
+        array = (mp_float_t *)ndarray->array->items;
+        fill_array_iterable(array, o_x);
     }
-    i = 0;
-    while ((x_item = mp_iternext(x_iterable)) != MP_OBJ_STOP_ITERATION) {
-        x = mp_obj_get_float(x_item);
+    mp_float_t x, y;
+    for(size_t i=0; i < ndarray->len; i++) {
+        x = array[i];
         y = p[0];
         for(uint8_t j=0; j < plen-1; j++) {
             y *= x;
             y += p[j+1];
         }
-        outf[i++] = y;
+        array[i] = y;
     }
     m_del(mp_float_t, p, plen);
-    return MP_OBJ_FROM_PTR(out);
+    return MP_OBJ_FROM_PTR(ndarray);
 }
 
 mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {
-    if((n_args != 2) && (n_args != 3)) {
-        mp_raise_ValueError("number of arguments must be 2, or 3");
+    if(n_args != 3) {
+        mp_raise_ValueError("number of arguments must be 3");
     }
-    if(!object_is_nditerable(args[0])) {
-        mp_raise_ValueError("input data must be an iterable");
+    if(!MP_OBJ_IS_TYPE(args[0], &ulab_ndarray_type) && !MP_OBJ_IS_TYPE(args[0], &mp_type_tuple) &&
+      !MP_OBJ_IS_TYPE(args[0], &mp_type_list) && !MP_OBJ_IS_TYPE(args[0], &mp_type_range) &&
+      !MP_OBJ_IS_TYPE(args[1], &ulab_ndarray_type) && !MP_OBJ_IS_TYPE(args[1], &mp_type_tuple) &&
+      !MP_OBJ_IS_TYPE(args[1], &mp_type_list) && !MP_OBJ_IS_TYPE(args[1], &mp_type_range)) {
+        mp_raise_ValueError("input data must be a 1D iterable");
     }
-    uint16_t lenx = 0, leny = 0;
-    uint8_t deg = 0;
+    uint16_t lenx, leny;
+    uint8_t deg;
     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, 
-        // in which case we are between the rock and a hard place
-        leny = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));
-        deg = (uint8_t)mp_obj_get_int(args[1]);
-        if(leny < deg) {
-            mp_raise_ValueError("more degrees of freedom than data points");
-        }
-        lenx = leny;
-        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(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]));
-        leny = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));
+    leny = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[1]));
     if(lenx != leny) {
         mp_raise_ValueError("input vectors must be of equal length");
     }
@@ -122,7 +82,6 @@ mp_obj_t poly_polyfit(size_t  n_args, const mp_obj_t *args) {
     fill_array_iterable(x, args[0]);
     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
@@ -170,7 +129,7 @@ mp_obj_t poly_polyfit(size_t  n_args, const mp_obj_t *args) {
     // XT is no longer needed
     m_del(mp_float_t, XT, (deg+1)*leny);
     
-    ndarray_obj_t *beta = create_new_ndarray(deg+1, 1, NDARRAY_FLOAT);
+    ndarray_obj_t *beta = ndarray_new_linear_array(deg+1, NDARRAY_FLOAT);
     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);

code/ulab.c

@@ -18,13 +18,13 @@
 #include "py/objarray.h"
 
 #include "ndarray.h"
-#include "linalg.h"
 #include "vectorise.h"
+#include "linalg.h"
 #include "poly.h"
-#include "fft.h"
 #include "numerical.h"
+#include "fft.h"
 
-#define ULAB_VERSION 0.262
+#define ULAB_VERSION 1.0
 
 typedef struct _mp_obj_float_t {
     mp_obj_base_t base;
@@ -34,20 +34,12 @@ typedef struct _mp_obj_float_t {
 mp_obj_float_t ulab_version = {{&mp_type_float}, ULAB_VERSION};
 
 MP_DEFINE_CONST_FUN_OBJ_1(ndarray_shape_obj, ndarray_shape);
-MP_DEFINE_CONST_FUN_OBJ_1(ndarray_rawsize_obj, ndarray_rawsize);
+MP_DEFINE_CONST_FUN_OBJ_1(ndarray_strides_obj, ndarray_strides);
+MP_DEFINE_CONST_FUN_OBJ_2(ndarray_reshape_obj, ndarray_reshape);
+MP_DEFINE_CONST_FUN_OBJ_1(ndarray_transpose_obj, ndarray_transpose);
 MP_DEFINE_CONST_FUN_OBJ_KW(ndarray_flatten_obj, 1, ndarray_flatten);
-MP_DEFINE_CONST_FUN_OBJ_1(ndarray_asbytearray_obj, ndarray_asbytearray);
-
-MP_DEFINE_CONST_FUN_OBJ_1(linalg_transpose_obj, linalg_transpose);
-MP_DEFINE_CONST_FUN_OBJ_2(linalg_reshape_obj, linalg_reshape);
-MP_DEFINE_CONST_FUN_OBJ_KW(linalg_size_obj, 1, linalg_size);
-MP_DEFINE_CONST_FUN_OBJ_1(linalg_inv_obj, linalg_inv);
-MP_DEFINE_CONST_FUN_OBJ_2(linalg_dot_obj, linalg_dot);
-MP_DEFINE_CONST_FUN_OBJ_KW(linalg_zeros_obj, 0, linalg_zeros);
-MP_DEFINE_CONST_FUN_OBJ_KW(linalg_ones_obj, 0, linalg_ones);
-MP_DEFINE_CONST_FUN_OBJ_KW(linalg_eye_obj, 0, linalg_eye);
-MP_DEFINE_CONST_FUN_OBJ_1(linalg_det_obj, linalg_det);
-MP_DEFINE_CONST_FUN_OBJ_1(linalg_eig_obj, linalg_eig);
+MP_DEFINE_CONST_FUN_OBJ_1(ndarray_itemsize_obj, ndarray_itemsize);
+MP_DEFINE_CONST_FUN_OBJ_1(ndarray_info_obj, ndarray_info);
 
 MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acos_obj, vectorise_acos);
 MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acosh_obj, vectorise_acosh);
@@ -73,36 +65,34 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sqrt_obj, vectorise_sqrt);
 MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tan_obj, vectorise_tan);
 MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tanh_obj, vectorise_tanh);
 
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_linspace_obj, 2, numerical_linspace);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sum_obj, 1, numerical_sum);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_mean_obj, 1, numerical_mean);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_std_obj, 1, numerical_std);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_min_obj, 1, numerical_min);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_max_obj, 1, numerical_max);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argmin_obj, 1, numerical_argmin);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argmax_obj, 1, numerical_argmax);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_roll_obj, 2, numerical_roll);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_flip_obj, 1, numerical_flip);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_diff_obj, 1, numerical_diff);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sort_obj, 1, numerical_sort);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sort_inplace_obj, 1, numerical_sort_inplace);
-STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argsort_obj, 1, numerical_argsort);
+MP_DEFINE_CONST_FUN_OBJ_1(linalg_inv_obj, linalg_inv);
+MP_DEFINE_CONST_FUN_OBJ_2(linalg_dot_obj, linalg_dot);
+MP_DEFINE_CONST_FUN_OBJ_KW(linalg_zeros_obj, 0, linalg_zeros);
+MP_DEFINE_CONST_FUN_OBJ_KW(linalg_ones_obj, 0, linalg_ones);
+MP_DEFINE_CONST_FUN_OBJ_KW(linalg_eye_obj, 0, linalg_eye);
+MP_DEFINE_CONST_FUN_OBJ_1(linalg_det_obj, linalg_det);
+MP_DEFINE_CONST_FUN_OBJ_1(linalg_eig_obj, linalg_eig);
 
-STATIC MP_DEFINE_CONST_FUN_OBJ_2(poly_polyval_obj, poly_polyval);
-STATIC MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(poly_polyfit_obj, 2, 3, poly_polyfit);
+MP_DEFINE_CONST_FUN_OBJ_2(poly_polyval_obj, poly_polyval);
+MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(poly_polyfit_obj, 2, 3, poly_polyfit);
 
-STATIC MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj, 1, 2, fft_fft);
-STATIC MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_ifft_obj, 1, 2, fft_ifft);
-STATIC MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_spectrum_obj, 1, 2, fft_spectrum);
+MP_DEFINE_CONST_FUN_OBJ_KW(numerical_linspace_obj, 2, numerical_linspace);
+MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sum_obj, 1, numerical_sum);
+MP_DEFINE_CONST_FUN_OBJ_KW(numerical_mean_obj, 1, numerical_mean);
+MP_DEFINE_CONST_FUN_OBJ_KW(numerical_std_obj, 1, numerical_std);
+
+
+MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj, 1, 2, fft_fft);
+MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_ifft_obj, 1, 2, fft_ifft);
+MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_spectrum_obj, 1, 2, fft_spectrum);
 
 STATIC const mp_rom_map_elem_t ulab_ndarray_locals_dict_table[] = {
     { MP_ROM_QSTR(MP_QSTR_shape), MP_ROM_PTR(&ndarray_shape_obj) },
-    { MP_ROM_QSTR(MP_QSTR_rawsize), MP_ROM_PTR(&ndarray_rawsize_obj) },
+    { MP_ROM_QSTR(MP_QSTR_strides), MP_ROM_PTR(&ndarray_strides_obj) },
+    { MP_ROM_QSTR(MP_QSTR_reshape), MP_ROM_PTR(&ndarray_reshape_obj) },
+    { MP_ROM_QSTR(MP_QSTR_transpose), MP_ROM_PTR(&ndarray_transpose_obj) },
     { MP_ROM_QSTR(MP_QSTR_flatten), MP_ROM_PTR(&ndarray_flatten_obj) },
-    { MP_ROM_QSTR(MP_QSTR_asbytearray), MP_ROM_PTR(&ndarray_asbytearray_obj) },
-    { MP_ROM_QSTR(MP_QSTR_transpose), MP_ROM_PTR(&linalg_transpose_obj) },
-    { MP_ROM_QSTR(MP_QSTR_reshape), MP_ROM_PTR(&linalg_reshape_obj) },
-    { MP_ROM_QSTR(MP_QSTR_sort), MP_ROM_PTR(&numerical_sort_inplace_obj) },
+    { MP_ROM_QSTR(MP_QSTR_itemsize), MP_ROM_PTR(&ndarray_itemsize_obj) },
 };
 
 STATIC MP_DEFINE_CONST_DICT(ulab_ndarray_locals_dict, ulab_ndarray_locals_dict_table);
@@ -123,14 +113,7 @@ STATIC const mp_map_elem_t ulab_globals_table[] = {
     { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_ulab) },
     { MP_ROM_QSTR(MP_QSTR___version__), MP_ROM_PTR(&ulab_version) },
     { MP_OBJ_NEW_QSTR(MP_QSTR_array), (mp_obj_t)&ulab_ndarray_type },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_size), (mp_obj_t)&linalg_size_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_inv), (mp_obj_t)&linalg_inv_obj },
-    { MP_ROM_QSTR(MP_QSTR_dot), (mp_obj_t)&linalg_dot_obj },
-    { MP_ROM_QSTR(MP_QSTR_zeros), (mp_obj_t)&linalg_zeros_obj },
-    { MP_ROM_QSTR(MP_QSTR_ones), (mp_obj_t)&linalg_ones_obj },
-    { MP_ROM_QSTR(MP_QSTR_eye), (mp_obj_t)&linalg_eye_obj },
-    { MP_ROM_QSTR(MP_QSTR_det), (mp_obj_t)&linalg_det_obj },
-    { MP_ROM_QSTR(MP_QSTR_eig), (mp_obj_t)&linalg_eig_obj },    
+    { MP_OBJ_NEW_QSTR(MP_QSTR_ndinfo), (mp_obj_t)&ndarray_info_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_acos), (mp_obj_t)&vectorise_acos_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_acosh), (mp_obj_t)&vectorise_acosh_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_asin), (mp_obj_t)&vectorise_asin_obj },
@@ -154,25 +137,24 @@ STATIC const mp_map_elem_t ulab_globals_table[] = {
     { MP_OBJ_NEW_QSTR(MP_QSTR_sqrt), (mp_obj_t)&vectorise_sqrt_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_tan), (mp_obj_t)&vectorise_tan_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_tanh), (mp_obj_t)&vectorise_tanh_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_polyval), (mp_obj_t)&poly_polyval_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_polyfit), (mp_obj_t)&poly_polyfit_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_inv), (mp_obj_t)&linalg_inv_obj },
+    { MP_ROM_QSTR(MP_QSTR_dot), (mp_obj_t)&linalg_dot_obj },
+    { MP_ROM_QSTR(MP_QSTR_zeros), (mp_obj_t)&linalg_zeros_obj },
+    { MP_ROM_QSTR(MP_QSTR_ones), (mp_obj_t)&linalg_ones_obj },
+    { MP_ROM_QSTR(MP_QSTR_eye), (mp_obj_t)&linalg_eye_obj },
+    { MP_ROM_QSTR(MP_QSTR_det), (mp_obj_t)&linalg_det_obj },
+    { MP_ROM_QSTR(MP_QSTR_eig), (mp_obj_t)&linalg_eig_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_linspace), (mp_obj_t)&numerical_linspace_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_sum), (mp_obj_t)&numerical_sum_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_mean), (mp_obj_t)&numerical_mean_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_std), (mp_obj_t)&numerical_std_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_min), (mp_obj_t)&numerical_min_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_max), (mp_obj_t)&numerical_max_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_argmin), (mp_obj_t)&numerical_argmin_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_argmax), (mp_obj_t)&numerical_argmax_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_roll), (mp_obj_t)&numerical_roll_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_flip), (mp_obj_t)&numerical_flip_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_diff), (mp_obj_t)&numerical_diff_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_sort), (mp_obj_t)&numerical_sort_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_argsort), (mp_obj_t)&numerical_argsort_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_polyval), (mp_obj_t)&poly_polyval_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_polyfit), (mp_obj_t)&poly_polyfit_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_fft), (mp_obj_t)&fft_fft_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_ifft), (mp_obj_t)&fft_ifft_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_spectrum), (mp_obj_t)&fft_spectrum_obj },
     // class constants
+    { MP_ROM_QSTR(MP_QSTR_bool), MP_ROM_INT(NDARRAY_BOOL) },
     { MP_ROM_QSTR(MP_QSTR_uint8), MP_ROM_INT(NDARRAY_UINT8) },
     { MP_ROM_QSTR(MP_QSTR_int8), MP_ROM_INT(NDARRAY_INT8) },
     { MP_ROM_QSTR(MP_QSTR_uint16), MP_ROM_INT(NDARRAY_UINT16) },

code/vectorise.c

@@ -26,34 +26,41 @@ mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) {
     if(mp_obj_is_float(o_in) || mp_obj_is_integer(o_in)) {
         return mp_obj_new_float(f(mp_obj_get_float(o_in)));
     }
-    mp_float_t x;
     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);
-        mp_float_t *dataout = (mp_float_t *)ndarray->array->items;
+        ndarray_obj_t *ndarray;
+        if(source->len == source->array->len) { // this is a dense array
+            ndarray = ndarray_new_dense_ndarray(source->ndim, source->shape, NDARRAY_FLOAT);
+            mp_float_t *array = (mp_float_t *)ndarray->array->items;
             if(source->array->typecode == NDARRAY_UINT8) {
-            ITERATE_VECTOR(uint8_t, source, dataout);
+                ITERATE_VECTOR(uint8_t, source, array);
             } else if(source->array->typecode == NDARRAY_INT8) {
-            ITERATE_VECTOR(int8_t, source, dataout);
+                ITERATE_VECTOR(int8_t, source, array);
             } else if(source->array->typecode == NDARRAY_UINT16) {
-            ITERATE_VECTOR(uint16_t, source, dataout);
+                ITERATE_VECTOR(uint16_t, source, array);
             } else if(source->array->typecode == NDARRAY_INT16) {
-            ITERATE_VECTOR(int16_t, source, dataout);
+                ITERATE_VECTOR(int16_t, source, array);
             } else {
-            ITERATE_VECTOR(mp_float_t, source, dataout);
+                ITERATE_VECTOR(mp_float_t, source, array);
+            } 
+        } else {
+            ndarray = ndarray_copy_view(source, NDARRAY_FLOAT);
+            mp_float_t *dataout = (mp_float_t *)ndarray->array->items;
+            ITERATE_VECTOR(mp_float_t, ndarray, 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_IS_TYPE(o_in, &mp_type_range)) { // i.e., the input is a generic, one-dimensional iterable
             mp_obj_array_t *o = MP_OBJ_TO_PTR(o_in);
-            ndarray_obj_t *out = create_new_ndarray(1, o->len, NDARRAY_FLOAT);
+            ndarray_obj_t *out = ndarray_new_linear_array(o->len, NDARRAY_FLOAT);
             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;
+            mp_float_t x;
             while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
                 x = mp_obj_get_float(item);
-                dataout[i++] = f(x);
+                *dataout++ = f(x);
             }
         return MP_OBJ_FROM_PTR(out);
     }

code/vectorise.h

@@ -40,7 +40,7 @@ mp_obj_t vectorise_tanh(mp_obj_t );
 #define ITERATE_VECTOR(type, source, out) do {\
     type *input = (type *)(source)->array->items;\
     for(size_t i=0; i < (source)->array->len; i++) {\
-                (out)[i] = f(input[i]);\
+        *(out)++ = f(*input++);\
     }\
 } while(0)