added argsort, and fixed a number of smaller bugs

README.md

@@ -8,9 +8,16 @@ Documentation can be found under https://micropython-ulab.readthedocs.io/en/late
 The source for the manual is in https://github.com/v923z/micropython-ulab/blob/master/docs/ulab-manual.ipynb,
 while developer help is in https://github.com/v923z/micropython-ulab/blob/master/docs/ulab.ipynb.
 
+# Firmware
+
 Firmware for pyboard.v.1.1, and PYBD_SF6 is updated once in a while, and can be downloaded 
-from https://github.com/v923z/micropython-ulab/releases, otherwise, it can be 
-compiled from the source by following these steps:
+from https://github.com/v923z/micropython-ulab/releases.
+
+## Compiling
+
+If you want to try the latest version of `ulab`, or your hardware is 
+different to pyboard.v.1.1, or PYBD_SF6, the firmware can be compiled 
+from the source by following these steps:
 
 First, you have to clone the micropython repository by running 
 
@@ -38,7 +45,12 @@ If that was successful, you can try to run the make command in the port's direct
 ```
 make BOARD=PYBV11 USER_C_MODULES=../../../ulab all
 ```
-which will prepare the firmware for pyboard.v.11. Provided that you managed to compile the firmware, you would upload that either by running 
+which will prepare the firmware for pyboard.v.11. Similarly, 
+```
+make BOARD=PYBD_SF6 USER_C_MODULES=../../../ulab all
+```
+will compile for the SF6 member of the PYBD series. Provided that you managed to compile the firmware, you would upload that by running
+either
 ```
 dfu-util --alt 0 -D firmware.dfu
 ```

code/fft.c

@@ -19,7 +19,6 @@
 #include "ndarray.h"
 #include "fft.h"
 
-
 enum FFT_TYPE {
     FFT_FFT,
     FFT_IFFT,
@@ -37,8 +36,8 @@ void fft_kernel(mp_float_t *real, mp_float_t *imag, int n, int isign) {
     j = 0;
     for(int i = 0; i < n; i++) {
         if (j > i) {
-            SWAP(float, real[i], real[j]);
-            SWAP(float, imag[i], imag[j]);            
+            SWAP(mp_float_t, real[i], real[j]);
+            SWAP(mp_float_t, imag[i], imag[j]);
         }
         m = n >> 1;
         while (j >= m && m > 0) {
@@ -51,7 +50,7 @@ void fft_kernel(mp_float_t *real, mp_float_t *imag, int n, int isign) {
     mmax = 1;
     while (n > mmax) {
         istep = mmax << 1;
-        theta = -1.0*isign*6.28318530717959/istep;
+        theta = -2.0*isign*MP_PI/istep;
         wtemp = MICROPY_FLOAT_C_FUN(sin)(0.5 * theta);
         wpr = -2.0 * wtemp * wtemp;
         wpi = MICROPY_FLOAT_C_FUN(sin)(theta);

code/numerical.c

@@ -539,38 +539,34 @@ mp_obj_t numerical_diff(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
     return MP_OBJ_FROM_PTR(out);
 }
 
-mp_obj_t numerical_sort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
-
-    static const mp_arg_t allowed_args[] = {
-        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } },
-        { MP_QSTR_n, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 1 } },
-        { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_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_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");
     }
-
-    mp_obj_t out = ndarray_copy(args[0].u_obj);
-    ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(out);
+
+    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(args[2].u_obj == mp_const_none) { // flatten the array
+    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(args[2].u_obj) == -1) || 
-              (mp_obj_get_int(args[2].u_obj) == 1)) { // sort along the horizontal axis
+    } 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(args[2].u_obj) == 0) { // sort along vertical axis
+    } else if(mp_obj_get_int(axis) == 0) { // sort along vertical axis
         increment = ndarray->n;
         start_inc = 1;
         end = ndarray->m;
@@ -592,5 +588,104 @@ mp_obj_t numerical_sort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
             HEAPSORT(mp_float_t, ndarray);
         }
     }
-    return out;
+    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

@@ -30,6 +30,8 @@ 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 {\
@@ -47,7 +49,6 @@ mp_obj_t numerical_sort(size_t , const mp_obj_t *, mp_map_t *);
     }\
 } 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;\
@@ -61,9 +62,9 @@ mp_obj_t numerical_sort(size_t , const mp_obj_t *, mp_map_t *);
 } while(0)
 
 #define HEAPSORT(type, ndarray) do {\
-    uint8_t *array = (uint8_t *)(ndarray)->array->items;\
-    uint8_t tmp;\
-    for (; q > 0;) {\
+    type *array = (type *)(ndarray)->array->items;\
+    type tmp;\
+    for (;;) {\
         if (k > 0) {\
             tmp = array[start+(--k)*increment];\
         } else {\
@@ -92,4 +93,43 @@ mp_obj_t numerical_sort(size_t , const mp_obj_t *, mp_map_t *);
     }\
 } 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 {\
+    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];\
+        } 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;\
+            }\
+        }\
+        index_array[start+p*increment] = itmp;\
+    }\
+} while(0)
+
 #endif

code/ulab.c

@@ -85,6 +85,8 @@ 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);
 
 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);
@@ -100,6 +102,7 @@ STATIC const mp_rom_map_elem_t ulab_ndarray_locals_dict_table[] = {
     { 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) },
 };
 
 STATIC MP_DEFINE_CONST_DICT(ulab_ndarray_locals_dict, ulab_ndarray_locals_dict_table);
@@ -163,6 +166,7 @@ STATIC const mp_map_elem_t ulab_globals_table[] = {
     { 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 },

code/vectorise.c

@@ -60,27 +60,6 @@ mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) {
     return mp_const_none;
 }
 
-// _degrees won't compile for the unix port
-/*
-mp_float_t _degreesf(mp_float_t x) {
-    return(180*x/MP_PI);
-}
-
-MATH_FUN_1(degrees, _degrees);
-
-// _radians won't compile for the unix port
-mp_float_t _radiansf(mp_float_t x) {
-    return(MP_PI*x/180.0);
-}
-
-MATH_FUN_1(radians, _radians);
-
-STATIC mp_float_t _fabsf(mp_float_t x) {
-    return fabsf(x);
-}
-
-MATH_FUN_1(fabs, _fabs);
-*/
 MATH_FUN_1(acos, acos);
 MATH_FUN_1(acosh, acosh);
 MATH_FUN_1(asin, asin);

docs/manual/source/conf.py

@@ -22,7 +22,7 @@
 author = 'Zoltán Vörös'
 
 # The full version, including alpha/beta/rc tags
-release = '0.252'
+release = '0.26'
 
 
 # -- General configuration ---------------------------------------------------