jepler/circuitpython-ulab
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
circuitpython-ulab/docs/ulab.ipynb
Zoltán Vörös e0e840f6d5 added circuitpython-related stuff to the manual
2020-02-26 18:05:49 +01:00

5114 lines
203 KiB
Text
Raw Permalink Blame History

{
"cells": [
{
"cell_type": "code",
"execution_count": 612,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-16T17:00:37.079043Z",
"start_time": "2020-02-16T17:00:34.761337Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"%pylab inline"
]
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"source": [
"# Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {
"hidden": true
},
"source": [
"ulab is a C module for micropython. My goal was to implement a small subset of numpy. I chose those functions and methods that might be useful in the context of a microcontroller. This means low-level data processing of linear (array) and two-dimensional (matrix) data.\n",
"\n",
"The main points of ulab are \n",
"\n",
"- compact, iterable and slicable container of numerical data in 1, and 2 dimensions (arrays and matrices). In addition, these containers support all the relevant unary and binary operators (e.g., `len`, ==, +, *, etc.)\n",
"- vectorised computations on micropython iterables and numerical arrays/matrices (universal functions)\n",
"- basic linear algebra routines (matrix inversion, matrix reshaping, and transposition)\n",
"- polynomial fits to numerical data\n",
"- fast Fourier transforms\n",
"\n",
"The code itself is split into submodules. This should make exclusion of unnecessary functions, if storage space is a concern. Each section of the implementation part kicks out with a short discussion on what can be done with the particular submodule, and what are the tripping points at the C level. I hope that these musings can be used as a starting point for further discussion on the code.\n",
"\n",
"The code and its documentation can be found under https://github.com/v923z/micropython-ulab/. The MIT licence applies to all material."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Environmental settings and magic commands\n",
"\n",
"The entire C source code, as well as the documentation (mainly verbose comments on certain aspects of the implementation) are contained in this notebook. The code is exported to separate C files in `/ulab/`, and then compiled from this notebook. However, I would like to stress that the compilation does not require a jupyter notebook. It can be done from the command line by invoking the command in the [make](#make), or [Compiling the module](#Compiling-the-module). After all, the ipython kernel simply passes the `make` commands to the underlying operating system.\n",
"\n",
"Testing is done on the unix and stm32 ports of micropython, also directly from the notebook. This is why this section contains a couple of magic functions. But once again: the C module can be used without the notebook. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T05:53:50.706261Z",
"start_time": "2020-02-26T05:53:50.694240Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/v923z/sandbox/micropython/v1.11/micropython/ports/unix\n"
]
}
],
"source": [
"%cd ../../micropython/ports/unix/"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T05:50:26.530343Z",
"start_time": "2020-02-26T05:50:26.521798Z"
}
},
"outputs": [],
"source": [
"from IPython.core.magic import Magics, magics_class, line_cell_magic\n",
"from IPython.core.magic import cell_magic, register_cell_magic, register_line_magic\n",
"from IPython.core.magic_arguments import argument, magic_arguments, parse_argstring\n",
"import subprocess\n",
"import os"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## micropython magic command\n",
"\n",
"The following magic class takes the content of a cell, and depending on the arguments, either passes it to the unix, or the stm32 implementation. In the latter case, a pyboard must be connected to the computer, and must be initialised beforehand. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T05:50:28.995120Z",
"start_time": "2020-02-26T05:50:28.979972Z"
}
},
"outputs": [],
"source": [
"@magics_class\n",
"class PyboardMagic(Magics):\n",
" @cell_magic\n",
" @magic_arguments()\n",
" @argument('-skip')\n",
" @argument('-unix')\n",
" @argument('-file')\n",
" @argument('-data')\n",
" @argument('-time')\n",
" @argument('-memory')\n",
" def micropython(self, line='', cell=None):\n",
" args = parse_argstring(self.micropython, line)\n",
" if args.skip: # doesn't care about the cell's content\n",
" print('skipped execution')\n",
" return None # do not parse the rest\n",
" if args.unix: # tests the code on the unix port. Note that this works on unix only\n",
" with open('/dev/shm/micropython.py', 'w') as fout:\n",
" fout.write(cell)\n",
" proc = subprocess.Popen([\"./micropython\", \"/dev/shm/micropython.py\"], \n",
" stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
" print(proc.stdout.read().decode(\"utf-8\"))\n",
" print(proc.stderr.read().decode(\"utf-8\"))\n",
" return None\n",
" if args.file: # can be used to copy the cell content onto the pyboard's flash\n",
" spaces = \" \"\n",
" try:\n",
" with open(args.file, 'w') as fout:\n",
" fout.write(cell.replace('\\t', spaces))\n",
" printf('written cell to {}'.format(args.file))\n",
" except:\n",
" print('Failed to write to disc!')\n",
" return None # do not parse the rest\n",
" if args.data: # can be used to load data from the pyboard directly into kernel space\n",
" message = pyb.exec(cell)\n",
" if len(message) == 0:\n",
" print('pyboard >>>')\n",
" else:\n",
" print(message.decode('utf-8'))\n",
" # register new variable in user namespace\n",
" self.shell.user_ns[args.data] = string_to_matrix(message.decode(\"utf-8\"))\n",
" \n",
" if args.time: # measures the time of executions\n",
" pyb.exec('import utime')\n",
" message = pyb.exec('t = utime.ticks_us()\\n' + cell + '\\ndelta = utime.ticks_diff(utime.ticks_us(), t)' + \n",
" \"\\nprint('execution time: {:d} us'.format(delta))\")\n",
" print(message.decode('utf-8'))\n",
" \n",
" if args.memory: # prints out memory information \n",
" message = pyb.exec('from micropython import mem_info\\nprint(mem_info())\\n')\n",
" print(\"memory before execution:\\n========================\\n\", message.decode('utf-8'))\n",
" message = pyb.exec(cell)\n",
" print(\">>> \", message.decode('utf-8'))\n",
" message = pyb.exec('print(mem_info())')\n",
" print(\"memory after execution:\\n========================\\n\", message.decode('utf-8'))\n",
"\n",
" else:\n",
" message = pyb.exec(cell)\n",
" print(message.decode('utf-8'))\n",
"\n",
"ip = get_ipython()\n",
"ip.register_magics(PyboardMagic)"
]
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"source": [
"### pyboard initialisation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-09-24T18:07:26.532228Z",
"start_time": "2019-09-24T18:07:26.398434Z"
},
"hidden": true
},
"outputs": [],
"source": [
"import pyboard\n",
"pyb = pyboard.Pyboard('/dev/ttyACM0')\n",
"pyb.enter_raw_repl()"
]
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"source": [
"### pyboad detach"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"pyb.exit_raw_repl()\n",
"pyb.close()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"ExecuteTime": {
"end_time": "2019-09-27T05:56:41.168185Z",
"start_time": "2019-09-27T05:56:41.162486Z"
},
"hidden": true
},
"outputs": [
{
"data": {
"application/javascript": [
"\n",
" IPython.CodeCell.options_default.highlight_modes['magic_text/x-csrc'] = {'reg':[/^\\s*%%ccode/]};\n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"\n",
" (function () {\n",
" var defaults = IPython.CodeCell.config_defaults || IPython.CodeCell.options_default;\n",
" defaults.highlight_modes['magic_text/x-csrc'] = {'reg':[/^\\s*%%makefile/]};\n",
" })();\n",
" "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import IPython\n",
"\n",
"js = \"\"\"\n",
" (function () {\n",
" var defaults = IPython.CodeCell.config_defaults || IPython.CodeCell.options_default;\n",
" defaults.highlight_modes['magic_text/x-csrc'] = {'reg':[/^\\\\s*%%ccode/]};\n",
" })();\n",
" \"\"\"\n",
"cjs = \"\"\"\n",
" IPython.CodeCell.options_default.highlight_modes['magic_text/x-csrc'] = {'reg':[/^\\\\s*%%ccode/]};\n",
" \"\"\"\n",
"\n",
"IPython.core.display.display_javascript(cjs, raw=True)\n",
"\n",
"js = \"\"\"\n",
" (function () {\n",
" var defaults = IPython.CodeCell.config_defaults || IPython.CodeCell.options_default;\n",
" defaults.highlight_modes['magic_text/x-csrc'] = {'reg':[/^\\\\s*%%makefile/]};\n",
" })();\n",
" \"\"\"\n",
"IPython.core.display.display_javascript(js, raw=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Code magic\n",
"\n",
"The following cell magic simply writes a licence header, and the contents of the cell to the file given in the header of the cell. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T05:50:31.256457Z",
"start_time": "2020-02-26T05:50:31.239429Z"
}
},
"outputs": [],
"source": [
"@magics_class\n",
"class MyMagics(Magics):\n",
" \n",
" @cell_magic\n",
" def ccode(self, line, cell):\n",
" if line:\n",
" with open('../../../ulab/code/'+line, 'w') as cout:\n",
" cout.write(cell)\n",
" print('written %d bytes to %s'%(len(cell), line))\n",
" return None\n",
" \n",
" @cell_magic\n",
" def rcode(self, line, cell=None):\n",
" with open('../../../ulab/code/'+line, 'r') as fin:\n",
" content = fin.read()\n",
" self.shell.set_next_input('%%ccode {}\\n{}'.format(line, content), replace=True)\n",
" return None\n",
" \n",
"ip = get_ipython()\n",
"ip.register_magics(MyMagics)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Functions included in the `ulab` firmware\n",
"\n",
"This header file defines the functions that will be included in the compiled `ulab` firmware. This header should be included in all submodules. Variable names follow a common pattern: the first part is `ULAB`, the second is the submodule, in which a particular function can be found, and the third is the function's name. \n",
"\n",
"Set the value of the variable to 1, if you want to include the function, otherwise, to 0."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ulab.h"
]
},
{
"cell_type": "code",
"execution_count": 501,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-16T14:41:11.454860Z",
"start_time": "2020-02-16T14:41:11.442473Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 709 bytes to ulab.h\n"
]
}
],
"source": [
"%%ccode ulab.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#ifndef __ULAB__\n",
"#define __ULAB__\n",
"\n",
"// vectorise (all functions) takes approx. 3 kB of flash space\n",
"#define ULAB_VECTORISE_MODULE (1)\n",
"\n",
"// linalg adds around 6 kB\n",
"#define ULAB_LINALG_MODULE (1)\n",
"\n",
"// poly is approx. 2.5 kB\n",
"#define ULAB_POLY_MODULE (1)\n",
"\n",
"// numerical is about 12 kB\n",
"#define ULAB_NUMERICAL_MODULE (1)\n",
"\n",
"// FFT costs about 2 kB of flash space\n",
"#define ULAB_FFT_MODULE (1)\n",
"\n",
"// the filter module takes about 1 kB of flash space\n",
"#define ULAB_FILTER_MODULE (1)\n",
"\n",
"// user-defined modules\n",
"#define ULAB_EXTRAS_MODULE (0)\n",
"\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The ndarray type"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## General comments\n",
"\n",
"`ndarrays` are efficient containers of numerical data of the same type (i.e., signed/unsigned chars, signed/unsigned integers or floats). Beyond storing the actual data, the type definition has three additional members (on top of the `base` type). Namely, two `size_t` objects, `m`, and `n`, which give the dimensions of the matrix (obviously, if the `ndarray` is meant to be linear, either `m`, or `n` is equal to 1), as well as the byte size, `bytes`, i.e., the total number of bytes consumed by the data container. `bytes` is equal to `m*n` for `byte` types (`uint8`, and `int8`), to `2*m*n` for integers (`uint16`, and `int16`), and `4*m*n` for floats. \n",
"\n",
"The type definition is as follows:\n",
"\n",
"```c\n",
"typedef struct _ndarray_obj_t {\n",
" mp_obj_base_t base;\n",
" size_t m, n;\n",
" mp_obj_array_t *array;\n",
" size_t bytes;\n",
"} ndarray_obj_t;\n",
"```\n",
"\n",
"**NOTE: with a little bit of extra effort, mp_obj_array_t can be replaced by a single void array. We should, perhaps, consider the pros and cons of that. One patent advantage is that we could get rid of the verbatim copy of array_new function in ndarray.c. On the other hand, objarray.c facilities couldn't be used anymore.**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Additional structure members in numpy\n",
"\n",
"Also note that, in addition, `numpy` defines the following members:\n",
"\n",
"- `.ndim`: the number of dimensions of the array (in our case, it would be 1, or 2)\n",
"- `.size`: the number of elements in the array; it is the product of m, and n\n",
"- `.dtype`: the data type; in our case, it is basically stored in data->typecode\n",
"- `.itemsize`: the size of a single element in the array: this can be gotten by calling `mp_binary_get_size('@', data->typecode, NULL)`.\n",
"\n",
"One should, perhaps, consider, whether these are necessary fields. E.g., if `ndim` were defined, then \n",
"\n",
"```c\n",
"if((myarray->m == 1) || (myarray->n == 1)) {\n",
" ...\n",
"}\n",
"```\n",
"\n",
"could be replaced by \n",
"\n",
"```c\n",
"if(myarray->ndim == 1) {\n",
" ...\n",
"}\n",
"```\n",
"and \n",
"```c\n",
"if((myarray->m > 1) && (myarray->n > 1)) {\n",
" ...\n",
"}\n",
"```\n",
"would be equivalent to \n",
"```c\n",
"if(myarray->ndim == 2) {\n",
" ...\n",
"}\n",
"```\n",
"\n",
"One could also save the extra function call `mp_binary_get_size('@', data->typecode, NULL)`, if `itemsize` is available. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ndarray.h"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T06:20:59.074722Z",
"start_time": "2020-02-26T06:20:59.063733Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 5469 bytes to ndarray.h\n"
]
}
],
"source": [
"%%ccode ndarray.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#ifndef _NDARRAY_\n",
"#define _NDARRAY_\n",
"\n",
"#include \"py/objarray.h\"\n",
"#include \"py/binary.h\"\n",
"#include \"py/objstr.h\"\n",
"#include \"py/objlist.h\"\n",
"\n",
"#define PRINT_MAX 10\n",
"\n",
"#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT\n",
"#define FLOAT_TYPECODE 'f'\n",
"#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE\n",
"#define FLOAT_TYPECODE 'd'\n",
"#endif\n",
"\n",
"#if !CIRCUITPY\n",
"#define translate(x) x\n",
"#endif\n",
"\n",
"#define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }\n",
"\n",
"extern const mp_obj_type_t ulab_ndarray_type;\n",
"\n",
"enum NDARRAY_TYPE {\n",
" NDARRAY_UINT8 = 'B',\n",
" NDARRAY_INT8 = 'b',\n",
" NDARRAY_UINT16 = 'H', \n",
" NDARRAY_INT16 = 'h',\n",
" NDARRAY_FLOAT = FLOAT_TYPECODE,\n",
"};\n",
"\n",
"typedef struct _ndarray_obj_t {\n",
" mp_obj_base_t base;\n",
" size_t m, n;\n",
" size_t len;\n",
" mp_obj_array_t *array;\n",
" size_t bytes;\n",
"} ndarray_obj_t;\n",
"\n",
"mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t , size_t , mp_obj_iter_buf_t *);\n",
"\n",
"mp_float_t ndarray_get_float_value(void *, uint8_t , size_t );\n",
"void fill_array_iterable(mp_float_t *, mp_obj_t );\n",
"\n",
"void ndarray_print_row(const mp_print_t *, mp_obj_array_t *, size_t , size_t );\n",
"void ndarray_print(const mp_print_t *, mp_obj_t , mp_print_kind_t );\n",
"void ndarray_assign_elements(mp_obj_array_t *, mp_obj_t , uint8_t , size_t *);\n",
"ndarray_obj_t *create_new_ndarray(size_t , size_t , uint8_t );\n",
"\n",
"mp_obj_t ndarray_copy(mp_obj_t );\n",
"#ifdef CIRCUITPY\n",
"mp_obj_t ndarray_make_new(const mp_obj_type_t *type, size_t n_args, const mp_obj_t *args, mp_map_t *kw_args);\n",
"#else\n",
"mp_obj_t ndarray_make_new(const mp_obj_type_t *, size_t , size_t , const mp_obj_t *);\n",
"#endif\n",
"mp_obj_t ndarray_subscr(mp_obj_t , mp_obj_t , mp_obj_t );\n",
"mp_obj_t ndarray_getiter(mp_obj_t , mp_obj_iter_buf_t *);\n",
"mp_obj_t ndarray_binary_op(mp_binary_op_t , mp_obj_t , mp_obj_t );\n",
"mp_obj_t ndarray_unary_op(mp_unary_op_t , mp_obj_t );\n",
"\n",
"mp_obj_t ndarray_shape(mp_obj_t );\n",
"mp_obj_t ndarray_size(mp_obj_t );\n",
"mp_obj_t ndarray_itemsize(mp_obj_t );\n",
"mp_obj_t ndarray_flatten(size_t , const mp_obj_t *, mp_map_t *);\n",
"\n",
"mp_obj_t ndarray_reshape(mp_obj_t , mp_obj_t );\n",
"MP_DECLARE_CONST_FUN_OBJ_2(ndarray_reshape_obj);\n",
"\n",
"mp_obj_t ndarray_transpose(mp_obj_t );\n",
"MP_DECLARE_CONST_FUN_OBJ_1(ndarray_transpose_obj);\n",
"\n",
"mp_int_t ndarray_get_buffer(mp_obj_t obj, mp_buffer_info_t *bufinfo, mp_uint_t flags);\n",
"//void ndarray_attributes(mp_obj_t , qstr , mp_obj_t *);\n",
"\n",
"\n",
"#define CREATE_SINGLE_ITEM(outarray, type, typecode, value) do {\\\n",
" ndarray_obj_t *tmp = create_new_ndarray(1, 1, (typecode));\\\n",
" type *tmparr = (type *)tmp->array->items;\\\n",
" tmparr[0] = (type)(value);\\\n",
" (outarray) = MP_OBJ_FROM_PTR(tmp);\\\n",
"} while(0)\n",
"\n",
"/* \n",
" mp_obj_t row = mp_obj_new_list(n, NULL);\n",
" mp_obj_list_t *row_ptr = MP_OBJ_TO_PTR(row);\n",
" \n",
" should work outside the loop, but it doesn't. Go figure! \n",
"*/\n",
"\n",
"#define RUN_BINARY_LOOP(typecode, type_out, type_left, type_right, ol, or, op) do {\\\n",
" type_left *left = (type_left *)(ol)->array->items;\\\n",
" type_right *right = (type_right *)(or)->array->items;\\\n",
" uint8_t inc = 0;\\\n",
" if((or)->array->len > 1) inc = 1;\\\n",
" if(((op) == MP_BINARY_OP_ADD) || ((op) == MP_BINARY_OP_SUBTRACT) || ((op) == MP_BINARY_OP_MULTIPLY)) {\\\n",
" ndarray_obj_t *out = create_new_ndarray(ol->m, ol->n, typecode);\\\n",
" type_out *(odata) = (type_out *)out->array->items;\\\n",
" 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];}\\\n",
" 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];}\\\n",
" 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];}\\\n",
" return MP_OBJ_FROM_PTR(out);\\\n",
" } else if((op) == MP_BINARY_OP_TRUE_DIVIDE) {\\\n",
" ndarray_obj_t *out = create_new_ndarray(ol->m, ol->n, NDARRAY_FLOAT);\\\n",
" mp_float_t *odata = (mp_float_t *)out->array->items;\\\n",
" for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = (mp_float_t)left[i]/(mp_float_t)right[j];\\\n",
" return MP_OBJ_FROM_PTR(out);\\\n",
" } else if(((op) == MP_BINARY_OP_LESS) || ((op) == MP_BINARY_OP_LESS_EQUAL) || \\\n",
" ((op) == MP_BINARY_OP_MORE) || ((op) == MP_BINARY_OP_MORE_EQUAL)) {\\\n",
" mp_obj_t out_list = mp_obj_new_list(0, NULL);\\\n",
" size_t m = (ol)->m, n = (ol)->n;\\\n",
" for(size_t i=0, r=0; i < m; i++, r+=inc) {\\\n",
" mp_obj_t row = mp_obj_new_list(n, NULL);\\\n",
" mp_obj_list_t *row_ptr = MP_OBJ_TO_PTR(row);\\\n",
" for(size_t j=0, s=0; j < n; j++, s+=inc) {\\\n",
" row_ptr->items[j] = mp_const_false;\\\n",
" if((op) == MP_BINARY_OP_LESS) {\\\n",
" if(left[i*n+j] < right[r*n+s]) row_ptr->items[j] = mp_const_true;\\\n",
" } else if((op) == MP_BINARY_OP_LESS_EQUAL) {\\\n",
" if(left[i*n+j] <= right[r*n+s]) row_ptr->items[j] = mp_const_true;\\\n",
" } else if((op) == MP_BINARY_OP_MORE) {\\\n",
" if(left[i*n+j] > right[r*n+s]) row_ptr->items[j] = mp_const_true;\\\n",
" } else if((op) == MP_BINARY_OP_MORE_EQUAL) {\\\n",
" if(left[i*n+j] >= right[r*n+s]) row_ptr->items[j] = mp_const_true;\\\n",
" }\\\n",
" }\\\n",
" if(m == 1) return row;\\\n",
" mp_obj_list_append(out_list, row);\\\n",
" }\\\n",
" return out_list;\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"#endif\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ndarray.c"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T06:15:40.809455Z",
"start_time": "2020-02-26T06:15:40.797803Z"
},
"code_folding": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 43934 bytes to ndarray.c\n"
]
}
],
"source": [
"%%ccode ndarray.c\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#include <math.h>\n",
"#include <stdio.h>\n",
"#include <stdlib.h>\n",
"#include <string.h>\n",
"#include \"py/runtime.h\"\n",
"#include \"py/binary.h\"\n",
"#include \"py/obj.h\"\n",
"#include \"py/objtuple.h\"\n",
"#include \"ndarray.h\"\n",
"\n",
"// This function is copied verbatim from objarray.c\n",
"STATIC mp_obj_array_t *array_new(char typecode, size_t n) {\n",
" int typecode_size = mp_binary_get_size('@', typecode, NULL);\n",
" mp_obj_array_t *o = m_new_obj(mp_obj_array_t);\n",
" // this step could probably be skipped: we are never going to store a bytearray per se\n",
" #if MICROPY_PY_BUILTINS_BYTEARRAY && MICROPY_PY_ARRAY\n",
" o->base.type = (typecode == BYTEARRAY_TYPECODE) ? &mp_type_bytearray : &mp_type_array;\n",
" #elif MICROPY_PY_BUILTINS_BYTEARRAY\n",
" o->base.type = &mp_type_bytearray;\n",
" #else\n",
" o->base.type = &mp_type_array;\n",
" #endif\n",
" o->typecode = typecode;\n",
" o->free = 0;\n",
" o->len = n;\n",
" o->items = m_new(byte, typecode_size * o->len);\n",
" return o;\n",
"}\n",
"\n",
"mp_float_t ndarray_get_float_value(void *data, uint8_t typecode, size_t index) {\n",
" if(typecode == NDARRAY_UINT8) {\n",
" return (mp_float_t)((uint8_t *)data)[index];\n",
" } else if(typecode == NDARRAY_INT8) {\n",
" return (mp_float_t)((int8_t *)data)[index];\n",
" } else if(typecode == NDARRAY_UINT16) {\n",
" return (mp_float_t)((uint16_t *)data)[index];\n",
" } else if(typecode == NDARRAY_INT16) {\n",
" return (mp_float_t)((int16_t *)data)[index];\n",
" } else {\n",
" return (mp_float_t)((mp_float_t *)data)[index];\n",
" }\n",
"}\n",
"\n",
"void fill_array_iterable(mp_float_t *array, mp_obj_t iterable) {\n",
" mp_obj_iter_buf_t x_buf;\n",
" mp_obj_t x_item, x_iterable = mp_getiter(iterable, &x_buf);\n",
" size_t i=0;\n",
" while ((x_item = mp_iternext(x_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" array[i] = (mp_float_t)mp_obj_get_float(x_item);\n",
" i++;\n",
" }\n",
"}\n",
"\n",
"void ndarray_print_row(const mp_print_t *print, mp_obj_array_t *data, size_t n0, size_t n) {\n",
" mp_print_str(print, \"[\");\n",
" size_t i;\n",
" if(n < PRINT_MAX) { // if the array is short, print everything\n",
" mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0), PRINT_REPR);\n",
" for(i=1; i<n; i++) {\n",
" mp_print_str(print, \", \");\n",
" mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+i), PRINT_REPR);\n",
" }\n",
" } else {\n",
" mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0), PRINT_REPR);\n",
" for(i=1; i<3; i++) {\n",
" mp_print_str(print, \", \");\n",
" mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+i), PRINT_REPR);\n",
" }\n",
" mp_printf(print, \", ..., \");\n",
" mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+n-3), PRINT_REPR);\n",
" for(size_t i=1; i<3; i++) {\n",
" mp_print_str(print, \", \");\n",
" mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+n-3+i), PRINT_REPR);\n",
" }\n",
" }\n",
" mp_print_str(print, \"]\");\n",
"}\n",
"\n",
"void ndarray_print(const mp_print_t *print, mp_obj_t self_in, mp_print_kind_t kind) {\n",
" (void)kind;\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" mp_print_str(print, \"array(\");\n",
" \n",
" if(self->array->len == 0) {\n",
" mp_print_str(print, \"[]\");\n",
" } else {\n",
" if((self->m == 1) || (self->n == 1)) {\n",
" ndarray_print_row(print, self->array, 0, self->array->len);\n",
" } else {\n",
" // TODO: add vertical ellipses for the case, when self->m > PRINT_MAX\n",
" mp_print_str(print, \"[\");\n",
" ndarray_print_row(print, self->array, 0, self->n);\n",
" for(size_t i=1; i < self->m; i++) {\n",
" mp_print_str(print, \",\\n\\t \");\n",
" ndarray_print_row(print, self->array, i*self->n, self->n);\n",
" }\n",
" mp_print_str(print, \"]\");\n",
" }\n",
" }\n",
" if(self->array->typecode == NDARRAY_UINT8) {\n",
" mp_print_str(print, \", dtype=uint8)\");\n",
" } else if(self->array->typecode == NDARRAY_INT8) {\n",
" mp_print_str(print, \", dtype=int8)\");\n",
" } else if(self->array->typecode == NDARRAY_UINT16) {\n",
" mp_print_str(print, \", dtype=uint16)\");\n",
" } else if(self->array->typecode == NDARRAY_INT16) {\n",
" mp_print_str(print, \", dtype=int16)\");\n",
" } else if(self->array->typecode == NDARRAY_FLOAT) {\n",
" mp_print_str(print, \", dtype=float)\");\n",
" }\n",
"}\n",
"\n",
"void ndarray_assign_elements(mp_obj_array_t *data, mp_obj_t iterable, uint8_t typecode, size_t *idx) {\n",
" // assigns a single row in the matrix\n",
" mp_obj_t item;\n",
" while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" mp_binary_set_val_array(typecode, data->items, (*idx)++, item);\n",
" }\n",
"}\n",
"\n",
"ndarray_obj_t *create_new_ndarray(size_t m, size_t n, uint8_t typecode) {\n",
" // Creates the base ndarray with shape (m, n), and initialises the values to straight 0s\n",
" ndarray_obj_t *ndarray = m_new_obj(ndarray_obj_t);\n",
" ndarray->base.type = &ulab_ndarray_type;\n",
" ndarray->m = m;\n",
" ndarray->n = n;\n",
" mp_obj_array_t *array = array_new(typecode, m*n);\n",
" ndarray->bytes = m * n * mp_binary_get_size('@', typecode, NULL);\n",
" // this should set all elements to 0, irrespective of the of the typecode (all bits are zero)\n",
" // we could, perhaps, leave this step out, and initialise the array only, when needed\n",
" memset(array->items, 0, ndarray->bytes); \n",
" ndarray->array = array;\n",
" return ndarray;\n",
"}\n",
"\n",
"mp_obj_t ndarray_copy(mp_obj_t self_in) {\n",
" // returns a verbatim (shape and typecode) copy of self_in\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" ndarray_obj_t *out = create_new_ndarray(self->m, self->n, self->array->typecode);\n",
" memcpy(out->array->items, self->array->items, self->bytes);\n",
" return MP_OBJ_FROM_PTR(out);\n",
"}\n",
"\n",
"STATIC uint8_t ndarray_init_helper(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT } },\n",
" };\n",
" \n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" \n",
" uint8_t dtype = args[1].u_int;\n",
" return dtype;\n",
"}\n",
"\n",
"STATIC mp_obj_t ndarray_make_new_core(const mp_obj_type_t *type, size_t n_args, size_t n_kw, const mp_obj_t *args, mp_map_t *kw_args) {\n",
" uint8_t dtype = ndarray_init_helper(n_args, args, kw_args);\n",
"\n",
" size_t len1, len2=0, i=0;\n",
" mp_obj_t len_in = mp_obj_len_maybe(args[0]);\n",
" if (len_in == MP_OBJ_NULL) {\n",
" mp_raise_ValueError(translate(\"first argument must be an iterable\"));\n",
" } else {\n",
" // len1 is either the number of rows (for matrices), or the number of elements (row vectors)\n",
" len1 = MP_OBJ_SMALL_INT_VALUE(len_in);\n",
" }\n",
"\n",
" // We have to figure out, whether the first element of the iterable is an iterable itself\n",
" // Perhaps, there is a more elegant way of handling this\n",
" mp_obj_iter_buf_t iter_buf1;\n",
" mp_obj_t item1, iterable1 = mp_getiter(args[0], &iter_buf1);\n",
" while ((item1 = mp_iternext(iterable1)) != MP_OBJ_STOP_ITERATION) {\n",
" len_in = mp_obj_len_maybe(item1);\n",
" if(len_in != MP_OBJ_NULL) { // indeed, this seems to be an iterable\n",
" // Next, we have to check, whether all elements in the outer loop have the same length\n",
" if(i > 0) {\n",
" if(len2 != MP_OBJ_SMALL_INT_VALUE(len_in)) {\n",
" mp_raise_ValueError(translate(\"iterables are not of the same length\"));\n",
" }\n",
" }\n",
" len2 = MP_OBJ_SMALL_INT_VALUE(len_in);\n",
" i++;\n",
" }\n",
" }\n",
" // By this time, it should be established, what the shape is, so we can now create the array\n",
" ndarray_obj_t *self = create_new_ndarray((len2 == 0) ? 1 : len1, (len2 == 0) ? len1 : len2, dtype);\n",
" iterable1 = mp_getiter(args[0], &iter_buf1);\n",
" i = 0;\n",
" if(len2 == 0) { // the first argument is a single iterable\n",
" ndarray_assign_elements(self->array, iterable1, dtype, &i);\n",
" } else {\n",
" mp_obj_iter_buf_t iter_buf2;\n",
" mp_obj_t iterable2; \n",
"\n",
" while ((item1 = mp_iternext(iterable1)) != MP_OBJ_STOP_ITERATION) {\n",
" iterable2 = mp_getiter(item1, &iter_buf2);\n",
" ndarray_assign_elements(self->array, iterable2, dtype, &i);\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(self);\n",
"}\n",
"\n",
"#ifdef CIRCUITPY\n",
"mp_obj_t ndarray_make_new(const mp_obj_type_t *type, size_t n_args, const mp_obj_t *args, mp_map_t *kw_args) {\n",
" mp_arg_check_num(n_args, kw_args, 1, 2, true);\n",
" size_t n_kw = 0;\n",
" if (kw_args != 0) {\n",
" n_kw = kw_args->used;\n",
" }\n",
" mp_map_init_fixed_table(kw_args, n_kw, args + n_args);\n",
" return ndarray_make_new_core(type, n_args, n_kw, args, kw_args);\n",
"}\n",
"#else\n",
"mp_obj_t ndarray_make_new(const mp_obj_type_t *type, size_t n_args, size_t n_kw, const mp_obj_t *args) {\n",
" mp_arg_check_num(n_args, n_kw, 1, 2, true);\n",
" mp_map_t kw_args;\n",
" mp_map_init_fixed_table(&kw_args, n_kw, args + n_args);\n",
" return ndarray_make_new_core(type, n_args, n_kw, args, &kw_args);\n",
"}\n",
"#endif\n",
"\n",
"size_t slice_length(mp_bound_slice_t slice) {\n",
" int32_t len, correction = 1;\n",
" if(slice.step > 0) correction = -1;\n",
" len = (slice.stop - slice.start + (slice.step + correction)) / slice.step;\n",
" if(len < 0) return 0;\n",
" return (size_t)len;\n",
"}\n",
"\n",
"size_t true_length(mp_obj_t bool_list) {\n",
" // returns the number of Trues in a Boolean list\n",
" // I wonder, wouldn't this be faster, if we looped through bool_list->items instead?\n",
" mp_obj_iter_buf_t iter_buf;\n",
" mp_obj_t item, iterable = mp_getiter(bool_list, &iter_buf);\n",
" size_t trues = 0;\n",
" while((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if(!MP_OBJ_IS_TYPE(item, &mp_type_bool)) {\n",
" // numpy seems to be a little bit inconsistent in when an index is considered\n",
" // to be True/False. Bail out immediately, if the items are not True/False\n",
" return 0;\n",
" }\n",
" if(mp_obj_is_true(item)) {\n",
" trues++;\n",
" }\n",
" }\n",
" return trues;\n",
"}\n",
"\n",
"mp_bound_slice_t generate_slice(mp_uint_t n, mp_obj_t index) {\n",
" // micropython seems to have difficulties with negative steps\n",
" mp_bound_slice_t slice;\n",
" if(MP_OBJ_IS_TYPE(index, &mp_type_slice)) {\n",
" mp_seq_get_fast_slice_indexes(n, index, &slice);\n",
" } else if(MP_OBJ_IS_INT(index)) {\n",
" int32_t _index = mp_obj_get_int(index);\n",
" if(_index < 0) {\n",
" _index += n;\n",
" } \n",
" if((_index >= n) || (_index < 0)) {\n",
" mp_raise_msg(&mp_type_IndexError, translate(\"index is out of bounds\"));\n",
" }\n",
" slice.start = _index;\n",
" slice.stop = _index + 1;\n",
" slice.step = 1;\n",
" } else {\n",
" mp_raise_msg(&mp_type_IndexError, translate(\"indices must be integers, slices, or Boolean lists\"));\n",
" }\n",
" return slice;\n",
"}\n",
"\n",
"mp_bound_slice_t simple_slice(int16_t start, int16_t stop, int16_t step) {\n",
" mp_bound_slice_t slice;\n",
" slice.start = start;\n",
" slice.stop = stop;\n",
" slice.step = step;\n",
" return slice;\n",
"}\n",
"\n",
"void insert_binary_value(ndarray_obj_t *ndarray, size_t nd_index, ndarray_obj_t *values, size_t value_index) {\n",
" // there is probably a more elegant implementation...\n",
" mp_obj_t tmp = mp_binary_get_val_array(values->array->typecode, values->array->items, value_index);\n",
" if((values->array->typecode == NDARRAY_FLOAT) && (ndarray->array->typecode != NDARRAY_FLOAT)) {\n",
" // workaround: rounding seems not to work in the arm compiler\n",
" int32_t x = (int32_t)floorf(mp_obj_get_float(tmp)+0.5);\n",
" tmp = mp_obj_new_int(x);\n",
" }\n",
" mp_binary_set_val_array(ndarray->array->typecode, ndarray->array->items, nd_index, tmp); \n",
"}\n",
"\n",
"mp_obj_t insert_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n, \n",
" mp_bound_slice_t row, mp_bound_slice_t column, \n",
" mp_obj_t row_list, mp_obj_t column_list, \n",
" ndarray_obj_t *values) {\n",
" if((m != values->m) && (n != values->n)) {\n",
" if(values->array->len != 1) { // not a single item\n",
" mp_raise_ValueError(translate(\"could not broadast input array from shape\"));\n",
" }\n",
" }\n",
" size_t cindex, rindex;\n",
" // M, and N are used to manipulate how the source index is incremented in the loop\n",
" uint8_t M = 1, N = 1;\n",
" if(values->m == 1) {\n",
" M = 0;\n",
" }\n",
" if(values->n == 1) {\n",
" N = 0;\n",
" }\n",
" \n",
" if(row_list == mp_const_none) { // rows are indexed by a slice\n",
" rindex = row.start;\n",
" if(column_list == mp_const_none) { // columns are indexed by a slice\n",
" for(size_t i=0; i < m; i++) {\n",
" cindex = column.start;\n",
" for(size_t j=0; j < n; j++) {\n",
" insert_binary_value(ndarray, rindex*ndarray->n+cindex, values, i*M*n+j*N);\n",
" cindex += column.step;\n",
" }\n",
" rindex += row.step;\n",
" }\n",
" } else { // columns are indexed by a Boolean list\n",
" mp_obj_iter_buf_t column_iter_buf;\n",
" mp_obj_t column_item, column_iterable;\n",
" for(size_t i=0; i < m; i++) {\n",
" column_iterable = mp_getiter(column_list, &column_iter_buf);\n",
" size_t j = 0;\n",
" cindex = 0;\n",
" while((column_item = mp_iternext(column_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if(mp_obj_is_true(column_item)) {\n",
" insert_binary_value(ndarray, rindex*ndarray->n+cindex, values, i*M*n+j*N);\n",
" j++;\n",
" }\n",
" cindex++;\n",
" }\n",
" rindex += row.step;\n",
" }\n",
" }\n",
" } else { // rows are indexed by a Boolean list\n",
" mp_obj_iter_buf_t row_iter_buf;\n",
" mp_obj_t row_item, row_iterable;\n",
" row_iterable = mp_getiter(row_list, &row_iter_buf);\n",
" size_t i = 0;\n",
" rindex = 0;\n",
" if(column_list == mp_const_none) { // columns are indexed by a slice\n",
" while((row_item = mp_iternext(row_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if(mp_obj_is_true(row_item)) {\n",
" cindex = column.start;\n",
" for(size_t j=0; j < n; j++) {\n",
" insert_binary_value(ndarray, rindex*ndarray->n+cindex, values, i*M*n+j*N);\n",
" cindex += column.step;\n",
" }\n",
" i++;\n",
" }\n",
" rindex++;\n",
" } \n",
" } else { // columns are indexed by a list\n",
" mp_obj_iter_buf_t column_iter_buf;\n",
" mp_obj_t column_item, column_iterable;\n",
" size_t j = 0, cindex = 0;\n",
" while((row_item = mp_iternext(row_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if(mp_obj_is_true(row_item)) {\n",
" column_iterable = mp_getiter(column_list, &column_iter_buf); \n",
" while((column_item = mp_iternext(column_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if(mp_obj_is_true(column_item)) {\n",
" insert_binary_value(ndarray, rindex*ndarray->n+cindex, values, i*M*n+j*N);\n",
" j++;\n",
" }\n",
" cindex++;\n",
" }\n",
" i++;\n",
" }\n",
" rindex++;\n",
" }\n",
" }\n",
" }\n",
" return mp_const_none;\n",
"}\n",
"\n",
"mp_obj_t iterate_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n, \n",
" mp_bound_slice_t row, mp_bound_slice_t column, \n",
" mp_obj_t row_list, mp_obj_t column_list, \n",
" ndarray_obj_t *values) {\n",
" if((m == 0) || (n == 0)) {\n",
" mp_raise_msg(&mp_type_IndexError, translate(\"empty index range\"));\n",
" }\n",
"\n",
" if(values != NULL) {\n",
" return insert_slice_list(ndarray, m, n, row, column, row_list, column_list, values);\n",
" }\n",
" uint8_t _sizeof = mp_binary_get_size('@', ndarray->array->typecode, NULL);\n",
" ndarray_obj_t *out = create_new_ndarray(m, n, ndarray->array->typecode);\n",
" uint8_t *target = (uint8_t *)out->array->items;\n",
" uint8_t *source = (uint8_t *)ndarray->array->items;\n",
" size_t cindex, rindex; \n",
" if(row_list == mp_const_none) { // rows are indexed by a slice\n",
" rindex = row.start;\n",
" if(column_list == mp_const_none) { // columns are indexed by a slice\n",
" for(size_t i=0; i < m; i++) {\n",
" cindex = column.start;\n",
" for(size_t j=0; j < n; j++) {\n",
" memcpy(target+(i*n+j)*_sizeof, source+(rindex*ndarray->n+cindex)*_sizeof, _sizeof);\n",
" cindex += column.step;\n",
" }\n",
" rindex += row.step;\n",
" }\n",
" } else { // columns are indexed by a Boolean list\n",
" // TODO: the list must be exactly as long as the axis\n",
" mp_obj_iter_buf_t column_iter_buf;\n",
" mp_obj_t column_item, column_iterable;\n",
" for(size_t i=0; i < m; i++) {\n",
" column_iterable = mp_getiter(column_list, &column_iter_buf);\n",
" size_t j = 0;\n",
" cindex = 0;\n",
" while((column_item = mp_iternext(column_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if(mp_obj_is_true(column_item)) {\n",
" memcpy(target+(i*n+j)*_sizeof, source+(rindex*ndarray->n+cindex)*_sizeof, _sizeof);\n",
" j++;\n",
" }\n",
" cindex++;\n",
" }\n",
" rindex += row.step;\n",
" }\n",
" }\n",
" } else { // rows are indexed by a Boolean list\n",
" mp_obj_iter_buf_t row_iter_buf;\n",
" mp_obj_t row_item, row_iterable;\n",
" row_iterable = mp_getiter(row_list, &row_iter_buf);\n",
" size_t i = 0;\n",
" rindex = 0;\n",
" if(column_list == mp_const_none) { // columns are indexed by a slice\n",
" while((row_item = mp_iternext(row_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if(mp_obj_is_true(row_item)) {\n",
" cindex = column.start;\n",
" for(size_t j=0; j < n; j++) {\n",
" memcpy(target+(i*n+j)*_sizeof, source+(rindex*ndarray->n+cindex)*_sizeof, _sizeof);\n",
" cindex += column.step;\n",
" }\n",
" i++;\n",
" }\n",
" rindex++;\n",
" } \n",
" } else { // columns are indexed by a list\n",
" mp_obj_iter_buf_t column_iter_buf;\n",
" mp_obj_t column_item, column_iterable;\n",
" size_t j = 0, cindex = 0;\n",
" while((row_item = mp_iternext(row_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if(mp_obj_is_true(row_item)) {\n",
" column_iterable = mp_getiter(column_list, &column_iter_buf); \n",
" while((column_item = mp_iternext(column_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if(mp_obj_is_true(column_item)) {\n",
" memcpy(target+(i*n+j)*_sizeof, source+(rindex*ndarray->n+cindex)*_sizeof, _sizeof);\n",
" j++;\n",
" }\n",
" cindex++;\n",
" }\n",
" i++;\n",
" }\n",
" rindex++;\n",
" }\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(out);\n",
"}\n",
"\n",
"mp_obj_t ndarray_get_slice(ndarray_obj_t *ndarray, mp_obj_t index, ndarray_obj_t *values) {\n",
" mp_bound_slice_t row_slice = simple_slice(0, 0, 1), column_slice = simple_slice(0, 0, 1);\n",
"\n",
" size_t m = 0, n = 0;\n",
" if(MP_OBJ_IS_INT(index) && (ndarray->m == 1) && (values == NULL)) { \n",
" // we have a row vector, and don't want to assign\n",
" column_slice = generate_slice(ndarray->n, index);\n",
" if(slice_length(column_slice) == 1) { // we were asked for a single item\n",
" // subscribe returns an mp_obj_t, if and only, if the index is an integer, and we have a row vector\n",
" return mp_binary_get_val_array(ndarray->array->typecode, ndarray->array->items, column_slice.start);\n",
" }\n",
" }\n",
" \n",
" if(MP_OBJ_IS_INT(index) || MP_OBJ_IS_TYPE(index, &mp_type_slice)) {\n",
" if(ndarray->m == 1) { // we have a row vector\n",
" column_slice = generate_slice(ndarray->n, index);\n",
" row_slice = simple_slice(0, 1, 1);\n",
" } else { // we have a matrix\n",
" row_slice = generate_slice(ndarray->m, index);\n",
" column_slice = simple_slice(0, ndarray->n, 1); // take all columns\n",
" }\n",
" m = slice_length(row_slice);\n",
" n = slice_length(column_slice);\n",
" return iterate_slice_list(ndarray, m, n, row_slice, column_slice, mp_const_none, mp_const_none, values);\n",
" } else if(MP_OBJ_IS_TYPE(index, &mp_type_list)) {\n",
" n = true_length(index);\n",
" if(ndarray->m == 1) { // we have a flat array\n",
" // we might have to separate the n == 1 case\n",
" row_slice = simple_slice(0, 1, 1);\n",
" return iterate_slice_list(ndarray, 1, n, row_slice, column_slice, mp_const_none, index, values);\n",
" } else { // we have a matrix\n",
" return iterate_slice_list(ndarray, 1, n, row_slice, column_slice, mp_const_none, index, values);\n",
" }\n",
" }\n",
" else { // we certainly have a tuple, so let us deal with it\n",
" mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(index);\n",
" if(tuple->len != 2) {\n",
" mp_raise_msg(&mp_type_IndexError, translate(\"too many indices\"));\n",
" }\n",
" if(!(MP_OBJ_IS_TYPE(tuple->items[0], &mp_type_list) || \n",
" MP_OBJ_IS_TYPE(tuple->items[0], &mp_type_slice) || \n",
" MP_OBJ_IS_INT(tuple->items[0])) || \n",
" !(MP_OBJ_IS_TYPE(tuple->items[1], &mp_type_list) || \n",
" MP_OBJ_IS_TYPE(tuple->items[1], &mp_type_slice) || \n",
" MP_OBJ_IS_INT(tuple->items[1]))) {\n",
" mp_raise_msg(&mp_type_IndexError, translate(\"indices must be integers, slices, or Boolean lists\"));\n",
" }\n",
" if(MP_OBJ_IS_TYPE(tuple->items[0], &mp_type_list)) { // rows are indexed by Boolean list\n",
" m = true_length(tuple->items[0]);\n",
" if(MP_OBJ_IS_TYPE(tuple->items[1], &mp_type_list)) {\n",
" n = true_length(tuple->items[1]);\n",
" return iterate_slice_list(ndarray, m, n, row_slice, column_slice, \n",
" tuple->items[0], tuple->items[1], values);\n",
" } else { // the column is indexed by an integer, or a slice\n",
" column_slice = generate_slice(ndarray->n, tuple->items[1]);\n",
" n = slice_length(column_slice);\n",
" return iterate_slice_list(ndarray, m, n, row_slice, column_slice, \n",
" tuple->items[0], mp_const_none, values);\n",
" }\n",
" \n",
" } else { // rows are indexed by a slice, or an integer\n",
" row_slice = generate_slice(ndarray->m, tuple->items[0]);\n",
" m = slice_length(row_slice);\n",
" if(MP_OBJ_IS_TYPE(tuple->items[1], &mp_type_list)) { // columns are indexed by a Boolean list\n",
" n = true_length(tuple->items[1]);\n",
" return iterate_slice_list(ndarray, m, n, row_slice, column_slice, \n",
" mp_const_none, tuple->items[1], values);\n",
" } else { // columns are indexed by an integer, or a slice\n",
" column_slice = generate_slice(ndarray->n, tuple->items[1]);\n",
" n = slice_length(column_slice);\n",
" return iterate_slice_list(ndarray, m, n, row_slice, column_slice, \n",
" mp_const_none, mp_const_none, values); \n",
" \n",
" }\n",
" }\n",
" }\n",
"}\n",
"\n",
"mp_obj_t ndarray_subscr(mp_obj_t self_in, mp_obj_t index, mp_obj_t value) {\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" \n",
" if (value == MP_OBJ_SENTINEL) { // return value(s)\n",
" return ndarray_get_slice(self, index, NULL); \n",
" } else { // assignment to slices; the value must be an ndarray, or a scalar\n",
" if(!MP_OBJ_IS_TYPE(value, &ulab_ndarray_type) && \n",
" !MP_OBJ_IS_INT(value) && !mp_obj_is_float(value)) {\n",
" mp_raise_ValueError(translate(\"right hand side must be an ndarray, or a scalar\"));\n",
" } else {\n",
" ndarray_obj_t *values = NULL;\n",
" if(MP_OBJ_IS_INT(value)) {\n",
" values = create_new_ndarray(1, 1, self->array->typecode);\n",
" mp_binary_set_val_array(values->array->typecode, values->array->items, 0, value); \n",
" } else if(mp_obj_is_float(value)) {\n",
" values = create_new_ndarray(1, 1, NDARRAY_FLOAT);\n",
" mp_binary_set_val_array(NDARRAY_FLOAT, values->array->items, 0, value);\n",
" } else {\n",
" values = MP_OBJ_TO_PTR(value);\n",
" }\n",
" return ndarray_get_slice(self, index, values);\n",
" }\n",
" } \n",
" return mp_const_none;\n",
"}\n",
"\n",
"// itarray iterator\n",
"\n",
"mp_obj_t ndarray_getiter(mp_obj_t o_in, mp_obj_iter_buf_t *iter_buf) {\n",
" return mp_obj_new_ndarray_iterator(o_in, 0, iter_buf);\n",
"}\n",
"\n",
"typedef struct _mp_obj_ndarray_it_t {\n",
" mp_obj_base_t base;\n",
" mp_fun_1_t iternext;\n",
" mp_obj_t ndarray;\n",
" size_t cur;\n",
"} mp_obj_ndarray_it_t;\n",
"\n",
"mp_obj_t ndarray_iternext(mp_obj_t self_in) {\n",
" mp_obj_ndarray_it_t *self = MP_OBJ_TO_PTR(self_in);\n",
" ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(self->ndarray);\n",
" // TODO: in numpy, ndarrays are iterated with respect to the first axis. \n",
" size_t iter_end = 0;\n",
" if(ndarray->m == 1) {\n",
" iter_end = ndarray->array->len;\n",
" } else {\n",
" iter_end = ndarray->m;\n",
" }\n",
" if(self->cur < iter_end) {\n",
" if(ndarray->n == ndarray->array->len) { // we have a linear array\n",
" // read the current value\n",
" mp_obj_t value;\n",
" value = mp_binary_get_val_array(ndarray->array->typecode, ndarray->array->items, self->cur);\n",
" self->cur++;\n",
" return value;\n",
" } else { // we have a matrix, return the number of rows\n",
" ndarray_obj_t *value = create_new_ndarray(1, ndarray->n, ndarray->array->typecode);\n",
" // copy the memory content here\n",
" uint8_t *tmp = (uint8_t *)ndarray->array->items;\n",
" size_t strip_size = ndarray->n * mp_binary_get_size('@', ndarray->array->typecode, NULL);\n",
" memcpy(value->array->items, &tmp[self->cur*strip_size], strip_size);\n",
" self->cur++;\n",
" return value;\n",
" }\n",
" } else {\n",
" return MP_OBJ_STOP_ITERATION;\n",
" }\n",
"}\n",
"\n",
"mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t ndarray, size_t cur, mp_obj_iter_buf_t *iter_buf) {\n",
" assert(sizeof(mp_obj_ndarray_it_t) <= sizeof(mp_obj_iter_buf_t));\n",
" mp_obj_ndarray_it_t *o = (mp_obj_ndarray_it_t*)iter_buf;\n",
" o->base.type = &mp_type_polymorph_iter;\n",
" o->iternext = ndarray_iternext;\n",
" o->ndarray = ndarray;\n",
" o->cur = cur;\n",
" return MP_OBJ_FROM_PTR(o);\n",
"}\n",
"\n",
"mp_obj_t ndarray_shape(mp_obj_t self_in) {\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" mp_obj_t tuple[2] = {\n",
" mp_obj_new_int(self->m),\n",
" mp_obj_new_int(self->n)\n",
" };\n",
" return mp_obj_new_tuple(2, tuple);\n",
"}\n",
"\n",
"mp_obj_t ndarray_size(mp_obj_t self_in) {\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" return mp_obj_new_int(self->array->len);\n",
"}\n",
"\n",
"mp_obj_t ndarray_itemsize(mp_obj_t self_in) {\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" return MP_OBJ_NEW_SMALL_INT(mp_binary_get_size('@', self->array->typecode, NULL));\n",
"}\n",
"\n",
"mp_obj_t ndarray_flatten(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_order, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_QSTR(MP_QSTR_C)} },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(n_args - 1, pos_args + 1, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" mp_obj_t self_copy = ndarray_copy(pos_args[0]);\n",
" ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(self_copy);\n",
" \n",
" GET_STR_DATA_LEN(args[0].u_obj, order, len); \n",
" if((len != 1) || ((memcmp(order, \"C\", 1) != 0) && (memcmp(order, \"F\", 1) != 0))) {\n",
" mp_raise_ValueError(translate(\"flattening order must be either 'C', or 'F'\"));\n",
" }\n",
"\n",
" // if order == 'C', we simply have to set m, and n, there is nothing else to do\n",
" if(memcmp(order, \"F\", 1) == 0) {\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(pos_args[0]);\n",
" uint8_t _sizeof = mp_binary_get_size('@', self->array->typecode, NULL);\n",
" // get the data of self_in: we won't need a temporary buffer for the transposition\n",
" uint8_t *self_array = (uint8_t *)self->array->items;\n",
" uint8_t *array = (uint8_t *)ndarray->array->items;\n",
" size_t i=0;\n",
" for(size_t n=0; n < self->n; n++) {\n",
" for(size_t m=0; m < self->m; m++) {\n",
" memcpy(array+_sizeof*i, self_array+_sizeof*(m*self->n + n), _sizeof);\n",
" i++;\n",
" }\n",
" } \n",
" }\n",
" ndarray->n = ndarray->array->len;\n",
" ndarray->m = 1;\n",
" return self_copy;\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(ndarray_flatten_obj, 1, ndarray_flatten);\n",
"\n",
"// Binary operations\n",
"\n",
"mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) {\n",
"// if(op == MP_BINARY_OP_REVERSE_ADD) {\n",
" // return ndarray_binary_op(MP_BINARY_OP_ADD, rhs, lhs);\n",
" // } \n",
" // One of the operands is a scalar\n",
" // TODO: conform to numpy with the upcasting\n",
" // TODO: implement in-place operators\n",
" mp_obj_t RHS = MP_OBJ_NULL;\n",
" bool rhs_is_scalar = true;\n",
" if(MP_OBJ_IS_INT(rhs)) {\n",
" int32_t ivalue = mp_obj_get_int(rhs);\n",
" if((ivalue > 0) && (ivalue < 256)) {\n",
" CREATE_SINGLE_ITEM(RHS, uint8_t, NDARRAY_UINT8, ivalue);\n",
" } else if((ivalue > 255) && (ivalue < 65535)) {\n",
" CREATE_SINGLE_ITEM(RHS, uint16_t, NDARRAY_UINT16, ivalue);\n",
" } else if((ivalue < 0) && (ivalue > -128)) {\n",
" CREATE_SINGLE_ITEM(RHS, int8_t, NDARRAY_INT8, ivalue);\n",
" } else if((ivalue < -127) && (ivalue > -32767)) {\n",
" CREATE_SINGLE_ITEM(RHS, int16_t, NDARRAY_INT16, ivalue);\n",
" } else { // the integer value clearly does not fit the ulab types, so move on to float\n",
" CREATE_SINGLE_ITEM(RHS, mp_float_t, NDARRAY_FLOAT, ivalue);\n",
" }\n",
" } else if(mp_obj_is_float(rhs)) {\n",
" mp_float_t fvalue = mp_obj_get_float(rhs); \n",
" CREATE_SINGLE_ITEM(RHS, mp_float_t, NDARRAY_FLOAT, fvalue);\n",
" } else {\n",
" RHS = rhs;\n",
" rhs_is_scalar = false;\n",
" }\n",
" //else \n",
" if(MP_OBJ_IS_TYPE(lhs, &ulab_ndarray_type) && MP_OBJ_IS_TYPE(RHS, &ulab_ndarray_type)) { \n",
" // next, the ndarray stuff\n",
" ndarray_obj_t *ol = MP_OBJ_TO_PTR(lhs);\n",
" ndarray_obj_t *or = MP_OBJ_TO_PTR(RHS);\n",
" if(!rhs_is_scalar && ((ol->m != or->m) || (ol->n != or->n))) {\n",
" mp_raise_ValueError(translate(\"operands could not be broadcast together\"));\n",
" }\n",
" // At this point, the operands should have the same shape\n",
" switch(op) {\n",
" case MP_BINARY_OP_EQUAL:\n",
" // Two arrays are equal, if their shape, typecode, and elements are equal\n",
" if((ol->m != or->m) || (ol->n != or->n) || (ol->array->typecode != or->array->typecode)) {\n",
" return mp_const_false;\n",
" } else {\n",
" size_t i = ol->bytes;\n",
" uint8_t *l = (uint8_t *)ol->array->items;\n",
" uint8_t *r = (uint8_t *)or->array->items;\n",
" while(i) { // At this point, we can simply compare the bytes, the type is irrelevant\n",
" if(*l++ != *r++) {\n",
" return mp_const_false;\n",
" }\n",
" i--;\n",
" }\n",
" return mp_const_true;\n",
" }\n",
" break;\n",
" case MP_BINARY_OP_LESS:\n",
" case MP_BINARY_OP_LESS_EQUAL:\n",
" case MP_BINARY_OP_MORE:\n",
" case MP_BINARY_OP_MORE_EQUAL:\n",
" case MP_BINARY_OP_ADD:\n",
" case MP_BINARY_OP_SUBTRACT:\n",
" case MP_BINARY_OP_TRUE_DIVIDE:\n",
" case MP_BINARY_OP_MULTIPLY:\n",
" // TODO: I believe, this part can be made significantly smaller (compiled size)\n",
" // by doing only the typecasting in the large ifs, and moving the loops outside\n",
" // These are the upcasting rules\n",
" // float always becomes float\n",
" // operation on identical types preserves type\n",
" // uint8 + int8 => int16\n",
" // uint8 + int16 => int16\n",
" // uint8 + uint16 => uint16\n",
" // int8 + int16 => int16\n",
" // int8 + uint16 => uint16\n",
" // uint16 + int16 => float\n",
" // The parameters of RUN_BINARY_LOOP are \n",
" // typecode of result, type_out, type_left, type_right, lhs operand, rhs operand, operator\n",
" if(ol->array->typecode == NDARRAY_UINT8) {\n",
" if(or->array->typecode == NDARRAY_UINT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_UINT8, uint8_t, uint8_t, uint8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_UINT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint8_t, uint16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ol, or, op);\n",
" }\n",
" } else if(ol->array->typecode == NDARRAY_INT8) {\n",
" if(or->array->typecode == NDARRAY_UINT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT8, int8_t, int8_t, int8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_UINT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int8_t, mp_float_t, ol, or, op);\n",
" } \n",
" } else if(ol->array->typecode == NDARRAY_UINT16) {\n",
" if(or->array->typecode == NDARRAY_UINT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, int8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_UINT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ol, or, op);\n",
" }\n",
" } else if(ol->array->typecode == NDARRAY_INT16) {\n",
" if(or->array->typecode == NDARRAY_UINT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, uint8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_UINT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int16_t, uint16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, mp_float_t, ol, or, op);\n",
" }\n",
" } else if(ol->array->typecode == NDARRAY_FLOAT) {\n",
" if(or->array->typecode == NDARRAY_UINT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT8) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int8_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_UINT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_INT16) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int16_t, ol, or, op);\n",
" } else if(or->array->typecode == NDARRAY_FLOAT) {\n",
" RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, mp_float_t, ol, or, op);\n",
" }\n",
" } else { // this should never happen\n",
" mp_raise_TypeError(translate(\"wrong input type\"));\n",
" }\n",
" // this instruction should never be reached, but we have to make the compiler happy\n",
" return MP_OBJ_NULL; \n",
" default:\n",
" return MP_OBJ_NULL; // op not supported \n",
" }\n",
" } else {\n",
" mp_raise_TypeError(translate(\"wrong operand type on the right hand side\"));\n",
" }\n",
"}\n",
"\n",
"mp_obj_t ndarray_unary_op(mp_unary_op_t op, mp_obj_t self_in) {\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" ndarray_obj_t *ndarray = NULL;\n",
" switch (op) {\n",
" case MP_UNARY_OP_LEN: \n",
" if(self->m > 1) {\n",
" return mp_obj_new_int(self->m);\n",
" } else {\n",
" return mp_obj_new_int(self->n);\n",
" }\n",
" break;\n",
" \n",
" case MP_UNARY_OP_INVERT:\n",
" if(self->array->typecode == NDARRAY_FLOAT) {\n",
" mp_raise_ValueError(translate(\"operation is not supported for given type\"));\n",
" }\n",
" // we can invert the content byte by byte, there is no need to distinguish \n",
" // between different typecodes\n",
" ndarray = MP_OBJ_TO_PTR(ndarray_copy(self_in));\n",
" uint8_t *array = (uint8_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->bytes; i++) array[i] = ~array[i];\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
" break;\n",
" \n",
" case MP_UNARY_OP_NEGATIVE:\n",
" ndarray = MP_OBJ_TO_PTR(ndarray_copy(self_in));\n",
" if(self->array->typecode == NDARRAY_UINT8) {\n",
" uint8_t *array = (uint8_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) array[i] = -array[i];\n",
" } else if(self->array->typecode == NDARRAY_INT8) {\n",
" int8_t *array = (int8_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) array[i] = -array[i];\n",
" } else if(self->array->typecode == NDARRAY_UINT16) { \n",
" uint16_t *array = (uint16_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) array[i] = -array[i];\n",
" } else if(self->array->typecode == NDARRAY_INT16) {\n",
" int16_t *array = (int16_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) array[i] = -array[i];\n",
" } else {\n",
" mp_float_t *array = (mp_float_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) array[i] = -array[i];\n",
" }\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
" break;\n",
"\n",
" case MP_UNARY_OP_POSITIVE:\n",
" return ndarray_copy(self_in);\n",
"\n",
" case MP_UNARY_OP_ABS:\n",
" if((self->array->typecode == NDARRAY_UINT8) || (self->array->typecode == NDARRAY_UINT16)) {\n",
" return ndarray_copy(self_in);\n",
" }\n",
" ndarray = MP_OBJ_TO_PTR(ndarray_copy(self_in));\n",
" if(self->array->typecode == NDARRAY_INT8) {\n",
" int8_t *array = (int8_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) {\n",
" if(array[i] < 0) array[i] = -array[i];\n",
" }\n",
" } else if(self->array->typecode == NDARRAY_INT16) {\n",
" int16_t *array = (int16_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) {\n",
" if(array[i] < 0) array[i] = -array[i];\n",
" }\n",
" } else {\n",
" mp_float_t *array = (mp_float_t *)ndarray->array->items;\n",
" for(size_t i=0; i < self->array->len; i++) {\n",
" if(array[i] < 0) array[i] = -array[i];\n",
" } \n",
" }\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
" break;\n",
" default: return MP_OBJ_NULL; // operator not supported\n",
" }\n",
"}\n",
"\n",
"mp_obj_t ndarray_transpose(mp_obj_t self_in) {\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" // the size of a single item in the array\n",
" uint8_t _sizeof = mp_binary_get_size('@', self->array->typecode, NULL);\n",
" \n",
" // NOTE: \n",
" // if the matrices are square, we can simply swap items, but \n",
" // generic matrices can't be transposed in place, so we have to \n",
" // declare a temporary variable\n",
" \n",
" // NOTE: \n",
" // In the old matrix, the coordinate (m, n) is m*self->n + n\n",
" // We have to assign this to the coordinate (n, m) in the new \n",
" // matrix, i.e., to n*self->m + m (since the new matrix has self->m columns)\n",
" \n",
" // one-dimensional arrays can be transposed by simply swapping the dimensions\n",
" if((self->m != 1) && (self->n != 1)) {\n",
" uint8_t *c = (uint8_t *)self->array->items;\n",
" // self->bytes is the size of the bytearray, irrespective of the typecode\n",
" uint8_t *tmp = m_new(uint8_t, self->bytes);\n",
" for(size_t m=0; m < self->m; m++) {\n",
" for(size_t n=0; n < self->n; n++) {\n",
" memcpy(tmp+_sizeof*(n*self->m + m), c+_sizeof*(m*self->n + n), _sizeof);\n",
" }\n",
" }\n",
" memcpy(self->array->items, tmp, self->bytes);\n",
" m_del(uint8_t, tmp, self->bytes);\n",
" } \n",
" SWAP(size_t, self->m, self->n);\n",
" return mp_const_none;\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_1(ndarray_transpose_obj, ndarray_transpose);\n",
"\n",
"mp_obj_t ndarray_reshape(mp_obj_t self_in, mp_obj_t shape) {\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" if(!MP_OBJ_IS_TYPE(shape, &mp_type_tuple) || (MP_OBJ_SMALL_INT_VALUE(mp_obj_len_maybe(shape)) != 2)) {\n",
" mp_raise_ValueError(translate(\"shape must be a 2-tuple\"));\n",
" }\n",
"\n",
" mp_obj_iter_buf_t iter_buf;\n",
" mp_obj_t item, iterable = mp_getiter(shape, &iter_buf);\n",
" uint16_t m, n;\n",
" item = mp_iternext(iterable);\n",
" m = mp_obj_get_int(item);\n",
" item = mp_iternext(iterable);\n",
" n = mp_obj_get_int(item);\n",
" if(m*n != self->m*self->n) {\n",
" // TODO: the proper error message would be \"cannot reshape array of size %d into shape (%d, %d)\"\n",
" mp_raise_ValueError(translate(\"cannot reshape array (incompatible input/output shape)\"));\n",
" }\n",
" self->m = m;\n",
" self->n = n;\n",
" return MP_OBJ_FROM_PTR(self);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_2(ndarray_reshape_obj, ndarray_reshape);\n",
"\n",
"mp_int_t ndarray_get_buffer(mp_obj_t self_in, mp_buffer_info_t *bufinfo, mp_uint_t flags) {\n",
" ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);\n",
" // buffer_p.get_buffer() returns zero for success, while mp_get_buffer returns true for success\n",
" return !mp_get_buffer(self->array, bufinfo, flags);\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Properties"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### ndarray_properties.h"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T06:10:22.005503Z",
"start_time": "2020-02-26T06:10:21.992913Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 1561 bytes to ndarray_properties.h\n"
]
}
],
"source": [
"%%ccode ndarray_properties.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2020 Jeff Epler for Adafruit Industries\n",
" * 2020 Zoltán Vörös \n",
"*/\n",
"\n",
"#ifndef _NDARRAY_PROPERTIES_\n",
"#define _NDARRAY_PROPERTIES_\n",
"\n",
"#include \"py/runtime.h\"\n",
"#include \"py/binary.h\"\n",
"#include \"py/obj.h\"\n",
"#include \"py/objarray.h\"\n",
"\n",
"#include \"ndarray.h\"\n",
"\n",
"#if CIRCUITPY\n",
"typedef struct _mp_obj_property_t {\n",
" mp_obj_base_t base;\n",
" mp_obj_t proxy[3]; // getter, setter, deleter\n",
"} mp_obj_property_t;\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_1(ndarray_get_shape_obj, ndarray_shape);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(ndarray_get_size_obj, ndarray_size);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(ndarray_get_itemsize_obj, ndarray_itemsize);\n",
"\n",
"STATIC const mp_obj_property_t ndarray_shape_obj = {\n",
" .base.type = &mp_type_property,\n",
" .proxy = {(mp_obj_t)&ndarray_get_shape_obj,\n",
" MP_ROM_NONE,\n",
" MP_ROM_NONE },\n",
"};\n",
"\n",
"STATIC const mp_obj_property_t ndarray_size_obj = {\n",
" .base.type = &mp_type_property,\n",
" .proxy = {(mp_obj_t)&ndarray_get_size_obj,\n",
" MP_ROM_NONE,\n",
" MP_ROM_NONE },\n",
"};\n",
"\n",
"STATIC const mp_obj_property_t ndarray_itemsize_obj = {\n",
" .base.type = &mp_type_property,\n",
" .proxy = {(mp_obj_t)&ndarray_get_itemsize_obj,\n",
" MP_ROM_NONE,\n",
" MP_ROM_NONE },\n",
"};\n",
"#else\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_1(ndarray_size_obj, ndarray_size);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(ndarray_itemsize_obj, ndarray_itemsize);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(ndarray_shape_obj, ndarray_shape);\n",
"\n",
"#endif\n",
"\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Linear algebra\n",
"\n",
"This module contains very basic matrix operators, such as transposing, reshaping, inverting, and matrix multiplication. The actual inversion is factored out into a helper function, so that the routine can be re-used in other modules. (The `polyfit` function in `poly.c` uses that.) Also note that inversion is based on the notion of a *small number* (epsilon). During the computation of the inverse, a number is treated as 0, if its absolute value is smaller than epsilon. This precaution is required, otherwise, one might run into singular matrices. \n",
"\n",
"As in `numpy`, `inv` is not a class method, but it should be applied only on `ndarray`s. This is why one has to check the argument type at the beginning of the functions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## linalg.h"
]
},
{
"cell_type": "code",
"execution_count": 257,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-14T07:00:02.772291Z",
"start_time": "2020-02-14T07:00:02.765336Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 634 bytes to linalg.h\n"
]
}
],
"source": [
"%%ccode linalg.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#ifndef _LINALG_\n",
"#define _LINALG_\n",
"\n",
"#include \"ulab.h\"\n",
"#include \"ndarray.h\"\n",
"\n",
"#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT\n",
"#define epsilon 1.2e-7\n",
"#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE\n",
"#define epsilon 2.3e-16\n",
"#endif\n",
"\n",
"#define JACOBI_MAX 20\n",
"\n",
"#if ULAB_LINALG_MODULE || ULAB_POLY_MODULE\n",
"bool linalg_invert_matrix(mp_float_t *, size_t );\n",
"#endif\n",
"\n",
"#if ULAB_LINALG_MODULE\n",
"\n",
"mp_obj_module_t ulab_linalg_module;\n",
"\n",
"#endif\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## linalg.c"
]
},
{
"cell_type": "code",
"execution_count": 282,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-14T14:00:26.902481Z",
"start_time": "2020-02-14T14:00:26.886650Z"
},
"code_folding": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 17193 bytes to linalg.c\n"
]
}
],
"source": [
"%%ccode linalg.c\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#include <stdlib.h>\n",
"#include <string.h>\n",
"#include <math.h>\n",
"#include \"py/obj.h\"\n",
"#include \"py/runtime.h\"\n",
"#include \"py/misc.h\"\n",
"#include \"linalg.h\"\n",
"\n",
"#if ULAB_LINALG_MODULE\n",
"\n",
"mp_obj_t linalg_size(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
"\n",
" if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"size is defined for ndarrays only\"));\n",
" } else {\n",
" ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);\n",
" if(args[1].u_obj == mp_const_none) {\n",
" return mp_obj_new_int(ndarray->array->len);\n",
" } else if(MP_OBJ_IS_INT(args[1].u_obj)) {\n",
" uint8_t ax = mp_obj_get_int(args[1].u_obj);\n",
" if(ax == 0) {\n",
" if(ndarray->m == 1) {\n",
" return mp_obj_new_int(ndarray->n);\n",
" } else {\n",
" return mp_obj_new_int(ndarray->m); \n",
" }\n",
" } else if(ax == 1) {\n",
" if(ndarray->m == 1) {\n",
" mp_raise_ValueError(translate(\"tuple index out of range\"));\n",
" } else {\n",
" return mp_obj_new_int(ndarray->n);\n",
" }\n",
" } else {\n",
" mp_raise_ValueError(translate(\"tuple index out of range\"));\n",
" }\n",
" } else {\n",
" mp_raise_TypeError(translate(\"wrong argument type\"));\n",
" }\n",
" }\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(linalg_size_obj, 1, linalg_size);\n",
"\n",
"bool linalg_invert_matrix(mp_float_t *data, size_t N) {\n",
" // returns true, of the inversion was successful, \n",
" // false, if the matrix is singular\n",
" \n",
" // initially, this is the unit matrix: the contents of this matrix is what \n",
" // will be returned after all the transformations\n",
" mp_float_t *unit = m_new(mp_float_t, N*N);\n",
"\n",
" mp_float_t elem = 1.0;\n",
" // initialise the unit matrix\n",
" memset(unit, 0, sizeof(mp_float_t)*N*N);\n",
" for(size_t m=0; m < N; m++) {\n",
" memcpy(&unit[m*(N+1)], &elem, sizeof(mp_float_t));\n",
" }\n",
" for(size_t m=0; m < N; m++){\n",
" // this could be faster with ((c < epsilon) && (c > -epsilon))\n",
" if(MICROPY_FLOAT_C_FUN(fabs)(data[m*(N+1)]) < epsilon) {\n",
" m_del(mp_float_t, unit, N*N);\n",
" return false;\n",
" }\n",
" for(size_t n=0; n < N; n++){\n",
" if(m != n){\n",
" elem = data[N*n+m] / data[m*(N+1)];\n",
" for(size_t k=0; k < N; k++){\n",
" data[N*n+k] -= elem * data[N*m+k];\n",
" unit[N*n+k] -= elem * unit[N*m+k];\n",
" }\n",
" }\n",
" }\n",
" }\n",
" for(size_t m=0; m < N; m++){ \n",
" elem = data[m*(N+1)];\n",
" for(size_t n=0; n < N; n++){\n",
" data[N*m+n] /= elem;\n",
" unit[N*m+n] /= elem;\n",
" }\n",
" }\n",
" memcpy(data, unit, sizeof(mp_float_t)*N*N);\n",
" m_del(mp_float_t, unit, N*N);\n",
" return true;\n",
"}\n",
"\n",
"mp_obj_t linalg_inv(mp_obj_t o_in) {\n",
" // since inv is not a class method, we have to inspect the input argument first\n",
" if(!MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"only ndarrays can be inverted\"));\n",
" }\n",
" ndarray_obj_t *o = MP_OBJ_TO_PTR(o_in);\n",
" if(!MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"only ndarray objects can be inverted\"));\n",
" }\n",
" if(o->m != o->n) {\n",
" mp_raise_ValueError(translate(\"only square matrices can be inverted\"));\n",
" }\n",
" ndarray_obj_t *inverted = create_new_ndarray(o->m, o->n, NDARRAY_FLOAT);\n",
" mp_float_t *data = (mp_float_t *)inverted->array->items;\n",
" mp_obj_t elem;\n",
" for(size_t m=0; m < o->m; m++) { // rows first\n",
" for(size_t n=0; n < o->n; n++) { // columns next\n",
" // this could, perhaps, be done in single line... \n",
" // On the other hand, we probably spend little time here\n",
" elem = mp_binary_get_val_array(o->array->typecode, o->array->items, m*o->n+n);\n",
" data[m*o->n+n] = (mp_float_t)mp_obj_get_float(elem);\n",
" }\n",
" }\n",
" \n",
" if(!linalg_invert_matrix(data, o->m)) {\n",
" // TODO: I am not sure this is needed here. Otherwise, \n",
" // how should we free up the unused RAM of inverted?\n",
" m_del(mp_float_t, inverted->array->items, o->n*o->n);\n",
" mp_raise_ValueError(translate(\"input matrix is singular\"));\n",
" }\n",
" return MP_OBJ_FROM_PTR(inverted);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_1(linalg_inv_obj, linalg_inv);\n",
"\n",
"mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) {\n",
" // TODO: should the results be upcast?\n",
" if(!MP_OBJ_IS_TYPE(_m1, &ulab_ndarray_type) || !MP_OBJ_IS_TYPE(_m2, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"arguments must be ndarrays\"));\n",
" }\n",
" ndarray_obj_t *m1 = MP_OBJ_TO_PTR(_m1);\n",
" ndarray_obj_t *m2 = MP_OBJ_TO_PTR(_m2); \n",
" if(m1->n != m2->m) {\n",
" mp_raise_ValueError(translate(\"matrix dimensions do not match\"));\n",
" }\n",
" // TODO: numpy uses upcasting here\n",
" ndarray_obj_t *out = create_new_ndarray(m1->m, m2->n, NDARRAY_FLOAT);\n",
" mp_float_t *outdata = (mp_float_t *)out->array->items;\n",
" mp_float_t sum, v1, v2;\n",
" for(size_t i=0; i < m1->m; i++) { // rows of m1\n",
" for(size_t j=0; j < m2->n; j++) { // columns of m2\n",
" sum = 0.0;\n",
" for(size_t k=0; k < m2->m; k++) {\n",
" // (i, k) * (k, j)\n",
" v1 = ndarray_get_float_value(m1->array->items, m1->array->typecode, i*m1->n+k);\n",
" v2 = ndarray_get_float_value(m2->array->items, m2->array->typecode, k*m2->n+j);\n",
" sum += v1 * v2;\n",
" }\n",
" outdata[j*m1->m+i] = sum;\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(out);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_2(linalg_dot_obj, linalg_dot);\n",
"\n",
"mp_obj_t linalg_zeros_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t kind) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_obj = MP_OBJ_NULL} } ,\n",
" { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" \n",
" uint8_t dtype = args[1].u_int;\n",
" if(!MP_OBJ_IS_INT(args[0].u_obj) && !MP_OBJ_IS_TYPE(args[0].u_obj, &mp_type_tuple)) {\n",
" mp_raise_TypeError(translate(\"input argument must be an integer or a 2-tuple\"));\n",
" }\n",
" ndarray_obj_t *ndarray = NULL;\n",
" if(MP_OBJ_IS_INT(args[0].u_obj)) {\n",
" size_t n = mp_obj_get_int(args[0].u_obj);\n",
" ndarray = create_new_ndarray(1, n, dtype);\n",
" } else if(MP_OBJ_IS_TYPE(args[0].u_obj, &mp_type_tuple)) {\n",
" mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(args[0].u_obj);\n",
" if(tuple->len != 2) {\n",
" mp_raise_TypeError(translate(\"input argument must be an integer or a 2-tuple\"));\n",
" }\n",
" ndarray = create_new_ndarray(mp_obj_get_int(tuple->items[0]), \n",
" mp_obj_get_int(tuple->items[1]), dtype);\n",
" }\n",
" if(kind == 1) {\n",
" mp_obj_t one = mp_obj_new_int(1);\n",
" for(size_t i=0; i < ndarray->array->len; i++) {\n",
" mp_binary_set_val_array(dtype, ndarray->array->items, i, one);\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
"}\n",
"\n",
"mp_obj_t linalg_zeros(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" return linalg_zeros_ones(n_args, pos_args, kw_args, 0);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(linalg_zeros_obj, 0, linalg_zeros);\n",
"\n",
"mp_obj_t linalg_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" return linalg_zeros_ones(n_args, pos_args, kw_args, 1);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(linalg_ones_obj, 0, linalg_ones);\n",
"\n",
"mp_obj_t linalg_eye(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_INT, {.u_int = 0} },\n",
" { MP_QSTR_M, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_k, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} }, \n",
" { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
"\n",
" size_t n = args[0].u_int, m;\n",
" int16_t k = args[2].u_int;\n",
" uint8_t dtype = args[3].u_int;\n",
" if(args[1].u_rom_obj == mp_const_none) {\n",
" m = n;\n",
" } else {\n",
" m = mp_obj_get_int(args[1].u_rom_obj);\n",
" }\n",
" \n",
" ndarray_obj_t *ndarray = create_new_ndarray(m, n, dtype);\n",
" mp_obj_t one = mp_obj_new_int(1);\n",
" size_t i = 0;\n",
" if((k >= 0) && (k < n)) {\n",
" while(k < n) {\n",
" mp_binary_set_val_array(dtype, ndarray->array->items, i*n+k, one);\n",
" k++;\n",
" i++;\n",
" }\n",
" } else if((k < 0) && (-k < m)) {\n",
" k = -k;\n",
" i = 0;\n",
" while(k < m) {\n",
" mp_binary_set_val_array(dtype, ndarray->array->items, k*n+i, one);\n",
" k++;\n",
" i++;\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(linalg_eye_obj, 0, linalg_eye);\n",
"\n",
"mp_obj_t linalg_det(mp_obj_t oin) {\n",
" if(!MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"function defined for ndarrays only\"));\n",
" }\n",
" ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);\n",
" if(in->m != in->n) {\n",
" mp_raise_ValueError(translate(\"input must be square matrix\"));\n",
" }\n",
" \n",
" mp_float_t *tmp = m_new(mp_float_t, in->n*in->n);\n",
" for(size_t i=0; i < in->array->len; i++){\n",
" tmp[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);\n",
" }\n",
" mp_float_t c;\n",
" for(size_t m=0; m < in->m-1; m++){\n",
" if(MICROPY_FLOAT_C_FUN(fabs)(tmp[m*(in->n+1)]) < epsilon) {\n",
" m_del(mp_float_t, tmp, in->n*in->n);\n",
" return mp_obj_new_float(0.0);\n",
" }\n",
" for(size_t n=0; n < in->n; n++){\n",
" if(m != n) {\n",
" c = tmp[in->n*n+m] / tmp[m*(in->n+1)];\n",
" for(size_t k=0; k < in->n; k++){\n",
" tmp[in->n*n+k] -= c * tmp[in->n*m+k];\n",
" }\n",
" }\n",
" }\n",
" }\n",
" mp_float_t det = 1.0;\n",
" \n",
" for(size_t m=0; m < in->m; m++){ \n",
" det *= tmp[m*(in->n+1)];\n",
" }\n",
" m_del(mp_float_t, tmp, in->n*in->n);\n",
" return mp_obj_new_float(det);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_1(linalg_det_obj, linalg_det);\n",
"\n",
"mp_obj_t linalg_eig(mp_obj_t oin) {\n",
" if(!MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"function defined for ndarrays only\"));\n",
" }\n",
" ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);\n",
" if(in->m != in->n) {\n",
" mp_raise_ValueError(translate(\"input must be square matrix\"));\n",
" }\n",
" mp_float_t *array = m_new(mp_float_t, in->array->len);\n",
" for(size_t i=0; i < in->array->len; i++) {\n",
" array[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);\n",
" }\n",
" // make sure the matrix is symmetric\n",
" for(size_t m=0; m < in->m; m++) {\n",
" for(size_t n=m+1; n < in->n; n++) {\n",
" // compare entry (m, n) to (n, m)\n",
" // TODO: this must probably be scaled!\n",
" if(epsilon < MICROPY_FLOAT_C_FUN(fabs)(array[m*in->n + n] - array[n*in->n + m])) {\n",
" mp_raise_ValueError(translate(\"input matrix is asymmetric\"));\n",
" }\n",
" }\n",
" }\n",
" \n",
" // if we got this far, then the matrix will be symmetric\n",
" \n",
" ndarray_obj_t *eigenvectors = create_new_ndarray(in->m, in->n, NDARRAY_FLOAT);\n",
" mp_float_t *eigvectors = (mp_float_t *)eigenvectors->array->items;\n",
" // start out with the unit matrix\n",
" for(size_t m=0; m < in->m; m++) {\n",
" eigvectors[m*(in->n+1)] = 1.0;\n",
" }\n",
" mp_float_t largest, w, t, c, s, tau, aMk, aNk, vm, vn;\n",
" size_t M, N;\n",
" size_t iterations = JACOBI_MAX*in->n*in->n;\n",
" do {\n",
" iterations--;\n",
" // find the pivot here\n",
" M = 0;\n",
" N = 0;\n",
" largest = 0.0;\n",
" for(size_t m=0; m < in->m-1; m++) { // -1: no need to inspect last row\n",
" for(size_t n=m+1; n < in->n; n++) {\n",
" w = MICROPY_FLOAT_C_FUN(fabs)(array[m*in->n + n]);\n",
" if((largest < w) && (epsilon < w)) {\n",
" M = m;\n",
" N = n;\n",
" largest = w;\n",
" }\n",
" }\n",
" }\n",
" if(M+N == 0) { // all entries are smaller than epsilon, there is not much we can do...\n",
" break;\n",
" }\n",
" // at this point, we have the pivot, and it is the entry (M, N)\n",
" // now we have to find the rotation angle\n",
" w = (array[N*in->n + N] - array[M*in->n + M]) / (2.0*array[M*in->n + N]);\n",
" // The following if/else chooses the smaller absolute value for the tangent \n",
" // of the rotation angle. Going with the smaller should be numerically stabler.\n",
" if(w > 0) {\n",
" t = MICROPY_FLOAT_C_FUN(sqrt)(w*w + 1.0) - w;\n",
" } else {\n",
" t = -1.0*(MICROPY_FLOAT_C_FUN(sqrt)(w*w + 1.0) + w);\n",
" }\n",
" s = t / MICROPY_FLOAT_C_FUN(sqrt)(t*t + 1.0); // the sine of the rotation angle\n",
" c = 1.0 / MICROPY_FLOAT_C_FUN(sqrt)(t*t + 1.0); // the cosine of the rotation angle\n",
" tau = (1.0-c)/s; // this is equal to the tangent of the half of the rotation angle\n",
" \n",
" // at this point, we have the rotation angles, so we can transform the matrix\n",
" // first the two diagonal elements\n",
" // a(M, M) = a(M, M) - t*a(M, N)\n",
" array[M*in->n + M] = array[M*in->n + M] - t * array[M*in->n + N];\n",
" // a(N, N) = a(N, N) + t*a(M, N)\n",
" array[N*in->n + N] = array[N*in->n + N] + t * array[M*in->n + N];\n",
" // after the rotation, the a(M, N), and a(N, M) entries should become zero\n",
" array[M*in->n + N] = array[N*in->n + M] = 0.0;\n",
" // then all other elements in the column\n",
" for(size_t k=0; k < in->m; k++) {\n",
" if((k == M) || (k == N)) {\n",
" continue;\n",
" }\n",
" aMk = array[M*in->n + k];\n",
" aNk = array[N*in->n + k];\n",
" // a(M, k) = a(M, k) - s*(a(N, k) + tau*a(M, k))\n",
" array[M*in->n + k] -= s*(aNk + tau*aMk);\n",
" // a(N, k) = a(N, k) + s*(a(M, k) - tau*a(N, k))\n",
" array[N*in->n + k] += s*(aMk - tau*aNk);\n",
" // a(k, M) = a(M, k)\n",
" array[k*in->n + M] = array[M*in->n + k];\n",
" // a(k, N) = a(N, k)\n",
" array[k*in->n + N] = array[N*in->n + k];\n",
" }\n",
" // now we have to update the eigenvectors\n",
" // the rotation matrix, R, multiplies from the right\n",
" // R is the unit matrix, except for the \n",
" // R(M,M) = R(N, N) = c\n",
" // R(N, M) = s\n",
" // (M, N) = -s\n",
" // entries. This means that only the Mth, and Nth columns will change\n",
" for(size_t m=0; m < in->m; m++) {\n",
" vm = eigvectors[m*in->n+M];\n",
" vn = eigvectors[m*in->n+N];\n",
" // the new value of eigvectors(m, M)\n",
" eigvectors[m*in->n+M] = c * vm - s * vn;\n",
" // the new value of eigvectors(m, N)\n",
" eigvectors[m*in->n+N] = s * vm + c * vn;\n",
" }\n",
" } while(iterations > 0);\n",
" \n",
" if(iterations == 0) { \n",
" // the computation did not converge; numpy raises LinAlgError\n",
" m_del(mp_float_t, array, in->array->len);\n",
" mp_raise_ValueError(translate(\"iterations did not converge\"));\n",
" }\n",
" ndarray_obj_t *eigenvalues = create_new_ndarray(1, in->n, NDARRAY_FLOAT);\n",
" mp_float_t *eigvalues = (mp_float_t *)eigenvalues->array->items;\n",
" for(size_t i=0; i < in->n; i++) {\n",
" eigvalues[i] = array[i*(in->n+1)];\n",
" }\n",
" m_del(mp_float_t, array, in->array->len);\n",
" \n",
" mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(mp_obj_new_tuple(2, NULL));\n",
" tuple->items[0] = MP_OBJ_FROM_PTR(eigenvalues);\n",
" tuple->items[1] = MP_OBJ_FROM_PTR(eigenvectors);\n",
" return tuple;\n",
" return MP_OBJ_FROM_PTR(eigenvalues);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_1(linalg_eig_obj, linalg_eig);\n",
"\n",
"STATIC const mp_rom_map_elem_t ulab_linalg_globals_table[] = {\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_linalg) },\n",
" { MP_ROM_QSTR(MP_QSTR_size), (mp_obj_t)&linalg_size_obj },\n",
" { MP_ROM_QSTR(MP_QSTR_inv), (mp_obj_t)&linalg_inv_obj },\n",
" { MP_ROM_QSTR(MP_QSTR_dot), (mp_obj_t)&linalg_dot_obj },\n",
" { MP_ROM_QSTR(MP_QSTR_zeros), (mp_obj_t)&linalg_zeros_obj },\n",
" { MP_ROM_QSTR(MP_QSTR_ones), (mp_obj_t)&linalg_ones_obj },\n",
" { MP_ROM_QSTR(MP_QSTR_eye), (mp_obj_t)&linalg_eye_obj },\n",
" { MP_ROM_QSTR(MP_QSTR_det), (mp_obj_t)&linalg_det_obj },\n",
" { MP_ROM_QSTR(MP_QSTR_eig), (mp_obj_t)&linalg_eig_obj }, \n",
"};\n",
"\n",
"STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_linalg_globals, ulab_linalg_globals_table);\n",
"\n",
"mp_obj_module_t ulab_linalg_module = {\n",
" .base = { &mp_type_module },\n",
" .globals = (mp_obj_dict_t*)&mp_module_ulab_linalg_globals,\n",
"};\n",
"\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Vectorising mathematical operations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## vectorise.h"
]
},
{
"cell_type": "code",
"execution_count": 633,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-16T18:26:39.726068Z",
"start_time": "2020-02-16T18:26:39.716040Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 736 bytes to vectorise.h\n"
]
}
],
"source": [
"%%ccode vectorise.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#ifndef _VECTORISE_\n",
"#define _VECTORISE_\n",
"\n",
"#include \"ulab.h\"\n",
"#include \"ndarray.h\"\n",
"\n",
"#if ULAB_VECTORISE_MODULE\n",
"\n",
"mp_obj_module_t ulab_vectorise_module;\n",
"\n",
"#define ITERATE_VECTOR(type, source, out) do {\\\n",
" type *input = (type *)(source)->array->items;\\\n",
" for(size_t i=0; i < (source)->array->len; i++) {\\\n",
" (out)[i] = f(input[i]);\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"#define MATH_FUN_1(py_name, c_name) \\\n",
" mp_obj_t vectorise_ ## py_name(mp_obj_t x_obj) { \\\n",
" return vectorise_generic_vector(x_obj, MICROPY_FLOAT_C_FUN(c_name)); \\\n",
"}\n",
" \n",
"#endif\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## vectorise.c"
]
},
{
"cell_type": "code",
"execution_count": 321,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-14T18:50:23.939743Z",
"start_time": "2020-02-14T18:50:23.934622Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 6477 bytes to vectorise.c\n"
]
}
],
"source": [
"%%ccode vectorise.c\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#include <math.h>\n",
"#include <stdio.h>\n",
"#include <stdlib.h>\n",
"#include \"py/runtime.h\"\n",
"#include \"py/binary.h\"\n",
"#include \"py/obj.h\"\n",
"#include \"py/objarray.h\"\n",
"#include \"vectorise.h\"\n",
"\n",
"#ifndef MP_PI\n",
"#define MP_PI MICROPY_FLOAT_CONST(3.14159265358979323846)\n",
"#endif\n",
" \n",
"#if ULAB_VECTORISE_MODULE\n",
"mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) {\n",
" // Return a single value, if o_in is not iterable\n",
" if(mp_obj_is_float(o_in) || MP_OBJ_IS_INT(o_in)) {\n",
" return mp_obj_new_float(f(mp_obj_get_float(o_in)));\n",
" }\n",
" mp_float_t x;\n",
" if(MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {\n",
" ndarray_obj_t *source = MP_OBJ_TO_PTR(o_in);\n",
" ndarray_obj_t *ndarray = create_new_ndarray(source->m, source->n, NDARRAY_FLOAT);\n",
" mp_float_t *dataout = (mp_float_t *)ndarray->array->items;\n",
" if(source->array->typecode == NDARRAY_UINT8) {\n",
" ITERATE_VECTOR(uint8_t, source, dataout);\n",
" } else if(source->array->typecode == NDARRAY_INT8) {\n",
" ITERATE_VECTOR(int8_t, source, dataout);\n",
" } else if(source->array->typecode == NDARRAY_UINT16) {\n",
" ITERATE_VECTOR(uint16_t, source, dataout);\n",
" } else if(source->array->typecode == NDARRAY_INT16) {\n",
" ITERATE_VECTOR(int16_t, source, dataout);\n",
" } else {\n",
" ITERATE_VECTOR(mp_float_t, source, dataout);\n",
" }\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
" } else if(MP_OBJ_IS_TYPE(o_in, &mp_type_tuple) || MP_OBJ_IS_TYPE(o_in, &mp_type_list) || \n",
" MP_OBJ_IS_TYPE(o_in, &mp_type_range)) { // i.e., the input is a generic iterable\n",
" mp_obj_array_t *o = MP_OBJ_TO_PTR(o_in);\n",
" ndarray_obj_t *out = create_new_ndarray(1, o->len, NDARRAY_FLOAT);\n",
" mp_float_t *dataout = (mp_float_t *)out->array->items;\n",
" mp_obj_iter_buf_t iter_buf;\n",
" mp_obj_t item, iterable = mp_getiter(o_in, &iter_buf);\n",
" size_t i=0;\n",
" while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" x = mp_obj_get_float(item);\n",
" dataout[i++] = f(x);\n",
" }\n",
" return MP_OBJ_FROM_PTR(out);\n",
" }\n",
" return mp_const_none;\n",
"}\n",
"\n",
"\n",
"MATH_FUN_1(acos, acos);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acos_obj, vectorise_acos);\n",
"\n",
"MATH_FUN_1(acosh, acosh);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acosh_obj, vectorise_acosh);\n",
"\n",
"MATH_FUN_1(asin, asin);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_asin_obj, vectorise_asin);\n",
"\n",
"MATH_FUN_1(asinh, asinh);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_asinh_obj, vectorise_asinh);\n",
"\n",
"MATH_FUN_1(atan, atan);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_atan_obj, vectorise_atan);\n",
"\n",
"MATH_FUN_1(atanh, atanh);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_atanh_obj, vectorise_atanh);\n",
"\n",
"MATH_FUN_1(ceil, ceil);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_ceil_obj, vectorise_ceil);\n",
"\n",
"MATH_FUN_1(cos, cos);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_cos_obj, vectorise_cos);\n",
"\n",
"MATH_FUN_1(cosh, cosh);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_cosh_obj, vectorise_cosh);\n",
"\n",
"MATH_FUN_1(erf, erf);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_erf_obj, vectorise_erf);\n",
"\n",
"MATH_FUN_1(erfc, erfc);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_erfc_obj, vectorise_erfc);\n",
"\n",
"MATH_FUN_1(exp, exp);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_exp_obj, vectorise_exp);\n",
"\n",
"MATH_FUN_1(expm1, expm1);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_expm1_obj, vectorise_expm1);\n",
"\n",
"MATH_FUN_1(floor, floor);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_floor_obj, vectorise_floor);\n",
"\n",
"MATH_FUN_1(gamma, tgamma);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_gamma_obj, vectorise_gamma);\n",
"\n",
"MATH_FUN_1(lgamma, lgamma);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_lgamma_obj, vectorise_lgamma);\n",
"\n",
"MATH_FUN_1(log, log);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log_obj, vectorise_log);\n",
"\n",
"MATH_FUN_1(log10, log10);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log10_obj, vectorise_log10);\n",
"\n",
"MATH_FUN_1(log2, log2);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log2_obj, vectorise_log2);\n",
"\n",
"MATH_FUN_1(sin, sin);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sin_obj, vectorise_sin);\n",
"\n",
"MATH_FUN_1(sinh, sinh);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sinh_obj, vectorise_sinh);\n",
"\n",
"MATH_FUN_1(sqrt, sqrt);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sqrt_obj, vectorise_sqrt);\n",
"\n",
"MATH_FUN_1(tan, tan);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tan_obj, vectorise_tan);\n",
"\n",
"MATH_FUN_1(tanh, tanh);\n",
"MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tanh_obj, vectorise_tanh);\n",
"\n",
"STATIC const mp_rom_map_elem_t ulab_vectorise_globals_table[] = {\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_vector) },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_acos), (mp_obj_t)&vectorise_acos_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_acosh), (mp_obj_t)&vectorise_acosh_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_asin), (mp_obj_t)&vectorise_asin_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_asinh), (mp_obj_t)&vectorise_asinh_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_atan), (mp_obj_t)&vectorise_atan_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_atanh), (mp_obj_t)&vectorise_atanh_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_ceil), (mp_obj_t)&vectorise_ceil_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_cos), (mp_obj_t)&vectorise_cos_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_erf), (mp_obj_t)&vectorise_erf_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_erfc), (mp_obj_t)&vectorise_erfc_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_exp), (mp_obj_t)&vectorise_exp_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_expm1), (mp_obj_t)&vectorise_expm1_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_floor), (mp_obj_t)&vectorise_floor_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_gamma), (mp_obj_t)&vectorise_gamma_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_lgamma), (mp_obj_t)&vectorise_lgamma_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_log), (mp_obj_t)&vectorise_log_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_log10), (mp_obj_t)&vectorise_log10_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_log2), (mp_obj_t)&vectorise_log2_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_sin), (mp_obj_t)&vectorise_sin_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_sinh), (mp_obj_t)&vectorise_sinh_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_sqrt), (mp_obj_t)&vectorise_sqrt_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_tan), (mp_obj_t)&vectorise_tan_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_tanh), (mp_obj_t)&vectorise_tanh_obj },\n",
"};\n",
"\n",
"STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_vectorise_globals, ulab_vectorise_globals_table);\n",
"\n",
"mp_obj_module_t ulab_vectorise_module = {\n",
" .base = { &mp_type_module },\n",
" .globals = (mp_obj_dict_t*)&mp_module_ulab_vectorise_globals,\n",
"};\n",
"\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Polynomials"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## poly.h"
]
},
{
"cell_type": "code",
"execution_count": 153,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-13T19:54:43.020113Z",
"start_time": "2020-02-13T19:54:43.015257Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 303 bytes to poly.h\n"
]
}
],
"source": [
"%%ccode poly.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#ifndef _POLY_\n",
"#define _POLY_\n",
"\n",
"#include \"ulab.h\"\n",
"\n",
"#if ULAB_POLY_MODULE\n",
"\n",
"mp_obj_module_t ulab_poly_module;\n",
"\n",
"#endif\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## poly.c"
]
},
{
"cell_type": "code",
"execution_count": 406,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-15T10:16:25.378564Z",
"start_time": "2020-02-15T10:16:25.351584Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 7603 bytes to poly.c\n"
]
}
],
"source": [
"%%ccode poly.c\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#include \"py/obj.h\"\n",
"#include \"py/runtime.h\"\n",
"#include \"py/objarray.h\"\n",
"#include \"ndarray.h\"\n",
"#include \"linalg.h\"\n",
"#include \"poly.h\"\n",
"\n",
"#if ULAB_POLY_MODULE\n",
"bool object_is_nditerable(mp_obj_t o_in) {\n",
" if(MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type) || \n",
" MP_OBJ_IS_TYPE(o_in, &mp_type_tuple) || \n",
" MP_OBJ_IS_TYPE(o_in, &mp_type_list) || \n",
" MP_OBJ_IS_TYPE(o_in, &mp_type_range)) {\n",
" return true;\n",
" }\n",
" return false;\n",
"}\n",
"\n",
"size_t get_nditerable_len(mp_obj_t o_in) {\n",
" if(MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {\n",
" ndarray_obj_t *in = MP_OBJ_TO_PTR(o_in);\n",
" return in->array->len;\n",
" } else {\n",
" return (size_t)mp_obj_get_int(mp_obj_len_maybe(o_in));\n",
" }\n",
"}\n",
"\n",
"mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {\n",
" // TODO: return immediately, if o_p is not an iterable\n",
" // TODO: there is a bug here: matrices won't work, \n",
" // because there is a single iteration loop\n",
" size_t m, n;\n",
" if(MP_OBJ_IS_TYPE(o_x, &ulab_ndarray_type)) {\n",
" ndarray_obj_t *ndx = MP_OBJ_TO_PTR(o_x);\n",
" m = ndx->m;\n",
" n = ndx->n;\n",
" } else {\n",
" mp_obj_array_t *ix = MP_OBJ_TO_PTR(o_x);\n",
" m = 1;\n",
" n = ix->len;\n",
" }\n",
" // polynomials are going to be of type float, except, when both \n",
" // the coefficients and the independent variable are integers\n",
" ndarray_obj_t *out = create_new_ndarray(m, n, NDARRAY_FLOAT);\n",
" mp_obj_iter_buf_t x_buf;\n",
" mp_obj_t x_item, x_iterable = mp_getiter(o_x, &x_buf);\n",
"\n",
" mp_obj_iter_buf_t p_buf;\n",
" mp_obj_t p_item, p_iterable;\n",
"\n",
" mp_float_t x, y;\n",
" mp_float_t *outf = (mp_float_t *)out->array->items;\n",
" uint8_t plen = mp_obj_get_int(mp_obj_len_maybe(o_p));\n",
" mp_float_t *p = m_new(mp_float_t, plen);\n",
" p_iterable = mp_getiter(o_p, &p_buf);\n",
" uint16_t i = 0; \n",
" while((p_item = mp_iternext(p_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" p[i] = mp_obj_get_float(p_item);\n",
" i++;\n",
" }\n",
" i = 0;\n",
" while ((x_item = mp_iternext(x_iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" x = mp_obj_get_float(x_item);\n",
" y = p[0];\n",
" for(uint8_t j=0; j < plen-1; j++) {\n",
" y *= x;\n",
" y += p[j+1];\n",
" }\n",
" outf[i++] = y;\n",
" }\n",
" m_del(mp_float_t, p, plen);\n",
" return MP_OBJ_FROM_PTR(out);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_2(poly_polyval_obj, poly_polyval);\n",
"\n",
"mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) {\n",
" if((n_args != 2) && (n_args != 3)) {\n",
" mp_raise_ValueError(translate(\"number of arguments must be 2, or 3\"));\n",
" }\n",
" if(!object_is_nditerable(args[0])) {\n",
" mp_raise_ValueError(translate(\"input data must be an iterable\"));\n",
" }\n",
" uint16_t lenx = 0, leny = 0;\n",
" uint8_t deg = 0;\n",
" mp_float_t *x, *XT, *y, *prod;\n",
"\n",
" if(n_args == 2) { // only the y values are supplied\n",
" // TODO: this is actually not enough: the first argument can very well be a matrix, \n",
" // in which case we are between the rock and a hard place\n",
" leny = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));\n",
" deg = (uint8_t)mp_obj_get_int(args[1]);\n",
" if(leny < deg) {\n",
" mp_raise_ValueError(translate(\"more degrees of freedom than data points\"));\n",
" }\n",
" lenx = leny;\n",
" x = m_new(mp_float_t, lenx); // assume uniformly spaced data points\n",
" for(size_t i=0; i < lenx; i++) {\n",
" x[i] = i;\n",
" }\n",
" y = m_new(mp_float_t, leny);\n",
" fill_array_iterable(y, args[0]);\n",
" } else if(n_args == 3) {\n",
" lenx = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));\n",
" leny = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));\n",
" if(lenx != leny) {\n",
" mp_raise_ValueError(translate(\"input vectors must be of equal length\"));\n",
" }\n",
" deg = (uint8_t)mp_obj_get_int(args[2]);\n",
" if(leny < deg) {\n",
" mp_raise_ValueError(translate(\"more degrees of freedom than data points\"));\n",
" }\n",
" x = m_new(mp_float_t, lenx);\n",
" fill_array_iterable(x, args[0]);\n",
" y = m_new(mp_float_t, leny);\n",
" fill_array_iterable(y, args[1]);\n",
" }\n",
" \n",
" // one could probably express X as a function of XT, \n",
" // and thereby save RAM, because X is used only in the product\n",
" XT = m_new(mp_float_t, (deg+1)*leny); // XT is a matrix of shape (deg+1, len) (rows, columns)\n",
" for(uint8_t i=0; i < leny; i++) { // column index\n",
" XT[i+0*lenx] = 1.0; // top row\n",
" for(uint8_t j=1; j < deg+1; j++) { // row index\n",
" XT[i+j*leny] = XT[i+(j-1)*leny]*x[i];\n",
" }\n",
" }\n",
" \n",
" prod = m_new(mp_float_t, (deg+1)*(deg+1)); // the product matrix is of shape (deg+1, deg+1)\n",
" mp_float_t sum;\n",
" for(uint16_t i=0; i < deg+1; i++) { // column index\n",
" for(uint16_t j=0; j < deg+1; j++) { // row index\n",
" sum = 0.0;\n",
" for(size_t k=0; k < lenx; k++) {\n",
" // (j, k) * (k, i) \n",
" // Note that the second matrix is simply the transpose of the first: \n",
" // X(k, i) = XT(i, k) = XT[k*lenx+i]\n",
" sum += XT[j*lenx+k]*XT[i*lenx+k]; // X[k*(deg+1)+i];\n",
" }\n",
" prod[j*(deg+1)+i] = sum;\n",
" }\n",
" }\n",
" if(!linalg_invert_matrix(prod, deg+1)) {\n",
" // Although X was a Vandermonde matrix, whose inverse is guaranteed to exist, \n",
" // we bail out here, if prod couldn't be inverted: if the values in x are not all \n",
" // distinct, prod is singular\n",
" m_del(mp_float_t, XT, (deg+1)*lenx);\n",
" m_del(mp_float_t, x, lenx);\n",
" m_del(mp_float_t, y, lenx);\n",
" m_del(mp_float_t, prod, (deg+1)*(deg+1));\n",
" mp_raise_ValueError(translate(\"could not invert Vandermonde matrix\"));\n",
" } \n",
" // at this point, we have the inverse of X^T * X\n",
" // y is a column vector; x is free now, we can use it for storing intermediate values\n",
" for(uint16_t i=0; i < deg+1; i++) { // row index\n",
" sum = 0.0;\n",
" for(uint16_t j=0; j < lenx; j++) { // column index\n",
" sum += XT[i*lenx+j]*y[j];\n",
" }\n",
" x[i] = sum;\n",
" }\n",
" // XT is no longer needed\n",
" m_del(mp_float_t, XT, (deg+1)*leny);\n",
" \n",
" ndarray_obj_t *beta = create_new_ndarray(deg+1, 1, NDARRAY_FLOAT);\n",
" mp_float_t *betav = (mp_float_t *)beta->array->items;\n",
" // x[0..(deg+1)] contains now the product X^T * y; we can get rid of y\n",
" m_del(float, y, leny);\n",
" \n",
" // now, we calculate beta, i.e., we apply prod = (X^T * X)^(-1) on x = X^T * y; x is a column vector now\n",
" for(uint16_t i=0; i < deg+1; i++) {\n",
" sum = 0.0;\n",
" for(uint16_t j=0; j < deg+1; j++) {\n",
" sum += prod[i*(deg+1)+j]*x[j];\n",
" }\n",
" betav[i] = sum;\n",
" }\n",
" m_del(mp_float_t, x, lenx);\n",
" m_del(mp_float_t, prod, (deg+1)*(deg+1));\n",
" for(uint8_t i=0; i < (deg+1)/2; i++) {\n",
" // We have to reverse the array, for the leading coefficient comes first. \n",
" SWAP(mp_float_t, betav[i], betav[deg-i]);\n",
" }\n",
" return MP_OBJ_FROM_PTR(beta);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(poly_polyfit_obj, 2, 3, poly_polyfit);\n",
"\n",
"STATIC const mp_rom_map_elem_t ulab_poly_globals_table[] = {\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_poly) },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_polyval), (mp_obj_t)&poly_polyval_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_polyfit), (mp_obj_t)&poly_polyfit_obj },\n",
"};\n",
"\n",
"STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_poly_globals, ulab_poly_globals_table);\n",
"\n",
"mp_obj_module_t ulab_poly_module = {\n",
" .base = { &mp_type_module },\n",
" .globals = (mp_obj_dict_t*)&mp_module_ulab_poly_globals,\n",
"};\n",
"\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fast Fourier transform"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## fft.h"
]
},
{
"cell_type": "code",
"execution_count": 149,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-13T19:50:01.802329Z",
"start_time": "2020-02-13T19:50:01.791668Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 432 bytes to fft.h\n"
]
}
],
"source": [
"%%ccode fft.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#ifndef _FFT_\n",
"#define _FFT_\n",
"#include \"ulab.h\"\n",
"\n",
"#ifndef MP_PI\n",
"#define MP_PI MICROPY_FLOAT_CONST(3.14159265358979323846)\n",
"#endif\n",
"\n",
"#define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }\n",
"\n",
"#if ULAB_FFT_MODULE\n",
"\n",
"mp_obj_module_t ulab_fft_module;\n",
"\n",
"#endif\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## fft.c"
]
},
{
"cell_type": "code",
"execution_count": 628,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-16T18:23:54.268743Z",
"start_time": "2020-02-16T18:23:54.248822Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 6360 bytes to fft.c\n"
]
}
],
"source": [
"%%ccode fft.c\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#include <math.h>\n",
"#include <stdio.h>\n",
"#include <stdlib.h>\n",
"#include <string.h>\n",
"#include \"py/runtime.h\"\n",
"#include \"py/builtin.h\"\n",
"#include \"py/binary.h\"\n",
"#include \"py/obj.h\"\n",
"#include \"py/objarray.h\"\n",
"#include \"ndarray.h\"\n",
"#include \"fft.h\"\n",
"\n",
"#if ULAB_FFT_MODULE\n",
"\n",
"enum FFT_TYPE {\n",
" FFT_FFT,\n",
" FFT_IFFT,\n",
" FFT_SPECTRUM,\n",
"};\n",
"\n",
"void fft_kernel(mp_float_t *real, mp_float_t *imag, int n, int isign) {\n",
" // This is basically a modification of four1 from Numerical Recipes\n",
" // The main difference is that this function takes two arrays, one \n",
" // for the real, and one for the imaginary parts. \n",
" int j, m, mmax, istep;\n",
" mp_float_t tempr, tempi;\n",
" mp_float_t wtemp, wr, wpr, wpi, wi, theta;\n",
"\n",
" j = 0;\n",
" for(int i = 0; i < n; i++) {\n",
" if (j > i) {\n",
" SWAP(mp_float_t, real[i], real[j]);\n",
" SWAP(mp_float_t, imag[i], imag[j]);\n",
" }\n",
" m = n >> 1;\n",
" while (j >= m && m > 0) {\n",
" j -= m;\n",
" m >>= 1;\n",
" }\n",
" j += m;\n",
" }\n",
"\n",
" mmax = 1;\n",
" while (n > mmax) {\n",
" istep = mmax << 1;\n",
" theta = -2.0*isign*MP_PI/istep;\n",
" wtemp = MICROPY_FLOAT_C_FUN(sin)(0.5 * theta);\n",
" wpr = -2.0 * wtemp * wtemp;\n",
" wpi = MICROPY_FLOAT_C_FUN(sin)(theta);\n",
" wr = 1.0;\n",
" wi = 0.0;\n",
" for(m = 0; m < mmax; m++) {\n",
" for(int i = m; i < n; i += istep) {\n",
" j = i + mmax;\n",
" tempr = wr * real[j] - wi * imag[j];\n",
" tempi = wr * imag[j] + wi * real[j];\n",
" real[j] = real[i] - tempr;\n",
" imag[j] = imag[i] - tempi;\n",
" real[i] += tempr;\n",
" imag[i] += tempi;\n",
" }\n",
" wtemp = wr;\n",
" wr = wr*wpr - wi*wpi + wr;\n",
" wi = wi*wpr + wtemp*wpi + wi;\n",
" }\n",
" mmax = istep;\n",
" }\n",
"}\n",
"\n",
"mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im, uint8_t type) {\n",
" if(!MP_OBJ_IS_TYPE(arg_re, &ulab_ndarray_type)) {\n",
" mp_raise_NotImplementedError(translate(\"FFT is defined for ndarrays only\"));\n",
" } \n",
" if(n_args == 2) {\n",
" if(!MP_OBJ_IS_TYPE(arg_im, &ulab_ndarray_type)) {\n",
" mp_raise_NotImplementedError(translate(\"FFT is defined for ndarrays only\"));\n",
" }\n",
" }\n",
" // Check if input is of length of power of 2\n",
" ndarray_obj_t *re = MP_OBJ_TO_PTR(arg_re);\n",
" uint16_t len = re->array->len;\n",
" if((len & (len-1)) != 0) {\n",
" mp_raise_ValueError(translate(\"input array length must be power of 2\"));\n",
" }\n",
" \n",
" ndarray_obj_t *out_re = create_new_ndarray(1, len, NDARRAY_FLOAT);\n",
" mp_float_t *data_re = (mp_float_t *)out_re->array->items;\n",
" \n",
" if(re->array->typecode == NDARRAY_FLOAT) { \n",
" // By treating this case separately, we can save a bit of time.\n",
" // I don't know if it is worthwhile, though...\n",
" memcpy((mp_float_t *)out_re->array->items, (mp_float_t *)re->array->items, re->bytes);\n",
" } else {\n",
" for(size_t i=0; i < len; i++) {\n",
" *data_re++ = ndarray_get_float_value(re->array->items, re->array->typecode, i);\n",
" }\n",
" data_re -= len;\n",
" }\n",
" ndarray_obj_t *out_im = create_new_ndarray(1, len, NDARRAY_FLOAT);\n",
" mp_float_t *data_im = (mp_float_t *)out_im->array->items;\n",
"\n",
" if(n_args == 2) {\n",
" ndarray_obj_t *im = MP_OBJ_TO_PTR(arg_im);\n",
" if (re->array->len != im->array->len) {\n",
" mp_raise_ValueError(translate(\"real and imaginary parts must be of equal length\"));\n",
" }\n",
" if(im->array->typecode == NDARRAY_FLOAT) {\n",
" memcpy((mp_float_t *)out_im->array->items, (mp_float_t *)im->array->items, im->bytes);\n",
" } else {\n",
" for(size_t i=0; i < len; i++) {\n",
" *data_im++ = ndarray_get_float_value(im->array->items, im->array->typecode, i);\n",
" }\n",
" data_im -= len;\n",
" }\n",
" }\n",
"\n",
" if((type == FFT_FFT) || (type == FFT_SPECTRUM)) {\n",
" fft_kernel(data_re, data_im, len, 1);\n",
" if(type == FFT_SPECTRUM) {\n",
" for(size_t i=0; i < len; i++) {\n",
" *data_re = MICROPY_FLOAT_C_FUN(sqrt)(*data_re * *data_re + *data_im * *data_im);\n",
" data_re++;\n",
" data_im++;\n",
" }\n",
" }\n",
" } else { // inverse transform\n",
" fft_kernel(data_re, data_im, len, -1);\n",
" // TODO: numpy accepts the norm keyword argument\n",
" for(size_t i=0; i < len; i++) {\n",
" *data_re++ /= len;\n",
" *data_im++ /= len;\n",
" }\n",
" }\n",
" if(type == FFT_SPECTRUM) {\n",
" return MP_OBJ_TO_PTR(out_re);\n",
" } else {\n",
" mp_obj_t tuple[2];\n",
" tuple[0] = out_re;\n",
" tuple[1] = out_im;\n",
" return mp_obj_new_tuple(2, tuple);\n",
" }\n",
"}\n",
"\n",
"mp_obj_t fft_fft(size_t n_args, const mp_obj_t *args) {\n",
" if(n_args == 2) {\n",
" return fft_fft_ifft_spectrum(n_args, args[0], args[1], FFT_FFT);\n",
" } else {\n",
" return fft_fft_ifft_spectrum(n_args, args[0], mp_const_none, FFT_FFT); \n",
" }\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj, 1, 2, fft_fft);\n",
"\n",
"mp_obj_t fft_ifft(size_t n_args, const mp_obj_t *args) {\n",
" if(n_args == 2) {\n",
" return fft_fft_ifft_spectrum(n_args, args[0], args[1], FFT_IFFT);\n",
" } else {\n",
" return fft_fft_ifft_spectrum(n_args, args[0], mp_const_none, FFT_IFFT);\n",
" }\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_ifft_obj, 1, 2, fft_ifft);\n",
"\n",
"mp_obj_t fft_spectrum(size_t n_args, const mp_obj_t *args) {\n",
" if(n_args == 2) {\n",
" return fft_fft_ifft_spectrum(n_args, args[0], args[1], FFT_SPECTRUM);\n",
" } else {\n",
" return fft_fft_ifft_spectrum(n_args, args[0], mp_const_none, FFT_SPECTRUM);\n",
" }\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_spectrum_obj, 1, 2, fft_spectrum);\n",
"\n",
"STATIC const mp_rom_map_elem_t ulab_fft_globals_table[] = {\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_fft) },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_fft), (mp_obj_t)&fft_fft_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_ifft), (mp_obj_t)&fft_ifft_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_spectrum), (mp_obj_t)&fft_spectrum_obj },\n",
"};\n",
"\n",
"STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_fft_globals, ulab_fft_globals_table);\n",
"\n",
"mp_obj_module_t ulab_fft_module = {\n",
" .base = { &mp_type_module },\n",
" .globals = (mp_obj_dict_t*)&mp_module_ulab_fft_globals,\n",
"};\n",
"\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerical\n",
"\n",
"## General comments\n",
"\n",
"This section contains miscellaneous functions that did not fit in the other submodules. These include `linspace`, `min/max`, `argmin/argmax`, `sum`, `mean`, `std`. These latter functions work with iterables, or ndarrays. When the ndarray is two-dimensional, an `axis` keyword can be supplied, in which case, the function returns a vector, otherwise a scalar.\n",
"\n",
"Since the return values of `mean`, and `std` are most probably floats, these functions return ndarrays of type float, while `min/max` and `clip` do not change the type, and `argmin/argmax` return `uint8`, if the values are smaller than 255, otherwise, `uint16`.\n",
"\n",
"### roll\n",
"\n",
"Note that at present, arrays are always rolled to the left, even when the user specifies right. The reason for that is inner working of `memcpy`: one can shift contiguous chunks to the left only. If one tries to shift to the right, then the same value will be written into the new array over and over again."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## numerical.h"
]
},
{
"cell_type": "code",
"execution_count": 191,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-13T20:46:16.833471Z",
"start_time": "2020-02-13T20:46:16.825605Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 4904 bytes to numerical.h\n"
]
}
],
"source": [
"%%ccode numerical.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#ifndef _NUMERICAL_\n",
"#define _NUMERICAL_\n",
"\n",
"#include \"ulab.h\"\n",
"#include \"ndarray.h\"\n",
"\n",
"#if ULAB_NUMERICAL_MODULE\n",
"\n",
"mp_obj_module_t ulab_numerical_module;\n",
"\n",
"// TODO: implement minimum/maximum, and cumsum\n",
"//mp_obj_t numerical_minimum(mp_obj_t , mp_obj_t );\n",
"//mp_obj_t numerical_maximum(mp_obj_t , mp_obj_t );\n",
"//mp_obj_t numerical_cumsum(size_t , const mp_obj_t *, mp_map_t *);\n",
"\n",
"#define RUN_ARGMIN(in, out, typein, typeout, len, start, increment, op, pos) do {\\\n",
" typein *array = (typein *)(in)->array->items;\\\n",
" typeout *outarray = (typeout *)(out)->array->items;\\\n",
" size_t best_index = 0;\\\n",
" if(((op) == NUMERICAL_MAX) || ((op) == NUMERICAL_ARGMAX)) {\\\n",
" for(size_t i=1; i < (len); i++) {\\\n",
" if(array[(start)+i*(increment)] > array[(start)+best_index*(increment)]) best_index = i;\\\n",
" }\\\n",
" if((op) == NUMERICAL_MAX) outarray[(pos)] = array[(start)+best_index*(increment)];\\\n",
" else outarray[(pos)] = best_index;\\\n",
" } else{\\\n",
" for(size_t i=1; i < (len); i++) {\\\n",
" if(array[(start)+i*(increment)] < array[(start)+best_index*(increment)]) best_index = i;\\\n",
" }\\\n",
" if((op) == NUMERICAL_MIN) outarray[(pos)] = array[(start)+best_index*(increment)];\\\n",
" else outarray[(pos)] = best_index;\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"#define RUN_SUM(ndarray, type, optype, len, start, increment) do {\\\n",
" type *array = (type *)(ndarray)->array->items;\\\n",
" type value;\\\n",
" for(size_t j=0; j < (len); j++) {\\\n",
" value = array[(start)+j*(increment)];\\\n",
" sum += value;\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"#define RUN_STD(ndarray, type, len, start, increment) do {\\\n",
" type *array = (type *)(ndarray)->array->items;\\\n",
" mp_float_t value;\\\n",
" for(size_t j=0; j < (len); j++) {\\\n",
" sum += array[(start)+j*(increment)];\\\n",
" }\\\n",
" sum /= (len);\\\n",
" for(size_t j=0; j < (len); j++) {\\\n",
" value = (array[(start)+j*(increment)] - sum);\\\n",
" sum_sq += value * value;\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"#define CALCULATE_DIFF(in, out, type, M, N, inn, increment) do {\\\n",
" type *source = (type *)(in)->array->items;\\\n",
" type *target = (type *)(out)->array->items;\\\n",
" for(size_t i=0; i < (M); i++) {\\\n",
" for(size_t j=0; j < (N); j++) {\\\n",
" for(uint8_t k=0; k < n+1; k++) {\\\n",
" target[i*(N)+j] -= stencil[k]*source[i*(inn)+j+k*(increment)];\\\n",
" }\\\n",
" }\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"#define HEAPSORT(type, ndarray) do {\\\n",
" type *array = (type *)(ndarray)->array->items;\\\n",
" type tmp;\\\n",
" for (;;) {\\\n",
" if (k > 0) {\\\n",
" tmp = array[start+(--k)*increment];\\\n",
" } else {\\\n",
" q--;\\\n",
" if(q == 0) {\\\n",
" break;\\\n",
" }\\\n",
" tmp = array[start+q*increment];\\\n",
" array[start+q*increment] = array[start];\\\n",
" }\\\n",
" p = k;\\\n",
" c = k + k + 1;\\\n",
" while (c < q) {\\\n",
" if((c + 1 < q) && (array[start+(c+1)*increment] > array[start+c*increment])) {\\\n",
" c++;\\\n",
" }\\\n",
" if(array[start+c*increment] > tmp) {\\\n",
" array[start+p*increment] = array[start+c*increment];\\\n",
" p = c;\\\n",
" c = p + p + 1;\\\n",
" } else {\\\n",
" break;\\\n",
" }\\\n",
" }\\\n",
" array[start+p*increment] = tmp;\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"// This is pretty similar to HEAPSORT above; perhaps, the two could be combined somehow\n",
"// On the other hand, since this is a macro, it doesn't really matter\n",
"// Keep in mind that initially, index_array[start+s*increment] = s\n",
"#define HEAP_ARGSORT(type, ndarray, index_array) do {\\\n",
" type *array = (type *)(ndarray)->array->items;\\\n",
" type tmp;\\\n",
" uint16_t itmp;\\\n",
" for (;;) {\\\n",
" if (k > 0) {\\\n",
" k--;\\\n",
" tmp = array[start+index_array[start+k*increment]*increment];\\\n",
" itmp = index_array[start+k*increment];\\\n",
" } else {\\\n",
" q--;\\\n",
" if(q == 0) {\\\n",
" break;\\\n",
" }\\\n",
" tmp = array[start+index_array[start+q*increment]*increment];\\\n",
" itmp = index_array[start+q*increment];\\\n",
" index_array[start+q*increment] = index_array[start];\\\n",
" }\\\n",
" p = k;\\\n",
" c = k + k + 1;\\\n",
" while (c < q) {\\\n",
" if((c + 1 < q) && (array[start+index_array[start+(c+1)*increment]*increment] > array[start+index_array[start+c*increment]*increment])) {\\\n",
" c++;\\\n",
" }\\\n",
" if(array[start+index_array[start+c*increment]*increment] > tmp) {\\\n",
" index_array[start+p*increment] = index_array[start+c*increment];\\\n",
" p = c;\\\n",
" c = p + p + 1;\\\n",
" } else {\\\n",
" break;\\\n",
" }\\\n",
" }\\\n",
" index_array[start+p*increment] = itmp;\\\n",
" }\\\n",
"} while(0)\n",
"\n",
"#endif\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## numerical.c\n",
"\n",
"### Argument parsing\n",
"\n",
"\n",
"Since most of these functions operate on matrices along an axis, it might make sense to factor out the parsing of arguments and keyword arguments. The void function `numerical_parse_args` fills in the pointer for the matrix/array, and the axis."
]
},
{
"cell_type": "code",
"execution_count": 188,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-13T20:45:21.296834Z",
"start_time": "2020-02-13T20:45:21.286578Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 30538 bytes to numerical.c\n"
]
}
],
"source": [
"%%ccode numerical.c\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project, \n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#include <math.h>\n",
"#include <stdlib.h>\n",
"#include <string.h>\n",
"#include \"py/obj.h\"\n",
"#include \"py/objint.h\"\n",
"#include \"py/runtime.h\"\n",
"#include \"py/builtin.h\"\n",
"#include \"py/misc.h\"\n",
"#include \"numerical.h\"\n",
"\n",
"#if ULAB_NUMERICAL_MODULE\n",
"\n",
"enum NUMERICAL_FUNCTION_TYPE {\n",
" NUMERICAL_MIN,\n",
" NUMERICAL_MAX,\n",
" NUMERICAL_ARGMIN,\n",
" NUMERICAL_ARGMAX,\n",
" NUMERICAL_SUM,\n",
" NUMERICAL_MEAN,\n",
" NUMERICAL_STD,\n",
"};\n",
"\n",
"mp_obj_t numerical_linspace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_num, MP_ARG_INT, {.u_int = 50} },\n",
" { MP_QSTR_endpoint, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_true} },\n",
" { MP_QSTR_retstep, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_false} },\n",
" { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(2, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
"\n",
" uint16_t len = args[2].u_int;\n",
" if(len < 2) {\n",
" mp_raise_ValueError(translate(\"number of points must be at least 2\"));\n",
" }\n",
" mp_float_t value, step;\n",
" value = mp_obj_get_float(args[0].u_obj);\n",
" uint8_t typecode = args[5].u_int;\n",
" if(args[3].u_obj == mp_const_true) step = (mp_obj_get_float(args[1].u_obj)-value)/(len-1);\n",
" else step = (mp_obj_get_float(args[1].u_obj)-value)/len;\n",
" ndarray_obj_t *ndarray = create_new_ndarray(1, len, typecode);\n",
" if(typecode == NDARRAY_UINT8) {\n",
" uint8_t *array = (uint8_t *)ndarray->array->items;\n",
" for(size_t i=0; i < len; i++, value += step) array[i] = (uint8_t)value;\n",
" } else if(typecode == NDARRAY_INT8) {\n",
" int8_t *array = (int8_t *)ndarray->array->items;\n",
" for(size_t i=0; i < len; i++, value += step) array[i] = (int8_t)value;\n",
" } else if(typecode == NDARRAY_UINT16) {\n",
" uint16_t *array = (uint16_t *)ndarray->array->items;\n",
" for(size_t i=0; i < len; i++, value += step) array[i] = (uint16_t)value;\n",
" } else if(typecode == NDARRAY_INT16) {\n",
" int16_t *array = (int16_t *)ndarray->array->items;\n",
" for(size_t i=0; i < len; i++, value += step) array[i] = (int16_t)value;\n",
" } else {\n",
" mp_float_t *array = (mp_float_t *)ndarray->array->items;\n",
" for(size_t i=0; i < len; i++, value += step) array[i] = value;\n",
" }\n",
" if(args[4].u_obj == mp_const_false) {\n",
" return MP_OBJ_FROM_PTR(ndarray);\n",
" } else {\n",
" mp_obj_t tuple[2];\n",
" tuple[0] = ndarray;\n",
" tuple[1] = mp_obj_new_float(step);\n",
" return mp_obj_new_tuple(2, tuple);\n",
" }\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_linspace_obj, 2, numerical_linspace);\n",
"\n",
"void axis_sorter(ndarray_obj_t *ndarray, mp_obj_t axis, size_t *m, size_t *n, size_t *N, \n",
" size_t *increment, size_t *len, size_t *start_inc) {\n",
" if(axis == mp_const_none) { // flatten the array\n",
" *m = 1;\n",
" *n = 1;\n",
" *len = ndarray->array->len;\n",
" *N = 1;\n",
" *increment = 1;\n",
" *start_inc = ndarray->array->len;\n",
" } else if((mp_obj_get_int(axis) == 1)) { // along the horizontal axis\n",
" *m = ndarray->m;\n",
" *n = 1;\n",
" *len = ndarray->n;\n",
" *N = ndarray->m;\n",
" *increment = 1;\n",
" *start_inc = ndarray->n;\n",
" } else { // along vertical axis\n",
" *m = 1;\n",
" *n = ndarray->n;\n",
" *len = ndarray->m;\n",
" *N = ndarray->n;\n",
" *increment = ndarray->n;\n",
" *start_inc = 1;\n",
" } \n",
"}\n",
"\n",
"mp_obj_t numerical_sum_mean_std_iterable(mp_obj_t oin, uint8_t optype, size_t ddof) {\n",
" mp_float_t value, sum = 0.0, sq_sum = 0.0;\n",
" mp_obj_iter_buf_t iter_buf;\n",
" mp_obj_t item, iterable = mp_getiter(oin, &iter_buf);\n",
" mp_int_t len = mp_obj_get_int(mp_obj_len(oin));\n",
" while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" value = mp_obj_get_float(item);\n",
" sum += value;\n",
" }\n",
" if(optype == NUMERICAL_SUM) {\n",
" return mp_obj_new_float(sum);\n",
" } else if(optype == NUMERICAL_MEAN) {\n",
" return mp_obj_new_float(sum/len);\n",
" } else { // this should be the case of the standard deviation\n",
" // TODO: note that we could get away with a single pass, if we used the Weldorf algorithm\n",
" // That should save a fair amount of time, because we would have to extract the values only once\n",
" iterable = mp_getiter(oin, &iter_buf);\n",
" sum /= len; // this is now the mean!\n",
" while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" value = mp_obj_get_float(item) - sum;\n",
" sq_sum += value * value;\n",
" }\n",
" return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(sq_sum/(len-ddof)));\n",
" }\n",
"}\n",
"\n",
"STATIC mp_obj_t numerical_sum_mean_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {\n",
" size_t m, n, increment, start, start_inc, N, len; \n",
" axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);\n",
" ndarray_obj_t *results = create_new_ndarray(m, n, NDARRAY_FLOAT);\n",
" mp_float_t sum, sq_sum;\n",
" mp_float_t *farray = (mp_float_t *)results->array->items;\n",
" for(size_t j=0; j < N; j++) { // result index\n",
" start = j * start_inc;\n",
" sum = sq_sum = 0.0;\n",
" if(ndarray->array->typecode == NDARRAY_UINT8) {\n",
" RUN_SUM(ndarray, uint8_t, optype, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_INT8) {\n",
" RUN_SUM(ndarray, int8_t, optype, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_UINT16) {\n",
" RUN_SUM(ndarray, uint16_t, optype, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_INT16) {\n",
" RUN_SUM(ndarray, int16_t, optype, len, start, increment);\n",
" } else { // this will be mp_float_t, no need to check\n",
" RUN_SUM(ndarray, mp_float_t, optype, len, start, increment);\n",
" }\n",
" if(optype == NUMERICAL_SUM) {\n",
" farray[j] = sum;\n",
" } else { // this is the case of the mean\n",
" farray[j] = sum / len;\n",
" }\n",
" }\n",
" if(results->array->len == 1) {\n",
" return mp_obj_new_float(farray[0]);\n",
" }\n",
" return MP_OBJ_FROM_PTR(results);\n",
"}\n",
"\n",
"mp_obj_t numerical_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, size_t ddof) {\n",
" size_t m, n, increment, start, start_inc, N, len; \n",
" mp_float_t sum, sum_sq;\n",
" \n",
" axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);\n",
" if(ddof > len) {\n",
" mp_raise_ValueError(translate(\"ddof must be smaller than length of data set\"));\n",
" }\n",
" ndarray_obj_t *results = create_new_ndarray(m, n, NDARRAY_FLOAT);\n",
" mp_float_t *farray = (mp_float_t *)results->array->items;\n",
" for(size_t j=0; j < N; j++) { // result index\n",
" start = j * start_inc;\n",
" sum = 0.0;\n",
" sum_sq = 0.0;\n",
" if(ndarray->array->typecode == NDARRAY_UINT8) {\n",
" RUN_STD(ndarray, uint8_t, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_INT8) {\n",
" RUN_STD(ndarray, int8_t, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_UINT16) {\n",
" RUN_STD(ndarray, uint16_t, len, start, increment);\n",
" } else if(ndarray->array->typecode == NDARRAY_INT16) {\n",
" RUN_STD(ndarray, int16_t, len, start, increment);\n",
" } else { // this will be mp_float_t, no need to check\n",
" RUN_STD(ndarray, mp_float_t, len, start, increment);\n",
" }\n",
" farray[j] = MICROPY_FLOAT_C_FUN(sqrt)(sum_sq/(len - ddof));\n",
" }\n",
" if(results->array->len == 1) {\n",
" return mp_obj_new_float(farray[0]);\n",
" }\n",
" return MP_OBJ_FROM_PTR(results);\n",
"}\n",
"\n",
"mp_obj_t numerical_argmin_argmax_iterable(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {\n",
" size_t idx = 0, best_idx = 0;\n",
" mp_obj_iter_buf_t iter_buf;\n",
" mp_obj_t iterable = mp_getiter(oin, &iter_buf);\n",
" mp_obj_t best_obj = MP_OBJ_NULL;\n",
" mp_obj_t item;\n",
" mp_uint_t op = MP_BINARY_OP_LESS;\n",
" if((optype == NUMERICAL_ARGMAX) || (optype == NUMERICAL_MAX)) op = MP_BINARY_OP_MORE;\n",
" while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {\n",
" if ((best_obj == MP_OBJ_NULL) || (mp_binary_op(op, item, best_obj) == mp_const_true)) {\n",
" best_obj = item;\n",
" best_idx = idx;\n",
" }\n",
" idx++;\n",
" }\n",
" if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {\n",
" return MP_OBJ_NEW_SMALL_INT(best_idx);\n",
" } else {\n",
" return best_obj;\n",
" } \n",
"}\n",
"\n",
"mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {\n",
" size_t m, n, increment, start, start_inc, N, len;\n",
" axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);\n",
" ndarray_obj_t *results;\n",
" if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {\n",
" // we could save some RAM by taking NDARRAY_UINT8, if the dimensions \n",
" // are smaller than 256, but the code would become more involving \n",
" // (we would also need extra flash space)\n",
" results = create_new_ndarray(m, n, NDARRAY_UINT16);\n",
" } else {\n",
" results = create_new_ndarray(m, n, ndarray->array->typecode);\n",
" }\n",
" \n",
" for(size_t j=0; j < N; j++) { // result index\n",
" start = j * start_inc;\n",
" if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {\n",
" if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {\n",
" RUN_ARGMIN(ndarray, results, uint8_t, uint8_t, len, start, increment, optype, j);\n",
" } else {\n",
" RUN_ARGMIN(ndarray, results, uint8_t, uint16_t, len, start, increment, optype, j); \n",
" }\n",
" } else if((ndarray->array->typecode == NDARRAY_UINT16) || (ndarray->array->typecode == NDARRAY_INT16)) {\n",
" RUN_ARGMIN(ndarray, results, uint16_t, uint16_t, len, start, increment, optype, j);\n",
" } else {\n",
" if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {\n",
" RUN_ARGMIN(ndarray, results, mp_float_t, mp_float_t, len, start, increment, optype, j);\n",
" } else {\n",
" RUN_ARGMIN(ndarray, results, mp_float_t, uint16_t, len, start, increment, optype, j); \n",
" }\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(results);\n",
"}\n",
"\n",
"STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t optype) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} } ,\n",
" { MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" \n",
" mp_obj_t oin = args[0].u_obj;\n",
" mp_obj_t axis = args[1].u_obj;\n",
" if((axis != mp_const_none) && (mp_obj_get_int(axis) != 0) && (mp_obj_get_int(axis) != 1)) {\n",
" // this seems to pass with False, and True...\n",
" mp_raise_ValueError(translate(\"axis must be None, 0, or 1\"));\n",
" }\n",
" \n",
" if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) || \n",
" MP_OBJ_IS_TYPE(oin, &mp_type_range)) {\n",
" switch(optype) {\n",
" case NUMERICAL_MIN:\n",
" case NUMERICAL_ARGMIN:\n",
" case NUMERICAL_MAX:\n",
" case NUMERICAL_ARGMAX:\n",
" return numerical_argmin_argmax_iterable(oin, axis, optype);\n",
" case NUMERICAL_SUM:\n",
" case NUMERICAL_MEAN:\n",
" return numerical_sum_mean_std_iterable(oin, optype, 0);\n",
" default: // we should never reach this point, but whatever\n",
" return mp_const_none;\n",
" }\n",
" } else if(MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {\n",
" ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);\n",
" switch(optype) {\n",
" case NUMERICAL_MIN:\n",
" case NUMERICAL_MAX:\n",
" case NUMERICAL_ARGMIN:\n",
" case NUMERICAL_ARGMAX:\n",
" return numerical_argmin_argmax_ndarray(ndarray, axis, optype);\n",
" case NUMERICAL_SUM:\n",
" case NUMERICAL_MEAN:\n",
" return numerical_sum_mean_ndarray(ndarray, axis, optype);\n",
" default:\n",
" mp_raise_NotImplementedError(translate(\"operation is not implemented on ndarrays\"));\n",
" }\n",
" } else {\n",
" mp_raise_TypeError(translate(\"input must be tuple, list, range, or ndarray\"));\n",
" }\n",
" return mp_const_none;\n",
"}\n",
"\n",
"mp_obj_t numerical_min(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MIN);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_min_obj, 1, numerical_min);\n",
"\n",
"mp_obj_t numerical_max(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MAX);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_max_obj, 1, numerical_max);\n",
"\n",
"mp_obj_t numerical_argmin(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" return numerical_function(n_args, pos_args, kw_args, NUMERICAL_ARGMIN);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argmin_obj, 1, numerical_argmin);\n",
"\n",
"mp_obj_t numerical_argmax(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" return numerical_function(n_args, pos_args, kw_args, NUMERICAL_ARGMAX);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argmax_obj, 1, numerical_argmax);\n",
"\n",
"mp_obj_t numerical_sum(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" return numerical_function(n_args, pos_args, kw_args, NUMERICAL_SUM);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sum_obj, 1, numerical_sum);\n",
"\n",
"mp_obj_t numerical_mean(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MEAN);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_mean_obj, 1, numerical_mean);\n",
"\n",
"mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } } ,\n",
" { MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_ddof, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" \n",
" mp_obj_t oin = args[0].u_obj;\n",
" mp_obj_t axis = args[1].u_obj;\n",
" size_t ddof = args[2].u_int;\n",
" if((axis != mp_const_none) && (mp_obj_get_int(axis) != 0) && (mp_obj_get_int(axis) != 1)) {\n",
" // this seems to pass with False, and True...\n",
" mp_raise_ValueError(translate(\"axis must be None, 0, or 1\"));\n",
" }\n",
" if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) || MP_OBJ_IS_TYPE(oin, &mp_type_range)) {\n",
" return numerical_sum_mean_std_iterable(oin, NUMERICAL_STD, ddof);\n",
" } else if(MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {\n",
" ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(oin);\n",
" return numerical_std_ndarray(ndarray, axis, ddof);\n",
" } else {\n",
" mp_raise_TypeError(translate(\"input must be tuple, list, range, or ndarray\"));\n",
" }\n",
" return mp_const_none;\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_std_obj, 1, numerical_std);\n",
"\n",
"mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(2, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" \n",
" mp_obj_t oin = args[0].u_obj;\n",
" int16_t shift = mp_obj_get_int(args[1].u_obj);\n",
" if((args[2].u_obj != mp_const_none) && \n",
" (mp_obj_get_int(args[2].u_obj) != 0) && \n",
" (mp_obj_get_int(args[2].u_obj) != 1)) {\n",
" mp_raise_ValueError(translate(\"axis must be None, 0, or 1\"));\n",
" }\n",
"\n",
" ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);\n",
" uint8_t _sizeof = mp_binary_get_size('@', in->array->typecode, NULL);\n",
" size_t len;\n",
" int16_t _shift;\n",
" uint8_t *array = (uint8_t *)in->array->items;\n",
" // TODO: transpose the matrix, if axis == 0. Though, that is hard on the RAM...\n",
" if(shift < 0) {\n",
" _shift = -shift;\n",
" } else {\n",
" _shift = shift;\n",
" }\n",
" if((args[2].u_obj == mp_const_none) || (mp_obj_get_int(args[2].u_obj) == 1)) { // shift horizontally\n",
" uint16_t M;\n",
" if(args[2].u_obj == mp_const_none) {\n",
" len = in->array->len;\n",
" M = 1;\n",
" } else {\n",
" len = in->n;\n",
" M = in->m;\n",
" }\n",
" _shift = _shift % len;\n",
" if(shift < 0) _shift = len - _shift;\n",
" // TODO: if(shift > len/2), we should move in the opposite direction. That would save RAM\n",
" _shift *= _sizeof;\n",
" uint8_t *tmp = m_new(uint8_t, _shift);\n",
" for(size_t m=0; m < M; m++) {\n",
" memmove(tmp, &array[m*len*_sizeof], _shift);\n",
" memmove(&array[m*len*_sizeof], &array[m*len*_sizeof+_shift], len*_sizeof-_shift);\n",
" memmove(&array[(m+1)*len*_sizeof-_shift], tmp, _shift);\n",
" }\n",
" m_del(uint8_t, tmp, _shift);\n",
" return mp_const_none;\n",
" } else {\n",
" len = in->m;\n",
" // temporary buffer\n",
" uint8_t *_data = m_new(uint8_t, _sizeof*len);\n",
" \n",
" _shift = _shift % len;\n",
" if(shift < 0) _shift = len - _shift;\n",
" _shift *= _sizeof;\n",
" uint8_t *tmp = m_new(uint8_t, _shift);\n",
"\n",
" for(size_t n=0; n < in->n; n++) {\n",
" for(size_t m=0; m < len; m++) {\n",
" // this loop should fill up the temporary buffer\n",
" memmove(&_data[m*_sizeof], &array[(m*in->n+n)*_sizeof], _sizeof);\n",
" }\n",
" // now, the actual shift\n",
" memmove(tmp, _data, _shift);\n",
" memmove(_data, &_data[_shift], len*_sizeof-_shift);\n",
" memmove(&_data[len*_sizeof-_shift], tmp, _shift);\n",
" for(size_t m=0; m < len; m++) {\n",
" // this loop should dump the content of the temporary buffer into data\n",
" memmove(&array[(m*in->n+n)*_sizeof], &_data[m*_sizeof], _sizeof);\n",
" } \n",
" }\n",
" m_del(uint8_t, tmp, _shift);\n",
" m_del(uint8_t, _data, _sizeof*len);\n",
" return mp_const_none;\n",
" }\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_roll_obj, 2, numerical_roll);\n",
"\n",
"mp_obj_t numerical_flip(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" \n",
" if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"flip argument must be an ndarray\"));\n",
" }\n",
" if((args[1].u_obj != mp_const_none) && \n",
" (mp_obj_get_int(args[1].u_obj) != 0) && \n",
" (mp_obj_get_int(args[1].u_obj) != 1)) {\n",
" mp_raise_ValueError(translate(\"axis must be None, 0, or 1\"));\n",
" }\n",
"\n",
" ndarray_obj_t *in = MP_OBJ_TO_PTR(args[0].u_obj);\n",
" mp_obj_t oout = ndarray_copy(args[0].u_obj);\n",
" ndarray_obj_t *out = MP_OBJ_TO_PTR(oout);\n",
" uint8_t _sizeof = mp_binary_get_size('@', in->array->typecode, NULL);\n",
" uint8_t *array_in = (uint8_t *)in->array->items;\n",
" uint8_t *array_out = (uint8_t *)out->array->items; \n",
" size_t len;\n",
" if((args[1].u_obj == mp_const_none) || (mp_obj_get_int(args[1].u_obj) == 1)) { // flip horizontally\n",
" uint16_t M = in->m;\n",
" len = in->n;\n",
" if(args[1].u_obj == mp_const_none) { // flip flattened array\n",
" len = in->array->len;\n",
" M = 1;\n",
" }\n",
" for(size_t m=0; m < M; m++) {\n",
" for(size_t n=0; n < len; n++) {\n",
" memcpy(array_out+_sizeof*(m*len+n), array_in+_sizeof*((m+1)*len-n-1), _sizeof);\n",
" }\n",
" }\n",
" } else { // flip vertically\n",
" for(size_t m=0; m < in->m; m++) {\n",
" for(size_t n=0; n < in->n; n++) {\n",
" memcpy(array_out+_sizeof*(m*in->n+n), array_in+_sizeof*((in->m-m-1)*in->n+n), _sizeof);\n",
" }\n",
" }\n",
" }\n",
" return out;\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_flip_obj, 1, numerical_flip);\n",
"\n",
"mp_obj_t numerical_diff(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_n, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 1 } },\n",
" { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = -1 } },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" \n",
" if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"diff argument must be an ndarray\"));\n",
" }\n",
" \n",
" ndarray_obj_t *in = MP_OBJ_TO_PTR(args[0].u_obj);\n",
" size_t increment, N, M;\n",
" if((args[2].u_int == -1) || (args[2].u_int == 1)) { // differentiate along the horizontal axis\n",
" increment = 1;\n",
" } else if(args[2].u_int == 0) { // differtiate along vertical axis\n",
" increment = in->n;\n",
" } else {\n",
" mp_raise_ValueError(translate(\"axis must be -1, 0, or 1\"));\n",
" }\n",
" if((args[1].u_int < 0) || (args[1].u_int > 9)) {\n",
" mp_raise_ValueError(translate(\"n must be between 0, and 9\"));\n",
" }\n",
" uint8_t n = args[1].u_int;\n",
" int8_t *stencil = m_new(int8_t, n+1);\n",
" stencil[0] = 1;\n",
" for(uint8_t i=1; i < n+1; i++) {\n",
" stencil[i] = -stencil[i-1]*(n-i+1)/i;\n",
" }\n",
"\n",
" ndarray_obj_t *out;\n",
" \n",
" if(increment == 1) { // differentiate along the horizontal axis \n",
" if(n >= in->n) {\n",
" out = create_new_ndarray(in->m, 0, in->array->typecode);\n",
" m_del(uint8_t, stencil, n);\n",
" return MP_OBJ_FROM_PTR(out);\n",
" }\n",
" N = in->n - n;\n",
" M = in->m;\n",
" } else { // differentiate along vertical axis\n",
" if(n >= in->m) {\n",
" out = create_new_ndarray(0, in->n, in->array->typecode);\n",
" m_del(uint8_t, stencil, n);\n",
" return MP_OBJ_FROM_PTR(out);\n",
" }\n",
" M = in->m - n;\n",
" N = in->n;\n",
" }\n",
" out = create_new_ndarray(M, N, in->array->typecode);\n",
" if(in->array->typecode == NDARRAY_UINT8) {\n",
" CALCULATE_DIFF(in, out, uint8_t, M, N, in->n, increment);\n",
" } else if(in->array->typecode == NDARRAY_INT8) {\n",
" CALCULATE_DIFF(in, out, int8_t, M, N, in->n, increment);\n",
" } else if(in->array->typecode == NDARRAY_UINT16) {\n",
" CALCULATE_DIFF(in, out, uint16_t, M, N, in->n, increment);\n",
" } else if(in->array->typecode == NDARRAY_INT16) {\n",
" CALCULATE_DIFF(in, out, int16_t, M, N, in->n, increment);\n",
" } else {\n",
" CALCULATE_DIFF(in, out, mp_float_t, M, N, in->n, increment);\n",
" }\n",
" m_del(int8_t, stencil, n);\n",
" return MP_OBJ_FROM_PTR(out);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_diff_obj, 1, numerical_diff);\n",
"\n",
"mp_obj_t numerical_sort_helper(mp_obj_t oin, mp_obj_t axis, uint8_t inplace) {\n",
" if(!MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"sort argument must be an ndarray\"));\n",
" }\n",
"\n",
" ndarray_obj_t *ndarray;\n",
" mp_obj_t out;\n",
" if(inplace == 1) {\n",
" ndarray = MP_OBJ_TO_PTR(oin);\n",
" } else {\n",
" out = ndarray_copy(oin);\n",
" ndarray = MP_OBJ_TO_PTR(out);\n",
" }\n",
" size_t increment, start_inc, end, N;\n",
" if(axis == mp_const_none) { // flatten the array\n",
" ndarray->m = 1;\n",
" ndarray->n = ndarray->array->len;\n",
" increment = 1;\n",
" start_inc = ndarray->n;\n",
" end = ndarray->n;\n",
" N = ndarray->n;\n",
" } else if((mp_obj_get_int(axis) == -1) || \n",
" (mp_obj_get_int(axis) == 1)) { // sort along the horizontal axis\n",
" increment = 1;\n",
" start_inc = ndarray->n;\n",
" end = ndarray->array->len;\n",
" N = ndarray->n;\n",
" } else if(mp_obj_get_int(axis) == 0) { // sort along vertical axis\n",
" increment = ndarray->n;\n",
" start_inc = 1;\n",
" end = ndarray->m;\n",
" N = ndarray->m;\n",
" } else {\n",
" mp_raise_ValueError(translate(\"axis must be -1, 0, None, or 1\"));\n",
" }\n",
" \n",
" size_t q, k, p, c;\n",
"\n",
" for(size_t start=0; start < end; start+=start_inc) {\n",
" q = N; \n",
" k = (q >> 1);\n",
" if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {\n",
" HEAPSORT(uint8_t, ndarray);\n",
" } else if((ndarray->array->typecode == NDARRAY_INT16) || (ndarray->array->typecode == NDARRAY_INT16)) {\n",
" HEAPSORT(uint16_t, ndarray);\n",
" } else {\n",
" HEAPSORT(mp_float_t, ndarray);\n",
" }\n",
" }\n",
" if(inplace == 1) {\n",
" return mp_const_none;\n",
" } else {\n",
" return out;\n",
" }\n",
"}\n",
"\n",
"// numpy function\n",
"mp_obj_t numerical_sort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_int = -1 } },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
"\n",
" return numerical_sort_helper(args[0].u_obj, args[1].u_obj, 0);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sort_obj, 1, numerical_sort);\n",
"\n",
"// method of an ndarray\n",
"mp_obj_t numerical_sort_inplace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_int = -1 } },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
"\n",
" return numerical_sort_helper(args[0].u_obj, args[1].u_obj, 1);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sort_inplace_obj, 1, numerical_sort_inplace);\n",
"\n",
"mp_obj_t numerical_argsort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_int = -1 } },\n",
" };\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
" if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"argsort argument must be an ndarray\"));\n",
" }\n",
"\n",
" ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);\n",
" size_t increment, start_inc, end, N, m, n;\n",
" if(args[1].u_obj == mp_const_none) { // flatten the array\n",
" m = 1;\n",
" n = ndarray->array->len;\n",
" ndarray->m = m;\n",
" ndarray->n = n;\n",
" increment = 1;\n",
" start_inc = ndarray->n;\n",
" end = ndarray->n;\n",
" N = n;\n",
" } else if((mp_obj_get_int(args[1].u_obj) == -1) || \n",
" (mp_obj_get_int(args[1].u_obj) == 1)) { // sort along the horizontal axis\n",
" m = ndarray->m;\n",
" n = ndarray->n;\n",
" increment = 1;\n",
" start_inc = n;\n",
" end = ndarray->array->len;\n",
" N = n;\n",
" } else if(mp_obj_get_int(args[1].u_obj) == 0) { // sort along vertical axis\n",
" m = ndarray->m;\n",
" n = ndarray->n;\n",
" increment = n;\n",
" start_inc = 1;\n",
" end = m;\n",
" N = m;\n",
" } else {\n",
" mp_raise_ValueError(translate(\"axis must be -1, 0, None, or 1\"));\n",
" }\n",
"\n",
" // at the expense of flash, we could save RAM by creating \n",
" // an NDARRAY_UINT16 ndarray only, if needed, otherwise, NDARRAY_UINT8\n",
" ndarray_obj_t *indices = create_new_ndarray(m, n, NDARRAY_UINT16);\n",
" uint16_t *index_array = (uint16_t *)indices->array->items;\n",
" // initialise the index array\n",
" // if array is flat: 0 to indices->n\n",
" // if sorting vertically, identical indices are arranged row-wise\n",
" // if sorting horizontally, identical indices are arranged colunn-wise\n",
" for(uint16_t start=0; start < end; start+=start_inc) {\n",
" for(uint16_t s=0; s < N; s++) {\n",
" index_array[start+s*increment] = s;\n",
" }\n",
" }\n",
"\n",
" size_t q, k, p, c;\n",
" for(size_t start=0; start < end; start+=start_inc) {\n",
" q = N; \n",
" k = (q >> 1);\n",
" if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {\n",
" HEAP_ARGSORT(uint8_t, ndarray, index_array);\n",
" } else if((ndarray->array->typecode == NDARRAY_INT16) || (ndarray->array->typecode == NDARRAY_INT16)) {\n",
" HEAP_ARGSORT(uint16_t, ndarray, index_array);\n",
" } else {\n",
" HEAP_ARGSORT(mp_float_t, ndarray, index_array);\n",
" }\n",
" }\n",
" return MP_OBJ_FROM_PTR(indices);\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argsort_obj, 1, numerical_argsort);\n",
"\n",
"STATIC const mp_rom_map_elem_t ulab_numerical_globals_table[] = {\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_linspace), (mp_obj_t)&numerical_linspace_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_sum), (mp_obj_t)&numerical_sum_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_mean), (mp_obj_t)&numerical_mean_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_std), (mp_obj_t)&numerical_std_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_min), (mp_obj_t)&numerical_min_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_max), (mp_obj_t)&numerical_max_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_argmin), (mp_obj_t)&numerical_argmin_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_argmax), (mp_obj_t)&numerical_argmax_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_roll), (mp_obj_t)&numerical_roll_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_flip), (mp_obj_t)&numerical_flip_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_diff), (mp_obj_t)&numerical_diff_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_sort), (mp_obj_t)&numerical_sort_obj },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_argsort), (mp_obj_t)&numerical_argsort_obj }, \n",
"};\n",
"\n",
"STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_numerical_globals, ulab_numerical_globals_table);\n",
"\n",
"mp_obj_module_t ulab_numerical_module = {\n",
" .base = { &mp_type_module },\n",
" .globals = (mp_obj_dict_t*)&mp_module_ulab_numerical_globals,\n",
"};\n",
"\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Filter functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## filter.h"
]
},
{
"cell_type": "code",
"execution_count": 273,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-14T13:53:29.334248Z",
"start_time": "2020-02-14T13:53:29.322553Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 348 bytes to filter.h\n"
]
}
],
"source": [
"%%ccode filter.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project,\n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2020 Jeff Epler for Adafruit Industries\n",
"*/\n",
"\n",
"#ifndef _FILTER_\n",
"#define _FILTER_\n",
"\n",
"#include \"ulab.h\"\n",
"#include \"ndarray.h\"\n",
"\n",
"#if ULAB_FILTER_MODULE\n",
"\n",
"mp_obj_module_t ulab_filter_module;\n",
"\n",
"#endif\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## filter.c"
]
},
{
"cell_type": "code",
"execution_count": 285,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-14T14:01:12.987357Z",
"start_time": "2020-02-14T14:01:12.982177Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 3506 bytes to filter.c\n"
]
}
],
"source": [
"%%ccode filter.c\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project,\n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2020 Jeff Epler for Adafruit Industries\n",
"*/\n",
"\n",
"#include <math.h>\n",
"#include <stdlib.h>\n",
"#include <string.h>\n",
"#include \"py/obj.h\"\n",
"#include \"py/runtime.h\"\n",
"#include \"py/misc.h\"\n",
"#include \"filter.h\"\n",
"\n",
"#if ULAB_FILTER_MODULE\n",
"mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {\n",
" static const mp_arg_t allowed_args[] = {\n",
" { MP_QSTR_a, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" { MP_QSTR_v, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },\n",
" };\n",
"\n",
" mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];\n",
" mp_arg_parse_all(2, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);\n",
"\n",
" if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type) || !MP_OBJ_IS_TYPE(args[1].u_obj, &ulab_ndarray_type)) {\n",
" mp_raise_TypeError(translate(\"convolve arguments must be ndarrays\"));\n",
" }\n",
"\n",
" ndarray_obj_t *a = MP_OBJ_TO_PTR(args[0].u_obj);\n",
" ndarray_obj_t *c = MP_OBJ_TO_PTR(args[1].u_obj);\n",
" int len_a = a->array->len;\n",
" int len_c = c->array->len;\n",
" // deal with linear arrays only\n",
" if(a->m*a->n != len_a || c->m*c->n != len_c) {\n",
" mp_raise_TypeError(translate(\"convolve arguments must be linear arrays\"));\n",
" }\n",
" if(len_a == 0 || len_c == 0) {\n",
" mp_raise_TypeError(translate(\"convolve arguments must not be empty\"));\n",
" }\n",
"\n",
" int len = len_a + len_c - 1; // convolve mode \"full\"\n",
" ndarray_obj_t *out = create_new_ndarray(1, len, NDARRAY_FLOAT);\n",
" mp_float_t *outptr = out->array->items;\n",
" int off = len_c-1;\n",
"\n",
" if(a->array->typecode == NDARRAY_FLOAT && c->array->typecode == NDARRAY_FLOAT) {\n",
" mp_float_t* a_items = (mp_float_t*)a->array->items;\n",
" mp_float_t* c_items = (mp_float_t*)c->array->items;\n",
" for(int k=-off; k<len-off; k++) {\n",
" mp_float_t accum = (mp_float_t)0;\n",
" int top_n = MIN(len_c, len_a - k);\n",
" int bot_n = MAX(-k, 0);\n",
" mp_float_t* a_ptr = a_items + bot_n + k;\n",
" mp_float_t* a_end = a_ptr + (top_n - bot_n);\n",
" mp_float_t* c_ptr = c_items + len_c - bot_n - 1;\n",
" for(; a_ptr != a_end;) {\n",
" accum += *a_ptr++ * *c_ptr--;\n",
" }\n",
" *outptr++ = accum;\n",
" }\n",
" } else {\n",
" for(int k=-off; k<len-off; k++) {\n",
" mp_float_t accum = (mp_float_t)0;\n",
" int top_n = MIN(len_c, len_a - k);\n",
" int bot_n = MAX(-k, 0);\n",
" for(int n=bot_n; n<top_n; n++) {\n",
" int idx_c = len_c - n - 1;\n",
" int idx_a = n+k;\n",
" mp_float_t ai = ndarray_get_float_value(a->array->items, a->array->typecode, idx_a);\n",
" mp_float_t ci = ndarray_get_float_value(c->array->items, c->array->typecode, idx_c);\n",
" accum += ai * ci;\n",
" }\n",
" *outptr++ = accum;\n",
" }\n",
" }\n",
"\n",
" return out;\n",
"}\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(filter_convolve_obj, 2, filter_convolve);\n",
"\n",
"STATIC const mp_rom_map_elem_t ulab_filter_globals_table[] = {\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_filter) },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_convolve), (mp_obj_t)&filter_convolve_obj },\n",
"};\n",
"\n",
"STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_filter_globals, ulab_filter_globals_table);\n",
"\n",
"mp_obj_module_t ulab_filter_module = {\n",
" .base = { &mp_type_module },\n",
" .globals = (mp_obj_dict_t*)&mp_module_ulab_filter_globals,\n",
"};\n",
"\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Extras\n",
"\n",
"This sub-module is a boilerplate for user-defined extra functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## extras.h"
]
},
{
"cell_type": "code",
"execution_count": 498,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-16T11:28:52.911632Z",
"start_time": "2020-02-16T11:28:52.906281Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 324 bytes to extras.h\n"
]
}
],
"source": [
"%%ccode extras.h\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project,\n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#ifndef _EXTRA_\n",
"#define _EXTRA_\n",
"\n",
"#include \"ulab.h\"\n",
"#include \"ndarray.h\"\n",
"\n",
"#if ULAB_EXTRAS_MODULE\n",
"\n",
"mp_obj_module_t ulab_extras_module;\n",
"\n",
"#endif\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## extras.c"
]
},
{
"cell_type": "code",
"execution_count": 499,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-16T14:40:15.969034Z",
"start_time": "2020-02-16T14:40:15.956352Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 721 bytes to extras.c\n"
]
}
],
"source": [
"%%ccode extras.c\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project,\n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#include <math.h>\n",
"#include <stdlib.h>\n",
"#include <string.h>\n",
"#include \"py/obj.h\"\n",
"#include \"py/runtime.h\"\n",
"#include \"py/misc.h\"\n",
"#include \"extras.h\"\n",
"\n",
"#if ULAB_EXTRAS_MODULE\n",
"\n",
"STATIC const mp_rom_map_elem_t ulab_filter_globals_table[] = {\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_extras) },\n",
"};\n",
"\n",
"STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_extras_globals, ulab_extras_globals_table);\n",
"\n",
"mp_obj_module_t ulab_filter_module = {\n",
" .base = { &mp_type_module },\n",
" .globals = (mp_obj_dict_t*)&mp_module_ulab_extras_globals,\n",
"};\n",
"\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ulab module\n",
"\n",
"This module simply brings all components together, and does not contain new function definitions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ulab.c"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T06:21:51.334700Z",
"start_time": "2020-02-26T06:21:51.329358Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"written 3457 bytes to ulab.c\n"
]
}
],
"source": [
"%%ccode ulab.c\n",
"\n",
"/*\n",
" * This file is part of the micropython-ulab project,\n",
" *\n",
" * https://github.com/v923z/micropython-ulab\n",
" *\n",
" * The MIT License (MIT)\n",
" *\n",
" * Copyright (c) 2019-2020 Zoltán Vörös\n",
"*/\n",
"\n",
"#include <math.h>\n",
"#include <stdio.h>\n",
"#include <stdlib.h>\n",
"#include <string.h>\n",
"#include \"py/runtime.h\"\n",
"#include \"py/binary.h\"\n",
"#include \"py/obj.h\"\n",
"#include \"py/objarray.h\"\n",
"\n",
"#include \"ulab.h\"\n",
"#include \"ndarray.h\"\n",
"#include \"ndarray_properties.h\"\n",
"#include \"linalg.h\"\n",
"#include \"vectorise.h\"\n",
"#include \"poly.h\"\n",
"#include \"fft.h\"\n",
"#include \"filter.h\"\n",
"#include \"numerical.h\"\n",
"#include \"extras.h\"\n",
"\n",
"STATIC MP_DEFINE_STR_OBJ(ulab_version_obj, \"0.33.2\");\n",
"\n",
"MP_DEFINE_CONST_FUN_OBJ_KW(ndarray_flatten_obj, 1, ndarray_flatten);\n",
"\n",
"STATIC const mp_rom_map_elem_t ulab_ndarray_locals_dict_table[] = {\n",
" { MP_ROM_QSTR(MP_QSTR_reshape), MP_ROM_PTR(&ndarray_reshape_obj) },\n",
" { MP_ROM_QSTR(MP_QSTR_transpose), MP_ROM_PTR(&ndarray_transpose_obj) },\n",
" { MP_ROM_QSTR(MP_QSTR_flatten), MP_ROM_PTR(&ndarray_flatten_obj) },\n",
" #if !CIRCUITPY\n",
" { MP_ROM_QSTR(MP_QSTR_shape), MP_ROM_PTR(&ndarray_shape_obj) },\n",
" { MP_ROM_QSTR(MP_QSTR_size), MP_ROM_PTR(&ndarray_size_obj) },\n",
" { MP_ROM_QSTR(MP_QSTR_itemsize), MP_ROM_PTR(&ndarray_itemsize_obj) },\n",
" #endif\n",
"// { MP_ROM_QSTR(MP_QSTR_sort), MP_ROM_PTR(&numerical_sort_inplace_obj) },\n",
"};\n",
"\n",
"STATIC MP_DEFINE_CONST_DICT(ulab_ndarray_locals_dict, ulab_ndarray_locals_dict_table);\n",
"\n",
"const mp_obj_type_t ulab_ndarray_type = {\n",
" { &mp_type_type },\n",
" .name = MP_QSTR_ndarray,\n",
" .print = ndarray_print,\n",
" .make_new = ndarray_make_new,\n",
" .subscr = ndarray_subscr,\n",
" .getiter = ndarray_getiter,\n",
" .unary_op = ndarray_unary_op,\n",
" .binary_op = ndarray_binary_op,\n",
" .buffer_p = { .get_buffer = ndarray_get_buffer, },\n",
" .locals_dict = (mp_obj_dict_t*)&ulab_ndarray_locals_dict,\n",
"};\n",
"\n",
"STATIC const mp_map_elem_t ulab_globals_table[] = {\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_ulab) },\n",
" { MP_ROM_QSTR(MP_QSTR___version__), MP_ROM_PTR(&ulab_version_obj) },\n",
" { MP_OBJ_NEW_QSTR(MP_QSTR_array), (mp_obj_t)&ulab_ndarray_type },\n",
" #if ULAB_LINALG_MODULE\n",
" { MP_ROM_QSTR(MP_QSTR_linalg), MP_ROM_PTR(&ulab_linalg_module) },\n",
" #endif\n",
" #if ULAB_VECTORISE_MODULE\n",
" { MP_ROM_QSTR(MP_QSTR_vector), MP_ROM_PTR(&ulab_vectorise_module) },\n",
" #endif\n",
" #if ULAB_NUMERICAL_MODULE\n",
" { MP_ROM_QSTR(MP_QSTR_numerical), MP_ROM_PTR(&ulab_numerical_module) },\n",
" #endif\n",
" #if ULAB_POLY_MODULE\n",
" { MP_ROM_QSTR(MP_QSTR_poly), MP_ROM_PTR(&ulab_poly_module) },\n",
" #endif\n",
" #if ULAB_FFT_MODULE\n",
" { MP_ROM_QSTR(MP_QSTR_fft), MP_ROM_PTR(&ulab_fft_module) },\n",
" #endif\n",
" #if ULAB_FILTER_MODULE\n",
" { MP_ROM_QSTR(MP_QSTR_filter), MP_ROM_PTR(&ulab_filter_module) },\n",
" #endif\n",
" #if ULAB_EXTRAS_MODULE\n",
" { MP_ROM_QSTR(MP_QSTR_extras), MP_ROM_PTR(&ulab_extras_module) },\n",
" #endif\n",
" // class constants\n",
" { MP_ROM_QSTR(MP_QSTR_uint8), MP_ROM_INT(NDARRAY_UINT8) },\n",
" { MP_ROM_QSTR(MP_QSTR_int8), MP_ROM_INT(NDARRAY_INT8) },\n",
" { MP_ROM_QSTR(MP_QSTR_uint16), MP_ROM_INT(NDARRAY_UINT16) },\n",
" { MP_ROM_QSTR(MP_QSTR_int16), MP_ROM_INT(NDARRAY_INT16) },\n",
" { MP_ROM_QSTR(MP_QSTR_float), MP_ROM_INT(NDARRAY_FLOAT) },\n",
"};\n",
"\n",
"STATIC MP_DEFINE_CONST_DICT (\n",
" mp_module_ulab_globals,\n",
" ulab_globals_table\n",
");\n",
"\n",
"const mp_obj_module_t ulab_user_cmodule = {\n",
" .base = { &mp_type_module },\n",
" .globals = (mp_obj_dict_t*)&mp_module_ulab_globals,\n",
"};\n",
"\n",
"MP_REGISTER_MODULE(MP_QSTR_ulab, ulab_user_cmodule, MODULE_ULAB_ENABLED);\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## makefile"
]
},
{
"cell_type": "code",
"execution_count": 640,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-16T18:44:45.794621Z",
"start_time": "2020-02-16T18:44:45.790112Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting ../../../ulab/code/micropython.mk\n"
]
}
],
"source": [
"%%writefile ../../../ulab/code/micropython.mk\n",
"\n",
"USERMODULES_DIR := $(USERMOD_DIR)\n",
"\n",
"# Add all C files to SRC_USERMOD.\n",
"SRC_USERMOD += $(USERMODULES_DIR)/ndarray.c\n",
"SRC_USERMOD += $(USERMODULES_DIR)/linalg.c\n",
"SRC_USERMOD += $(USERMODULES_DIR)/vectorise.c\n",
"SRC_USERMOD += $(USERMODULES_DIR)/poly.c\n",
"SRC_USERMOD += $(USERMODULES_DIR)/fft.c\n",
"SRC_USERMOD += $(USERMODULES_DIR)/numerical.c\n",
"SRC_USERMOD += $(USERMODULES_DIR)/filter.c\n",
"SRC_USERMOD += $(USERMODULES_DIR)/extras.c\n",
"SRC_USERMOD += $(USERMODULES_DIR)/ulab.c\n",
"\n",
"# We can add our module folder to include paths if needed\n",
"# This is not actually needed in this example.\n",
"CFLAGS_USERMOD += -I$(USERMODULES_DIR)\n",
"\n",
"CFLAGS_EXTRA = -DMODULE_ULAB_ENABLED=1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## make"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### unix port"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T05:53:58.238857Z",
"start_time": "2020-02-26T05:53:58.230416Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/v923z/sandbox/micropython/v1.11/micropython/ports/unix\n"
]
}
],
"source": [
"%cd ../../../micropython/ports/unix/"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T06:21:03.988272Z",
"start_time": "2020-02-26T06:21:03.736322Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Use make V=1 or set BUILD_VERBOSE in your environment to increase build verbosity.\n",
"rm -f micropython\n",
"rm -f micropython.map\n",
"rm -rf build-standard \n"
]
}
],
"source": [
"!make clean"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T06:21:59.627556Z",
"start_time": "2020-02-26T06:21:55.037569Z"
},
"scrolled": true
},
"outputs": [],
"source": [
"!make USER_C_MODULES=../../../ulab all"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### stm32 port"
]
},
{
"cell_type": "code",
"execution_count": 635,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-16T18:34:25.810862Z",
"start_time": "2020-02-16T18:34:25.805101Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/v923z/sandbox/micropython/v1.11/micropython/ports/stm32\n"
]
}
],
"source": [
"%cd ../../../micropython/ports/stm32/"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-16T18:35:25.758071Z",
"start_time": "2020-02-16T18:34:27.683963Z"
},
"scrolled": true
},
"outputs": [],
"source": [
"!rm -rf build-PYBV11\n",
"!make BOARD=PYBV11 USER_C_MODULES=../../../ulab all"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Change log"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {
"ExecuteTime": {
"end_time": "2020-02-26T06:23:10.596014Z",
"start_time": "2020-02-26T06:23:10.583737Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting ../../../ulab/docs/ulab-change-log.md\n"
]
}
],
"source": [
"%%writefile ../../../ulab/docs/ulab-change-log.md\n",
"\n",
"Sun, 16 Feb 2020\n",
"\n",
"version 0.33.2\n",
"\n",
" moved properties into ndarray_properties.h, implemented pointer arithmetic in fft.c to save some time\n",
"\n",
"Fri, 14 Feb 2020\n",
"\n",
"version 0.33.1\n",
"\n",
" added the __name__attribute to all sub-modules\n",
"\n",
"Thu, 13 Feb 2020\n",
"\n",
"version 0.33.0\n",
"\n",
" sub-modules are now proper sub-modules of ulab\n",
" \n",
"Tue, 11 Feb 2020\n",
"\n",
"version 0.32.0\n",
"\n",
" added itemsize, size and shape attributes to ndarrays, and removed rawsize\n",
"\n",
"Mon, 10 Feb 2020\n",
"\n",
"version 0.31.0\n",
"\n",
" removed asbytearray, and added buffer protocol to ndarrays, fixed bad error in filter.c\n",
" \n",
"Sun, 09 Feb 2020\n",
"\n",
"version 0.30.2\n",
"\n",
" fixed slice_length in ndarray.c\n",
"\n",
"Sat, 08 Feb 2020\n",
"\n",
"version 0.30.1\n",
"\n",
" fixed typecode error, added variable inspection, and replaced ternary operators in filter.c\n",
" \n",
"Fri, 07 Feb 2020\n",
"\n",
"version 0.30.0\n",
"\n",
" ulab functions can arbitrarily be excluded from the firmware via the ulab.h configuration file\n",
" \n",
"Thu, 06 Feb 2020\n",
"\n",
"version 0.27.0\n",
"\n",
" add convolve, the start of a 'filter' functionality group\n",
"\n",
"Wed, 29 Jan 2020\n",
"\n",
"version 0.26.7\n",
"\n",
" fixed indexing error in linalg.dot\n",
"\n",
"Mon, 20 Jan 2020\n",
"\n",
"version 0.26.6\n",
"\n",
" replaced MP_ROM_PTR(&mp_const_none_obj), so that module can be compiled for the nucleo board\n",
"\n",
"Tue, 7 Jan 2020\n",
"\n",
"version 0.26.5\n",
"\n",
" fixed glitch in numerical.c, numerical.h\n",
"\n",
"Mon, 6 Jan 2020\n",
"\n",
"version 0.26.4\n",
"\n",
" switched version constant to string\n",
"\n",
"Tue, 31 Dec 2019\n",
"\n",
"version 0.263\n",
"\n",
" changed declaration of ulab_ndarray_type to extern\n",
"\n",
"Fri, 29 Nov 2019\n",
"\n",
"version 0.262\n",
"\n",
" fixed error in macro in vectorise.h\n",
"\n",
"Thu, 28 Nov 2019\n",
"\n",
"version 0.261\n",
"\n",
" fixed bad indexing error in linalg.dot\n",
"\n",
"Tue, 6 Nov 2019\n",
"\n",
"version 0.26\n",
"\n",
" added in-place sorting (method of ndarray), and argsort\n",
" \n",
"Mon, 4 Nov 2019\n",
"\n",
"version 0.25\n",
"\n",
" added first implementation of sort, and fixed section on compiling the module in the manual\n",
"\n",
"Thu, 31 Oct 2019\n",
"\n",
"version 0.24\n",
"\n",
" added diff to numerical.c\n",
" \n",
"Tue, 29 Oct 2019\n",
"\n",
"version 0.23\n",
"\n",
" major revamp of subscription method\n",
"\n",
"Sat, 19 Oct 2019\n",
"\n",
"version 0.21\n",
"\n",
" fixed trivial bug in .rawsize()\n",
"\n",
"Sat, 19 Oct 2019\n",
"\n",
"version 0.22\n",
"\n",
" fixed small error in linalg_det, and implemented linalg_eig.\n",
"\n",
"\n",
"Thu, 17 Oct 2019\n",
"\n",
"version 0.21\n",
"\n",
" implemented uniform interface for fft, and spectrum, and added ifft.\n",
"\n",
"Wed, 16 Oct 2019\n",
"\n",
"version 0.20\n",
"\n",
" Added flip function to numerical.c, and moved the size function to linalg. In addition, \n",
" size is a function now, and not a method.\n",
"\n",
"Tue, 15 Oct 2019\n",
"\n",
"version 0.19\n",
"\n",
" fixed roll in numerical.c: it can now accept the axis=None keyword argument, added determinant to linalg.c\n",
"\n",
"Mon, 14 Oct 2019\n",
"\n",
"version 0.18\n",
"\n",
" fixed min/man function in numerical.c; it conforms to numpy behaviour\n",
"\n",
"Fri, 11 Oct 2019\n",
"\n",
"version 0.171\n",
"\n",
" found and fixed small bux in roll function\n",
"\n",
"Fri, 11 Oct 2019\n",
"\n",
"version 0.17\n",
"\n",
" universal function can now take arbitrary typecodes\n",
"\n",
"Fri, 11 Oct 2019\n",
"\n",
"version 0.161\n",
"\n",
" fixed bad error in iterator, and make_new_ndarray \n",
" \n",
"Thu, 10 Oct 2019\n",
"\n",
"varsion 0.16\n",
"\n",
" changed ndarray to array in ulab.c, so as to conform to numpy's notation\n",
" extended subscr method to include slices (partially works)\n",
" \n",
"Tue, 8 Oct 2019\n",
"\n",
"version 0.15\n",
"\n",
" added inv, neg, pos, and abs unary operators to ndarray.c\n",
" \n",
"Mon, 7 Oct 2019\n",
"\n",
"version 0.14\n",
"\n",
" made the internal binary_op function tighter, and added keyword arguments to linspace\n",
" \n",
"Sat, 4 Oct 2019\n",
"\n",
"version 0.13\n",
"\n",
" added the <, <=, >, >= binary operators to ndarray\n",
"\n",
"Fri, 4 Oct 2019\n",
"\n",
"version 0.12\n",
"\n",
" added .flatten to ndarray, ones, zeros, and eye to linalg\n",
"\n",
"Thu, 3 Oct 2019\n",
"\n",
"version 0.11\n",
" \n",
" binary operators are now based on macros"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {
"height": "calc(100% - 180px)",
"left": "10px",
"top": "150px",
"width": "382.797px"
},
"toc_section_display": true,
"toc_window_display": true
},
"toc-autonumbering": true,
"toc-showcode": false,
"toc-showmarkdowntxt": true,
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Reference in a new issue View git blame Copy permalink