ipython3.. code :: %pylab inline .. parsed-literal:: Populating the interactive namespace from numpy and matplotlib Introduction ============ 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. The main points of ulab are - 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.) - vectorised computations on micropython iterables and numerical arrays/matrices (universal functions) - basic linear algebra routines (matrix inversion, matrix reshaping, and transposition) - polynomial fits to numerical data - fast Fourier transforms 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. The code and its documentation can be found under https://github.com/v923z/micropython-ulab/. The MIT licence applies to all material. Environmental settings and magic commands ========================================= 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. 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. ipython3.. code :: %cd ../../micropython/ports/unix/ .. parsed-literal:: /home/v923z/sandbox/micropython/v1.11/micropython/ports/unix ipython3.. code :: from IPython.core.magic import Magics, magics_class, line_cell_magic from IPython.core.magic import cell_magic, register_cell_magic, register_line_magic from IPython.core.magic_arguments import argument, magic_arguments, parse_argstring import subprocess import os ipython3.. code :: def string_to_matrix(string): matrix = [] string = string.replace("array(\'d\', ", '').replace(')', '').replace('[', '').replace(']', '') for _str in string.split('\r\n'): if len(_str) > 0: matrix.append([float(n) for n in _str.split(',')]) return array(matrix) micropython magic command ------------------------- 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. ipython3.. code :: @magics_class class PyboardMagic(Magics): @cell_magic @magic_arguments() @argument('-skip') @argument('-unix') @argument('-file') @argument('-data') @argument('-time') @argument('-memory') def micropython(self, line='', cell=None): args = parse_argstring(self.micropython, line) if args.skip: # doesn't care about the cell's content print('skipped execution') return None # do not parse the rest if args.unix: # tests the code on the unix port. Note that this works on unix only with open('/dev/shm/micropython.py', 'w') as fout: fout.write(cell) proc = subprocess.Popen(["./micropython", "/dev/shm/micropython.py"], stdout=subprocess.PIPE, stderr=subprocess.PIPE) print(proc.stdout.read().decode("utf-8")) print(proc.stderr.read().decode("utf-8")) return None if args.file: # can be used to copy the cell content onto the pyboard's flash spaces = " " try: with open(args.file, 'w') as fout: fout.write(cell.replace('\t', spaces)) printf('written cell to {}'.format(args.file)) except: print('Failed to write to disc!') return None # do not parse the rest if args.data: # can be used to load data from the pyboard directly into kernel space message = pyb.exec(cell) if len(message) == 0: print('pyboard >>>') else: print(message.decode('utf-8')) # register new variable in user namespace self.shell.user_ns[args.data] = string_to_matrix(message.decode("utf-8")) if args.time: # measures the time of executions pyb.exec('import utime') message = pyb.exec('t = utime.ticks_us()\n' + cell + '\ndelta = utime.ticks_diff(utime.ticks_us(), t)' + "\nprint('execution time: {:d} us'.format(delta))") print(message.decode('utf-8')) if args.memory: # prints out memory information message = pyb.exec('from micropython import mem_info\nprint(mem_info())\n') print("memory before execution:\n========================\n", message.decode('utf-8')) message = pyb.exec(cell) print(">>> ", message.decode('utf-8')) message = pyb.exec('print(mem_info())') print("memory after execution:\n========================\n", message.decode('utf-8')) else: message = pyb.exec(cell) print(message.decode('utf-8')) ip = get_ipython() ip.register_magics(PyboardMagic) pyboard initialisation ~~~~~~~~~~~~~~~~~~~~~~ ipython3.. code :: import pyboard pyb = pyboard.Pyboard('/dev/ttyACM0') pyb.enter_raw_repl() pyboad detach ~~~~~~~~~~~~~ ipython3.. code :: pyb.exit_raw_repl() pyb.close() ipython3.. code :: import IPython js = """ (function () { var defaults = IPython.CodeCell.config_defaults || IPython.CodeCell.options_default; defaults.highlight_modes['magic_text/x-csrc'] = {'reg':[/^\\s*%%ccode/]}; })(); """ cjs = """ IPython.CodeCell.options_default.highlight_modes['magic_text/x-csrc'] = {'reg':[/^\\s*%%ccode/]}; """ IPython.core.display.display_javascript(cjs, raw=True) js = """ (function () { var defaults = IPython.CodeCell.config_defaults || IPython.CodeCell.options_default; defaults.highlight_modes['magic_text/x-csrc'] = {'reg':[/^\\s*%%makefile/]}; })(); """ IPython.core.display.display_javascript(js, raw=True) Code magic ---------- 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. ipython3.. code :: @magics_class class MyMagics(Magics): @cell_magic def ccode(self, line, cell): copyright = """/* * This file is part of the micropython-ulab project, * * https://github.com/v923z/micropython-ulab * * The MIT License (MIT) * * Copyright (c) 2019 Zoltán Vörös */ """ if line: with open('../../../ulab/code/'+line, 'w') as cout: cout.write(copyright) cout.write(cell) print('written %d bytes to %s'%(len(copyright) + len(cell), line)) return None ip = get_ipython() ip.register_magics(MyMagics) Notebook conversion =================== ipython3.. code :: %cd ../../../ulab/docs/ .. parsed-literal:: /home/v923z/sandbox/micropython/v1.11/ulab/docs ipython3.. code :: import nbformat as nb import nbformat.v4.nbbase as nb4 from nbconvert import RSTExporter def convert_notebook(nbfile, rstfile): (rst, resources) = rstexporter.from_filename(nbfile) with open(rstfile, 'w') as fout: fout.write(rst) rstexporter = RSTExporter() rstexporter.template_file = './templates/rst.tpl' convert_notebook('ulab.ipynb', './source/ulab.rst') .. parsed-literal:: /home/v923z/anaconda3/lib/python3.7/site-packages/nbconvert/filters/datatypefilter.py:41: UserWarning: Your element with mimetype(s) dict_keys(['application/javascript']) is not able to be represented. mimetypes=output.keys()) Compiling the module ==================== Detailed instructions on how to set up and compile a C module can be found in chapter 2 of https://micropython-usermod.readthedocs.io/en/latest/. First, on the command line, you should clone both the micropython, and the ``ulab`` repositories: .. code:: bash git clone https://github.com/micropython/micropython.git Then navigate to your micropython folder, and run .. code:: bash git clone https://github.com/v923z/micropython-ulab.git ulab Finally, in the ``mpconfigport.h`` header file of the port that you want to compile for, you have to define the variable ``MODULE_ULAB_ENABLED`` .. code:: make #define MODULE_ULAB_ENABLED (1) At this point, you should be able to run make in the port’s root folder: .. code:: bash make USER_C_MODULES=../../../ulab all (unix port) or .. code:: bash make BOARD=PYBV11 CROSS_COMPILE=/bin/arm-none-eabi- (pyboard). When compiling for the pyboard (or any other hardware platform), you might or might not have to set the cross-compiler’s path. If your installation of the cross-compiler is system-wide, you can drop the ``make`` argument ``CROSS_COMPILE``. The ndarray type ================ General comments ---------------- ``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. The type definition is as follows: .. code:: c typedef struct _ndarray_obj_t { mp_obj_base_t base; size_t m, n; mp_obj_array_t *data; size_t bytes; } ndarray_obj_t; **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.** Handling different types ------------------------ In order to make the code type-agnostic, we will resort to macros, where necessary. This will inevitably add to the firmware size, because, in effect, we unroll the code for each possible case. However, the source will be much more readable. Also note that by unrolling, we no longer need intermediate containers and we no longer need to dispatch type-conversion functions, which means that we should be able to gain in speed. Additional structure members in numpy ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Also note that, in addition, ``numpy`` defines the following members: - ``.ndim``: the number of dimensions of the array (in our case, it would be 1, or 2) - ``.size``: the number of elements in the array; it is the product of m, and n - ``.dtype``: the data type; in our case, it is basically stored in data->typecode - ``.itemsize``: the size of a single element in the array: this can be gotten by calling ``mp_binary_get_size('@', data->typecode, NULL)``. One should, perhaps, consider, whether these are necessary fields. E.g., if ``ndim`` were defined, then .. code:: c if((myarray->m == 1) || (myarray->n == 1)) { ... } could be replaced by .. code:: c if(myarray->ndim == 1) { ... } and .. code:: c if((myarray->m > 1) && (myarray->n > 1)) { ... } would be equivalent to .. code:: c if(myarray->ndim == 2) { ... } One could also save the extra function call ``mp_binary_get_size('@', data->typecode, NULL)``, if ``itemsize`` is available. Returning and accepting raw bytes ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ It might make sense to have a function that returns the raw content of the ``ndarray``. The rationale for this is that this would make direct use of calculation results a piece of cake. E.g., the DAC could be fed as .. code:: python length = 100 amp = 127 x = linspace(0, 2*pi, length) y = ndarray(128 + amp*sin(x), dtype=uint8) buf = y.bytearray() dac = DAC(1) dac.write_timed(buf, 400*length, mode=DAC.CIRCULAR) Likewise, having the option of writing raw data directly into the ``ndarray`` could simplify data analysis. E.g., ADC results could be processed as follows: .. code:: python length = 100 y = ndarray([0]*length, dtype=uint16) adc = pyb.ADC(pyb.Pin.board.X19) tim = pyb.Timer(6, freq=10) buf = y.bytearray() adc.read_timed(buf, tim) y.reshape((10, 10)) # or whatever Exposed functions and properties ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Most of the functions in ndarray.c are internal (i.e., not exposed to the python interpreter). Exception to this rule are the ``shape``, ``size``, and ``rawsize`` functions, and the ``.unary_op``, ``.binary_op``, and ``.iter_next`` class methods. Note that ``rawsize`` is is not standard in numpy, and is meant to return the total number of bytes used by the container. Since the RAM of a microcontroller is limited, I deemed this to be a reasonable addition for optimisation purposes, but could later be removed, if it turns out to be of no interest. As mentioned above, ``numpy`` defines a number of extra members in its ``ndarray``. It would be great, if we could return these members as properties of the ``ndarray`` instance. At the moment, ``shape`` is a function, as is ``rawsize``. Initialisation -------------- An ``ndarray`` can be initialised by passing an iterable (linear arrays), or an iterable of iterables (matrices) into the constructor. In addition, the constructor can take a keyword argument, ``dtype``, that will force type conversion. The default value is ``float``. ipython3.. code :: %%micropython -unix 1 from ulab import ndarray a = ndarray([1, 2, 3, 4]) print(a) a = ndarray([[1, 2, 3, 4], [2, 3, 4, 5]]) print('\n', a) a = ndarray([range(10), range(10)]) print('\n', a) .. parsed-literal:: ndarray([1.0, 2.0, 3.0, 4.0], dtype=float) ndarray([[1.0, 2.0, 3.0, 4.0], [2.0, 3.0, 4.0, 5.0]], dtype=float) ndarray([[0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0], [0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0]], dtype=float) Slicing and subscriptions ------------------------- Slicing and subscription work work only partially. This has to do with the fact that slices can be given only in the form ``myobject[1:10:2]``, and slice objects cannot be separated by commas, i.e., this won’t be interpreted: .. code:: python ndarray[1:10:2, 3:10:3] This is something that should be sorted out in the future. As a temporary solution, we could implement the ``getter`` and ``setter`` special methods that handle this issue. ipython3.. code :: %%micropython -unix 1 from ulab import ndarray # initialise a matrix a = ndarray([[1, 2, 3, 4], [6, 7, 8, 9]]) print('2D array: \n', a) # print out the second row print('second row of matrix: ', a[1]) #initialise an array a = ndarray([1, 2, 3, 4, 5, 6, 7, 8, 9]) print('\n1D array: ', a) # slice the array print('slize between 1, and 5: ', a[1:5]) .. parsed-literal:: 2D array: ndarray([[1.0, 2.0, 3.0, 4.0], [6.0, 7.0, 8.0, 9.0]], dtype=float) second row of matrix: ndarray([6.0, 7.0, 8.0, 9.0], dtype=float) 1D array: ndarray([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0], dtype=float) slize between 1, and 5: ndarray([2.0, 3.0, 4.0, 5.0], dtype=float) Iterators --------- ``ndarray`` objects can be iterated on, and just as in numpy, matrices are iterated along their first axis, and they return ``ndarray``\ s. ipython3.. code :: %%micropython -unix 1 from ulab import ndarray # initialise a matrix a = ndarray([[1, 2, 3, 4], [6, 7, 8, 9]]) print('2D array: \n', a) # print out the matrix' rows, one by one for i, _a in enumerate(a): print('\nrow %d: '%i, _a) .. parsed-literal:: 2D array: ndarray([[1.0, 2.0, 3.0, 4.0], [6.0, 7.0, 8.0, 9.0]], dtype=float) row 0: ndarray([1.0, 2.0, 3.0, 4.0], dtype=float) row 1: ndarray([6.0, 7.0, 8.0, 9.0], dtype=float) On the other hand, flat arrays return their elements: ipython3.. code :: %%micropython -unix 1 from ulab import ndarray, uint8 # initialise an array a = ndarray(range(10), dtype=uint8) print('1D array: ', a) # print out the array's elements, one by one for i, _a in enumerate(a): print('element %d: '%i, _a) .. parsed-literal:: 1D array: ndarray([0, 1, 2, ..., 7, 8, 9], dtype=uint8) element 0: 0 element 1: 1 element 2: 2 element 3: 3 element 4: 4 element 5: 5 element 6: 6 element 7: 7 element 8: 8 element 9: 9 Upcasting --------- The following section shows the upcasting rules of ``numpy``, and immediately after each case, the test for ``ulab``. uint8 ~~~~~ ipython3.. code :: a = array([100], dtype=uint8) b = array([101], dtype=uint8) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=uint8), array([255], dtype=uint8), array([116], dtype=uint8), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.uint8) b = ulab.ndarray([101], dtype=ulab.uint8) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201], dtype=uint8) ndarray([255], dtype=uint8) ndarray([116], dtype=uint8) ndarray([0.9900990128517151], dtype=float) ipython3.. code :: a = array([100], dtype=uint8) b = array([101], dtype=int8) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=int16), array([-1], dtype=int16), array([10100], dtype=int16), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.uint8) b = ulab.ndarray([101], dtype=ulab.int8) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201], dtype=int16) ndarray([-1], dtype=int16) ndarray([10100], dtype=int16) ndarray([0.9900990128517151], dtype=float) ipython3.. code :: a = array([100], dtype=uint8) b = array([101], dtype=uint16) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=uint16), array([65535], dtype=uint16), array([10100], dtype=uint16), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.uint8) b = ulab.ndarray([101], dtype=ulab.uint16) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201], dtype=uint16) ndarray([65535], dtype=uint16) ndarray([10100], dtype=uint16) ndarray([0.9900990128517151], dtype=float) ipython3.. code :: a = array([100], dtype=uint8) b = array([101], dtype=int16) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=int16), array([-1], dtype=int16), array([10100], dtype=int16), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.uint8) b = ulab.ndarray([101], dtype=ulab.int16) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201], dtype=int16) ndarray([-1], dtype=int16) ndarray([10100], dtype=int16) ndarray([0.9900990128517151], dtype=float) ipython3.. code :: a = array([100], dtype=uint8) b = array([101], dtype=float) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201.]), array([-1.]), array([10100.]), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.uint8) b = ulab.ndarray([101], dtype=ulab.float) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201.0], dtype=float) ndarray([-1.0], dtype=float) ndarray([10100.0], dtype=float) ndarray([0.9900990128517151], dtype=float) int8 ~~~~ ipython3.. code :: a = array([100], dtype=int8) b = array([101], dtype=uint8) a + b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=int16), array([-1], dtype=int16), array([10100], dtype=int16), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.int8) b = ulab.ndarray([101], dtype=ulab.uint8) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201], dtype=int16) ndarray([-1], dtype=int16) ndarray([10100], dtype=int16) ndarray([0.9900990128517151], dtype=float) ipython3.. code :: a = array([100, 101], dtype=int8) b = array([200, 101], dtype=int8) a+b, a-b, a*b, a/b .. parsed-literal:: (array([ 44, -54], dtype=int8), array([-100, 0], dtype=int8), array([ 32, -39], dtype=int8), array([-1.78571429, 1. ])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100, 101], dtype=ulab.int8) b = ulab.ndarray([200, 101], dtype=ulab.int8) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([44, -54], dtype=int8) ndarray([-100, 0], dtype=int8) ndarray([32, -39], dtype=int8) ndarray([-1.785714268684387, 1.0], dtype=float) ipython3.. code :: a = array([100], dtype=int8) b = array([200], dtype=uint16) a+b, a-b, a*b, a/b .. parsed-literal:: (array([300], dtype=int32), array([-100], dtype=int32), array([20000], dtype=int32), array([0.5])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.int8) b = ulab.ndarray([200], dtype=ulab.uint16) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([300], dtype=int16) ndarray([-100], dtype=int16) ndarray([20000], dtype=int16) ndarray([0], dtype=int16) ipython3.. code :: a = array([100], dtype=int8) b = array([200], dtype=int16) a+b, a-b, a*b, a/b .. parsed-literal:: (array([300], dtype=int16), array([-100], dtype=int16), array([20000], dtype=int16), array([0.5])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.int8) b = ulab.ndarray([200], dtype=ulab.int16) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([300], dtype=int16) ndarray([-100], dtype=int16) ndarray([20000], dtype=int16) ndarray([0.5], dtype=float) ipython3.. code :: a = array([100, 101], dtype=int8) b = array([200, 101], dtype=float) a+b, a-b, a*b, a/b .. parsed-literal:: (array([300., 202.]), array([-100., 0.]), array([20000., 10201.]), array([0.5, 1. ])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100, 101], dtype=ulab.int8) b = ulab.ndarray([200, 101], dtype=ulab.float) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([300.0, 202.0], dtype=float) ndarray([-100.0, 0.0], dtype=float) ndarray([20000.0, 10201.0], dtype=float) ndarray([0.5, 1.0], dtype=float) uint16 ~~~~~~ ipython3.. code :: a = array([100, 101], dtype=uint16) b = array([200, 101], dtype=uint8) a+b, a-b, a*b, a/b .. parsed-literal:: (array([300, 202], dtype=uint16), array([65436, 0], dtype=uint16), array([20000, 10201], dtype=uint16), array([0.5, 1. ])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100, 101], dtype=ulab.uint16) b = ulab.ndarray([200, 101], dtype=ulab.uint8) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([300, 202], dtype=uint16) ndarray([65436, 0], dtype=uint16) ndarray([20000, 10201], dtype=uint16) ndarray([0.5, 1.0], dtype=float) ipython3.. code :: a = array([100, 101], dtype=uint16) b = array([200, 101], dtype=int8) a+b, a-b, a*b, a/b .. parsed-literal:: (array([ 44, 202], dtype=int32), array([156, 0], dtype=int32), array([-5600, 10201], dtype=int32), array([-1.78571429, 1. ])) **This deviates from numpy’s behaviour, because we don’t have int32.** ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100, 101], dtype=ulab.uint16) b = ulab.ndarray([200, 101], dtype=ulab.int8) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([44, 202], dtype=uint16) ndarray([156, 0], dtype=uint16) ndarray([59936, 10201], dtype=uint16) ndarray([-1.785714268684387, 1.0], dtype=float) ipython3.. code :: a = array([100], dtype=uint16) b = array([101], dtype=uint16) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=uint16), array([65535], dtype=uint16), array([10100], dtype=uint16), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.uint16) b = ulab.ndarray([101], dtype=ulab.uint16) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201], dtype=uint16) ndarray([65535], dtype=uint16) ndarray([10100], dtype=uint16) ndarray([0.9900990128517151], dtype=float) ipython3.. code :: a = array([100], dtype=uint16) b = array([101], dtype=int16) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=int32), array([-1], dtype=int32), array([10100], dtype=int32), array([0.99009901])) **Again, in numpy, the result is an int32** ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.uint16) b = ulab.ndarray([101], dtype=ulab.int16) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201.0], dtype=float) ndarray([-1.0], dtype=float) ndarray([10100.0], dtype=float) ndarray([0.9900990128517151], dtype=float) ipython3.. code :: a = array([100], dtype=uint16) b = array([101], dtype=float) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201.]), array([-1.]), array([10100.]), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.uint16) b = ulab.ndarray([101], dtype=ulab.float) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201.0], dtype=float) ndarray([-1.0], dtype=float) ndarray([10100.0], dtype=float) ndarray([0.9900990128517151], dtype=float) int16 ~~~~~ ipython3.. code :: a = array([100], dtype=int16) b = array([101], dtype=uint8) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=int16), array([-1], dtype=int16), array([10100], dtype=int16), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.int16) b = ulab.ndarray([101], dtype=ulab.uint8) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201], dtype=int16) ndarray([-1], dtype=int16) ndarray([10100], dtype=int16) ndarray([0.9900990128517151], dtype=float) ipython3.. code :: a = array([100], dtype=int16) b = array([101], dtype=int8) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=int16), array([-1], dtype=int16), array([10100], dtype=int16), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.int16) b = ulab.ndarray([101], dtype=ulab.int8) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201], dtype=int16) ndarray([-1], dtype=int16) ndarray([10100], dtype=int16) ndarray([0.9900990128517151], dtype=float) ipython3.. code :: a = array([100], dtype=int16) b = array([101], dtype=uint16) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=int32), array([-1], dtype=int32), array([10100], dtype=int32), array([0.99009901])) **While the results are correct, here we have float instead of int32** ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.int16) b = ulab.ndarray([101], dtype=ulab.uint16) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201.0], dtype=float) ndarray([-1.0], dtype=float) ndarray([10100.0], dtype=float) ndarray([0.0], dtype=float) ipython3.. code :: a = array([100], dtype=int16) b = array([101], dtype=int16) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201], dtype=int16), array([-1], dtype=int16), array([10100], dtype=int16), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.int16) b = ulab.ndarray([101], dtype=ulab.int16) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201], dtype=int16) ndarray([-1], dtype=int16) ndarray([10100], dtype=int16) ndarray([0], dtype=int16) ipython3.. code :: a = array([100], dtype=int16) b = array([101], dtype=float) a+b, a-b, a*b, a/b .. parsed-literal:: (array([201.]), array([-1.]), array([10100.]), array([0.99009901])) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([100], dtype=ulab.int16) b = ulab.ndarray([101], dtype=ulab.float) print(a+b) print(a-b) print(a*b) print(a/b) .. parsed-literal:: ndarray([201.0], dtype=float) ndarray([-1.0], dtype=float) ndarray([10100.0], dtype=float) ndarray([0.9900990128517151], dtype=float) When in an operation the ``dtype`` of two arrays is different, the result’s ``dtype`` will be decided by the following upcasting rules: 1. Operations with two ``ndarray``\ s of the same ``dtype`` preserve their ``dtype``, even when the results overflow. 2. if either of the operands is a float, the results is also a float 3. - ``uint8`` + ``int8`` => ``int16``, - ``uint8`` + ``int16`` => ``int16`` - ``uint8`` + ``uint16`` => ``uint16`` - ``int8`` + ``int16`` => ``int16`` - ``int8`` + ``uint16`` => ``uint16`` (in numpy, it is a ``int32``) - ``uint16`` + ``int16`` => ``float`` (in numpy, it is a ``int32``) 4. When the right hand side of a binary operator is a micropython variable, ``mp_obj_int``, or ``mp_obj_float``, then the result will be promoted to ``dtype`` ``float``. This is necessary, because a micropython integer can be 31 bites wide. ``numpy`` is also inconsistent in how it represents ``dtype``: as an argument, it is denoted by the constants ``int8``, ``uint8``, etc., while a string will be returned, if the user asks for the type of an array. The upcasting rules are stipulated in a single C function, ``upcasting()``. upcasting rules with scalars ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When a ipython3.. code :: a = array([1, 2, 3], dtype=int8) b = a * 555 a *= -555 b.dtype, b, a.dtype, a .. parsed-literal:: (dtype('int16'), array([ 555, 1110, 1665], dtype=int16), dtype('int8'), array([-43, -86, 127], dtype=int8)) Binary operations ----------------- In the standard binary operations, the operands are either two ``ndarray``\ s or an ``ndarray``, and a number. From the C standpoint, these operations are probably the most difficult: the problem is that two operands, each with 5 possible C types are added, multiplied, subtracted, or divided, hence making the number of possible combinations large. In order to mitigate the situation, we make use of macros: this would make most of the code type-agnostic. Also, when an operation involves a scalar, and an ``ndarray``, we will turn the scalar into an ``ndarray`` of length 1. In this way, we can reduce the code size of the binary handler by almost a factor of two. ipython3.. code :: %%micropython -unix 1 from ulab import ndarray, float a = ndarray([1, 2, 3], dtype=float) print(a + a) print(a * 5.0) print(a / 2) print(a - 10) .. parsed-literal:: ndarray([2.0, 4.0, 6.0], dtype=float) ndarray([5.0, 10.0, 15.0], dtype=float) ndarray([0.5, 1.0, 1.5], dtype=float) ndarray([-9.0, -8.0, -7.0], dtype=float) in-place operators ~~~~~~~~~~~~~~~~~~ In-place operators preserve the type of the array’s type. Here are a couple of caveats: 1. overflow obviously occurs 2. float can be added only to float type 3. true divide fails, except when the array is of type float ipython3.. code :: a = array([1, 2, 3, 40], dtype=uint8) a += 220 a .. parsed-literal:: array([221, 222, 223, 4], dtype=uint8) ipython3.. code :: a = array([1, 2, 3, 40], dtype=int8) a += 220 a .. parsed-literal:: array([-35, -34, -33, 4], dtype=int8) ipython3.. code :: a = array([1, 2, 3, 40], dtype=uint16) a += 220 a .. parsed-literal:: array([221, 222, 223, 260], dtype=uint16) ipython3.. code :: a = array([1, 2, 3, 40], dtype=int16) a += 220 a .. parsed-literal:: array([221, 222, 223, 260], dtype=int16) ipython3.. code :: a = array([1, 2, 3, 40], dtype=float) a += 220 a .. parsed-literal:: array([221., 222., 223., 260.]) ipython3.. code :: a = array([1, 2, 3, 40], dtype=uint8) a += 220.0 a :: --------------------------------------------------------------------------- TypeError Traceback (most recent call last) in 1 a = array([1, 2, 3, 40], dtype=uint8) ----> 2 a += 220.0 3 a TypeError: Cannot cast ufunc add output from dtype('float64') to dtype('uint8') with casting rule 'same_kind' ipython3.. code :: a = array([1, 2, 3, 40], dtype=uint8) a /= 22 a :: --------------------------------------------------------------------------- TypeError Traceback (most recent call last) in 1 a = array([1, 2, 3, 40], dtype=uint8) ----> 2 a /= 22 3 a TypeError: No loop matching the specified signature and casting was found for ufunc true_divide ipython3.. code :: a = array([1, 2, 3, 40], dtype=float) a += 220 a .. parsed-literal:: array([221., 222., 223., 260.]) ipython3.. code :: a = array([1, 2, 3, 4], dtype=int8) b = array([5, 10, 15, 20], dtype=float) a /= b a :: --------------------------------------------------------------------------- TypeError Traceback (most recent call last) in 1 a = array([1, 2, 3, 4], dtype=int8) 2 b = array([5, 10, 15, 20], dtype=float) ----> 3 a /= b 4 a TypeError: ufunc 'true_divide' output (typecode 'd') could not be coerced to provided output parameter (typecode 'b') according to the casting rule ''same_kind'' ipython3.. code :: a = array([1, 2, 3, 4], dtype=int8) b = array([5, 10, 15, 20], dtype=int8) a /= b a :: --------------------------------------------------------------------------- TypeError Traceback (most recent call last) in 1 a = array([1, 2, 3, 4], dtype=int8) 2 b = array([5, 10, 15, 20], dtype=int8) ----> 3 a /= b 4 a TypeError: No loop matching the specified signature and casting was found for ufunc true_divide ipython3.. code :: a = array([1, 2, 3, 4], dtype=int8) b = array([5, 10, 15, 100], dtype=int16) a *= b a .. parsed-literal:: array([ 5, 20, 45, -112], dtype=int8) ipython3.. code :: a = array([1, 2, 3, 4], dtype=int8) a **= 2 a .. parsed-literal:: array([ 1, 4, 9, 16], dtype=int8) Comparison operators ~~~~~~~~~~~~~~~~~~~~ These return list(s) of Booleans. ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([1, 2, 3, 4, 5, 6, 7, 8]) a.reshape((1, 8)) print(a < 4) print(a <= 4) print(a > 4) print(a >= 4) .. parsed-literal:: [True, True, True, False, False, False, False, False] [True, True, True, True, False, False, False, False] [False, False, False, False, True, True, True, True] [False, False, False, True, True, True, True, True] ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([[1, 2, 3, 4], [5, 6, 7, 8]]) print(a < 4) print(a <= 4) print(a > 4) print(a >= 4) .. parsed-literal:: [[True, True, True, False], [False, False, False, False]] [[True, True, True, True], [False, False, False, False]] [[False, False, False, False], [True, True, True, True]] [[False, False, False, True], [True, True, True, True]] ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([1, 2, 3, 4, 4, 6, 7, 8]) a.transpose() b = ulab.ndarray([8, 7, 6, 5, 4, 3, 2, 1]) b.transpose() print(a < b) print(a <= b) print(a > b) print(a >= b) .. parsed-literal:: [True, True, True, True, False, False, False, False] [True, True, True, True, True, False, False, False] [False, False, False, False, False, True, True, True] [False, False, False, False, True, True, True, True] ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([[1, 2, 3, 4], [5, 6, 7, 8]]) b = ulab.ndarray([[8, 7, 1, 1], [4, 3, 2, 1]], dtype=ulab.int8) print(a < b) print(a <= b) print(a > b) print(a >= b) .. parsed-literal:: [[True, True, False, False], [False, False, False, False]] [[True, True, False, False], [False, False, False, False]] [[False, False, True, True], [True, True, True, True]] [[False, False, True, True], [True, True, True, True]] Simple running weighted average ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ With the subscription tools, a weighted running average can very easily be implemented as follows: ipython3.. code :: %%micropython -unix 1 from ulab import ndarray, mean, roll # These are the weights; the last entry is the most dominant weight = ndarray([1, 2, 3, 4, 5]) # initial array of samples samples = ndarray([0]*5) for i in range(5): # a new datum is inserted on the right hand side. This simply overwrites whatever was in the last slot samples[-1] = 2 print(mean(samples*weight)) # the data are shifted by one position to the left roll(samples, 1) .. parsed-literal:: 2.0 3.6 4.8 5.6 6.0 Unary operators --------------- At the moment, only ``len`` is implemented, which returns the number of elements for one-dimensional arrays, and the length of the first axis for matrices. One should consider other possibilities. ipython3.. code :: %%micropython -unix 1 from ulab import ndarray # initialise an array a = ndarray(range(10)) print('1D array: ', a) # print out the array's length print('length of array: ', len(a)) .. parsed-literal:: 1D array: ndarray([0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0], dtype=float) length of array: 10 ipython3.. code :: %%micropython -unix 1 from ulab import ndarray # initialise a matrix a = ndarray([range(10), range(10), range(10)]) print('2D array: \n', a) # print out the array's elements, one by one print('length of array: ', len(a)) .. parsed-literal:: 2D array: ndarray([[0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0], [0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0], [0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0]], dtype=float) length of array: 3 Class methods: shape, size, rawsize, flatten -------------------------------------------- ipython3.. code :: %%micropython -unix 1 from ulab import ndarray # initialise an array a = ndarray(range(10)) print('1D array: ', a) # print out the shape print('shape: ', a.shape()) #print out the size print('size 0: ', a.size(0), '\nsize 1: ', a.size(1), '\nsize 2: ', a.size(2)) #print out the raw size print('raw size: ', a.rawsize()) # initialise a matrix a = ndarray([range(10), range(10), range(10)]) print('\n2D array: \n', a) # print out the shape print('shape: ', a.shape()) #print out the size print('size 0: ', a.size(0), '\nsize 1: ', a.size(1), '\nsize 2: ', a.size(2)) #print out the raw size print('raw size: ', a.rawsize()) #flattened array a = ndarray([range(3), range(3), range(3)]) print('\n2D array: \n', a) print('flattened array: (C)', a.flatten(order='C')) print('flattened array: (F)', a.flatten(order='F')) .. parsed-literal:: 1D array: ndarray([0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0], dtype=float) shape: (10, 1) size 0: 10 size 1: 10 size 2: 1 raw size: (10, 1, 40, 10, 4) 2D array: ndarray([[0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0], [0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0], [0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0]], dtype=float) shape: (3, 10) size 0: 30 size 1: 3 size 2: 10 raw size: (3, 10, 120, 30, 4) 2D array: ndarray([[0.0, 1.0, 2.0], [0.0, 1.0, 2.0], [0.0, 1.0, 2.0]], dtype=float) flattened array: (C) ndarray([0.0, 1.0, 2.0, 0.0, 1.0, 2.0, 0.0, 1.0, 2.0], dtype=float) flattened array: (F) ndarray([0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0], dtype=float) ndarray.h --------- https://github.com/v923z/micropython-ulab/tree/master/code/ndarray.h .. code:: cpp #ifndef _NDARRAY_ #define _NDARRAY_ #include "py/objarray.h" #include "py/binary.h" #include "py/objstr.h" #include "py/objlist.h" #define PRINT_MAX 10 const mp_obj_type_t ulab_ndarray_type; #define CREATE_SINGLE_ITEM(outarray, type, typecode, value) do {\ ndarray_obj_t *tmp = create_new_ndarray(1, 1, (typecode));\ type *tmparr = (type *)tmp->array->items;\ tmparr[0] = (type)(value);\ (outarray) = MP_OBJ_FROM_PTR(tmp);\ } while(0) /* mp_obj_t row = mp_obj_new_list(n, NULL); mp_obj_list_t *row_ptr = MP_OBJ_TO_PTR(row); should work outside the loop, but it doesn't. Go figure! */ #define RUN_BINARY_LOOP(typecode, type_out, type_left, type_right, ol, or, op) do {\ type_left *left = (type_left *)(ol)->array->items;\ type_right *right = (type_right *)(or)->array->items;\ uint8_t inc = 0;\ if((or)->array->len > 1) inc = 1;\ if(((op) == MP_BINARY_OP_ADD) || ((op) == MP_BINARY_OP_SUBTRACT) || ((op) == MP_BINARY_OP_MULTIPLY)) {\ ndarray_obj_t *out = create_new_ndarray(ol->m, ol->n, typecode);\ type_out *(odata) = (type_out *)out->array->items;\ if((op) == MP_BINARY_OP_ADD) { for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = left[i] + right[j];}\ if((op) == MP_BINARY_OP_SUBTRACT) { for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = left[i] - right[j];}\ if((op) == MP_BINARY_OP_MULTIPLY) { for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = left[i] * right[j];}\ return MP_OBJ_FROM_PTR(out);\ } else if((op) == MP_BINARY_OP_TRUE_DIVIDE) {\ ndarray_obj_t *out = create_new_ndarray(ol->m, ol->n, NDARRAY_FLOAT);\ float *odata = (float *)out->array->items;\ for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = (float)left[i]/(float)right[j];\ return MP_OBJ_FROM_PTR(out);\ } else if(((op) == MP_BINARY_OP_LESS) || ((op) == MP_BINARY_OP_LESS_EQUAL) || \ ((op) == MP_BINARY_OP_MORE) || ((op) == MP_BINARY_OP_MORE_EQUAL)) {\ mp_obj_t out_list = mp_obj_new_list(0, NULL);\ size_t m = (ol)->m, n = (ol)->n;\ for(size_t i=0, r=0; i < m; i++, r+=inc) {\ mp_obj_t row = mp_obj_new_list(n, NULL);\ mp_obj_list_t *row_ptr = MP_OBJ_TO_PTR(row);\ for(size_t j=0, s=0; j < n; j++, s+=inc) {\ row_ptr->items[j] = mp_const_false;\ if((op) == MP_BINARY_OP_LESS) {\ if(left[i*n+j] < right[r*n+s]) row_ptr->items[j] = mp_const_true;\ } else if((op) == MP_BINARY_OP_LESS_EQUAL) {\ if(left[i*n+j] <= right[r*n+s]) row_ptr->items[j] = mp_const_true;\ } else if((op) == MP_BINARY_OP_MORE) {\ if(left[i*n+j] > right[r*n+s]) row_ptr->items[j] = mp_const_true;\ } else if((op) == MP_BINARY_OP_MORE_EQUAL) {\ if(left[i*n+j] >= right[r*n+s]) row_ptr->items[j] = mp_const_true;\ }\ }\ if(m == 1) return row;\ mp_obj_list_append(out_list, row);\ }\ return out_list;\ }\ } while(0) enum NDARRAY_TYPE { NDARRAY_UINT8 = 'B', NDARRAY_INT8 = 'b', NDARRAY_UINT16 = 'H', NDARRAY_INT16 = 'h', NDARRAY_FLOAT = 'f', }; typedef struct _ndarray_obj_t { mp_obj_base_t base; size_t m, n; size_t len; mp_obj_array_t *array; size_t bytes; } ndarray_obj_t; mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t , size_t , mp_obj_iter_buf_t *); float ndarray_get_float_value(void *, uint8_t , size_t ); void fill_array_iterable(float *, mp_obj_t ); void ndarray_print_row(const mp_print_t *, mp_obj_array_t *, size_t , size_t ); void ndarray_print(const mp_print_t *, mp_obj_t , mp_print_kind_t ); void ndarray_assign_elements(mp_obj_array_t *, mp_obj_t , uint8_t , size_t *); ndarray_obj_t *create_new_ndarray(size_t , size_t , uint8_t ); mp_obj_t ndarray_copy(mp_obj_t ); mp_obj_t ndarray_make_new(const mp_obj_type_t *, size_t , size_t , const mp_obj_t *); mp_obj_t ndarray_subscr(mp_obj_t , mp_obj_t , mp_obj_t ); mp_obj_t ndarray_getiter(mp_obj_t , mp_obj_iter_buf_t *); mp_obj_t ndarray_binary_op(mp_binary_op_t , mp_obj_t , mp_obj_t ); mp_obj_t ndarray_unary_op(mp_unary_op_t , mp_obj_t ); mp_obj_t ndarray_shape(mp_obj_t ); mp_obj_t ndarray_size(mp_obj_t , mp_obj_t ); mp_obj_t ndarray_rawsize(mp_obj_t ); mp_obj_t ndarray_flatten(size_t , const mp_obj_t *, mp_map_t *); #endif ndarray.c --------- https://github.com/v923z/micropython-ulab/tree/master/code/ndarray.c .. code:: cpp #include #include #include #include #include "py/runtime.h" #include "py/binary.h" #include "py/obj.h" #include "py/objtuple.h" #include "ndarray.h" // This function is copied verbatim from objarray.c STATIC mp_obj_array_t *array_new(char typecode, size_t n) { int typecode_size = mp_binary_get_size('@', typecode, NULL); mp_obj_array_t *o = m_new_obj(mp_obj_array_t); // this step could probably be skipped: we are never going to store a bytearray per se #if MICROPY_PY_BUILTINS_BYTEARRAY && MICROPY_PY_ARRAY o->base.type = (typecode == BYTEARRAY_TYPECODE) ? &mp_type_bytearray : &mp_type_array; #elif MICROPY_PY_BUILTINS_BYTEARRAY o->base.type = &mp_type_bytearray; #else o->base.type = &mp_type_array; #endif o->typecode = typecode; o->free = 0; o->len = n; o->items = m_new(byte, typecode_size * o->len); return o; } float ndarray_get_float_value(void *data, uint8_t typecode, size_t index) { if(typecode == NDARRAY_UINT8) { return (float)((uint8_t *)data)[index]; } else if(typecode == NDARRAY_INT8) { return (float)((int8_t *)data)[index]; } else if(typecode == NDARRAY_UINT16) { return (float)((uint16_t *)data)[index]; } else if(typecode == NDARRAY_INT16) { return (float)((int16_t *)data)[index]; } else { return (float)((float_t *)data)[index]; } } void fill_array_iterable(float *array, mp_obj_t iterable) { mp_obj_iter_buf_t x_buf; mp_obj_t x_item, x_iterable = mp_getiter(iterable, &x_buf); size_t i=0; while ((x_item = mp_iternext(x_iterable)) != MP_OBJ_STOP_ITERATION) { array[i] = (float)mp_obj_get_float(x_item); i++; } } void ndarray_print_row(const mp_print_t *print, mp_obj_array_t *data, size_t n0, size_t n) { mp_print_str(print, "["); size_t i; if(n < PRINT_MAX) { // if the array is short, print everything mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0), PRINT_REPR); for(i=1; itypecode, data->items, n0+i), PRINT_REPR); } } else { mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0), PRINT_REPR); for(i=1; i<3; i++) { mp_print_str(print, ", "); mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+i), PRINT_REPR); } mp_printf(print, ", ..., "); mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+n-3), PRINT_REPR); for(size_t i=1; i<3; i++) { mp_print_str(print, ", "); mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+n-3+i), PRINT_REPR); } } mp_print_str(print, "]"); } void ndarray_print(const mp_print_t *print, mp_obj_t self_in, mp_print_kind_t kind) { (void)kind; ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in); mp_print_str(print, "ndarray("); if((self->m == 1) || (self->n == 1)) { ndarray_print_row(print, self->array, 0, self->array->len); } else { // TODO: add vertical ellipses for the case, when self->m > PRINT_MAX mp_print_str(print, "["); ndarray_print_row(print, self->array, 0, self->n); for(size_t i=1; i < self->m; i++) { mp_print_str(print, ",\n\t "); ndarray_print_row(print, self->array, i*self->n, self->n); } mp_print_str(print, "]"); } // TODO: print typecode if(self->array->typecode == NDARRAY_UINT8) { printf(", dtype=uint8)"); } else if(self->array->typecode == NDARRAY_INT8) { printf(", dtype=int8)"); } else if(self->array->typecode == NDARRAY_UINT16) { printf(", dtype=uint16)"); } else if(self->array->typecode == NDARRAY_INT16) { printf(", dtype=int16)"); } else if(self->array->typecode == NDARRAY_FLOAT) { printf(", dtype=float)"); } } void ndarray_assign_elements(mp_obj_array_t *data, mp_obj_t iterable, uint8_t typecode, size_t *idx) { // assigns a single row in the matrix mp_obj_t item; while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) { mp_binary_set_val_array(typecode, data->items, (*idx)++, item); } } ndarray_obj_t *create_new_ndarray(size_t m, size_t n, uint8_t typecode) { // Creates the base ndarray with shape (m, n), and initialises the values to straight 0s ndarray_obj_t *ndarray = m_new_obj(ndarray_obj_t); ndarray->base.type = &ulab_ndarray_type; ndarray->m = m; ndarray->n = n; mp_obj_array_t *array = array_new(typecode, m*n); ndarray->bytes = m * n * mp_binary_get_size('@', typecode, NULL); // this should set all elements to 0, irrespective of the of the typecode (all bits are zero) // we could, perhaps, leave this step out, and initialise the array only, when needed memset(array->items, 0, ndarray->bytes); ndarray->array = array; return ndarray; } mp_obj_t ndarray_copy(mp_obj_t self_in) { // returns a verbatim (shape and typecode) copy of self_in ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in); ndarray_obj_t *out = create_new_ndarray(self->m, self->n, self->array->typecode); memcpy(out->array->items, self->array->items, self->bytes); return MP_OBJ_FROM_PTR(out); } STATIC uint8_t ndarray_init_helper(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { static const mp_arg_t allowed_args[] = { { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} }, { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT } }, }; mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); uint8_t dtype = args[1].u_int; return dtype; } mp_obj_t ndarray_make_new(const mp_obj_type_t *type, size_t n_args, size_t n_kw, const mp_obj_t *args) { mp_arg_check_num(n_args, n_kw, 1, 2, true); mp_map_t kw_args; mp_map_init_fixed_table(&kw_args, n_kw, args + n_args); uint8_t dtype = ndarray_init_helper(n_args, args, &kw_args); size_t len1, len2=0, i=0; mp_obj_t len_in = mp_obj_len_maybe(args[0]); if (len_in == MP_OBJ_NULL) { mp_raise_ValueError("first argument must be an iterable"); } else { len1 = MP_OBJ_SMALL_INT_VALUE(len_in); } // We have to figure out, whether the first element of the iterable is an iterable itself // Perhaps, there is a more elegant way of handling this mp_obj_iter_buf_t iter_buf1; mp_obj_t item1, iterable1 = mp_getiter(args[0], &iter_buf1); while ((item1 = mp_iternext(iterable1)) != MP_OBJ_STOP_ITERATION) { len_in = mp_obj_len_maybe(item1); if(len_in != MP_OBJ_NULL) { // indeed, this seems to be an iterable // Next, we have to check, whether all elements in the outer loop have the same length if(i > 0) { if(len2 != MP_OBJ_SMALL_INT_VALUE(len_in)) { mp_raise_ValueError("iterables are not of the same length"); } } len2 = MP_OBJ_SMALL_INT_VALUE(len_in); i++; } } // By this time, it should be established, what the shape is, so we can now create the array ndarray_obj_t *self = create_new_ndarray(len1, (len2 == 0) ? 1 : len2, dtype); iterable1 = mp_getiter(args[0], &iter_buf1); i = 0; if(len2 == 0) { // the first argument is a single iterable ndarray_assign_elements(self->array, iterable1, dtype, &i); } else { mp_obj_iter_buf_t iter_buf2; mp_obj_t iterable2; while ((item1 = mp_iternext(iterable1)) != MP_OBJ_STOP_ITERATION) { iterable2 = mp_getiter(item1, &iter_buf2); ndarray_assign_elements(self->array, iterable2, dtype, &i); } } return MP_OBJ_FROM_PTR(self); } mp_obj_t ndarray_subscr(mp_obj_t self_in, mp_obj_t index, mp_obj_t value) { // NOTE: this will work only on the flattened array! ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in); if (value == MP_OBJ_SENTINEL) { // simply return the values at index, no assignment if (MP_OBJ_IS_TYPE(index, &mp_type_slice)) { mp_bound_slice_t slice; mp_seq_get_fast_slice_indexes(self->array->len, index, &slice); // TODO: this won't work with index reversion!!! size_t len = (slice.stop - slice.start) / slice.step; ndarray_obj_t *out = create_new_ndarray(1, len, self->array->typecode); int _sizeof = mp_binary_get_size('@', self->array->typecode, NULL); uint8_t *indata = (uint8_t *)self->array->items; uint8_t *outdata = (uint8_t *)out->array->items; for(size_t i=0; i < len; i++) { memcpy(outdata+(i*_sizeof), indata+(slice.start+i*slice.step)*_sizeof, _sizeof); } return MP_OBJ_FROM_PTR(out); } // we have a single index, return either a single number (arrays), or an array (matrices) int16_t idx = mp_obj_get_int(index); if(idx < 0) { idx = self->m > 1 ? self->m + idx : self->n + idx; } if(self->m > 1) { // we do have a matrix if(idx >= self->m) { mp_raise_ValueError("index is out of range"); } if(self->n == 1) { // the matrix is actually a column vector return mp_binary_get_val_array(self->array->typecode, self->array->items, idx); } // return an array ndarray_obj_t *out = create_new_ndarray(1, self->n, self->array->typecode); int _sizeof = mp_binary_get_size('@', self->array->typecode, NULL); uint8_t *indata = (uint8_t *)self->array->items; uint8_t *outdata = (uint8_t *)out->array->items; memcpy(outdata, &indata[idx*self->n*_sizeof], self->n*_sizeof); return MP_OBJ_FROM_PTR(out); } // since self->m == 1, we have a flat array, hence, we've got to return a single number if(idx >= self->n) { mp_raise_ValueError("index is out of range"); } return mp_binary_get_val_array(self->array->typecode, self->array->items, idx); } else { int16_t idx = mp_obj_get_int(index); if((self->m == 1) || (self->n == 1)) { if(idx < 0) { idx = self->m > 1 ? self->m + idx : self->n + idx; } if((idx > self->m) && (idx > self->n)) { mp_raise_ValueError("index is out of range"); } mp_binary_set_val_array(self->array->typecode, self->array->items, idx, value); } else { // do not deal with assignment, bail out, if the array is two-dimensional mp_raise_NotImplementedError("subcript assignment is not implemented for 2D arrays"); } } return mp_const_none; } // itarray iterator mp_obj_t ndarray_getiter(mp_obj_t o_in, mp_obj_iter_buf_t *iter_buf) { return mp_obj_new_ndarray_iterator(o_in, 0, iter_buf); } typedef struct _mp_obj_ndarray_it_t { mp_obj_base_t base; mp_fun_1_t iternext; mp_obj_t ndarray; size_t cur; } mp_obj_ndarray_it_t; mp_obj_t ndarray_iternext(mp_obj_t self_in) { mp_obj_ndarray_it_t *self = MP_OBJ_TO_PTR(self_in); ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(self->ndarray); // TODO: in numpy, ndarrays are iterated with respect to the first axis. size_t iter_end = 0; if((ndarray->m == 1) || (ndarray->n ==1)) { iter_end = ndarray->array->len; } else { iter_end = ndarray->m; } if(self->cur < iter_end) { if(ndarray->m == ndarray->array->len) { // we have a linear array // read the current value mp_obj_t value; value = mp_binary_get_val_array(ndarray->array->typecode, ndarray->array->items, self->cur); self->cur++; return value; } else { // we have a matrix, return the number of rows ndarray_obj_t *value = create_new_ndarray(1, ndarray->n, ndarray->array->typecode); // copy the memory content here uint8_t *tmp = (uint8_t *)ndarray->array->items; size_t strip_size = ndarray->n * mp_binary_get_size('@', ndarray->array->typecode, NULL); memcpy(value->array->items, &tmp[self->cur*strip_size], strip_size); self->cur++; return value; } } else { return MP_OBJ_STOP_ITERATION; } } mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t ndarray, size_t cur, mp_obj_iter_buf_t *iter_buf) { assert(sizeof(mp_obj_ndarray_it_t) <= sizeof(mp_obj_iter_buf_t)); mp_obj_ndarray_it_t *o = (mp_obj_ndarray_it_t*)iter_buf; o->base.type = &mp_type_polymorph_iter; o->iternext = ndarray_iternext; o->ndarray = ndarray; o->cur = cur; return MP_OBJ_FROM_PTR(o); } mp_obj_t ndarray_shape(mp_obj_t self_in) { ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in); mp_obj_t tuple[2] = { mp_obj_new_int(self->m), mp_obj_new_int(self->n) }; return mp_obj_new_tuple(2, tuple); } mp_obj_t ndarray_size(mp_obj_t self_in, mp_obj_t axis) { ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in); uint8_t ax = mp_obj_get_int(axis); if(ax == 0) { return mp_obj_new_int(self->array->len); } else if(ax == 1) { return mp_obj_new_int(self->m); } else if(ax == 2) { return mp_obj_new_int(self->n); } else { return mp_const_none; } } mp_obj_t ndarray_rawsize(mp_obj_t self_in) { // returns a 5-tuple with the // // 1. number of rows // 2. number of columns // 3. length of the storage (should be equal to the product of 1. and 2.) // 4. length of the data storage in bytes // 5. datum size in bytes ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in); mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(mp_obj_new_tuple(5, NULL)); tuple->items[0] = MP_OBJ_NEW_SMALL_INT(self->m); tuple->items[1] = MP_OBJ_NEW_SMALL_INT(self->n); tuple->items[2] = MP_OBJ_NEW_SMALL_INT(self->bytes); tuple->items[3] = MP_OBJ_NEW_SMALL_INT(self->array->len); tuple->items[4] = MP_OBJ_NEW_SMALL_INT(mp_binary_get_size('@', self->array->typecode, NULL)); return tuple; } mp_obj_t ndarray_flatten(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { static const mp_arg_t allowed_args[] = { { MP_QSTR_order, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_QSTR(MP_QSTR_C)} }, }; mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; mp_arg_parse_all(n_args - 1, pos_args + 1, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); mp_obj_t self_copy = ndarray_copy(pos_args[0]); ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(self_copy); GET_STR_DATA_LEN(args[0].u_obj, order, len); if((len != 1) || ((memcmp(order, "C", 1) != 0) && (memcmp(order, "F", 1) != 0))) { mp_raise_ValueError("flattening order must be either 'C', or 'F'"); } // if order == 'C', we simply have to set m, and n, there is nothing else to do if(memcmp(order, "F", 1) == 0) { ndarray_obj_t *self = MP_OBJ_TO_PTR(pos_args[0]); uint8_t _sizeof = mp_binary_get_size('@', self->array->typecode, NULL); // get the data of self_in: we won't need a temporary buffer for the transposition uint8_t *self_array = (uint8_t *)self->array->items; uint8_t *array = (uint8_t *)ndarray->array->items; size_t i=0; for(size_t n=0; n < self->n; n++) { for(size_t m=0; m < self->m; m++) { memcpy(array+_sizeof*i, self_array+_sizeof*(m*self->n + n), _sizeof); i++; } } } ndarray->n = ndarray->array->len; ndarray->m = 1; return self_copy; } // Binary operations mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) { // One of the operands is a scalar // TODO: conform to numpy with the upcasting // TODO: implement in-place operators mp_obj_t RHS = MP_OBJ_NULL; bool rhs_is_scalar = true; if(mp_obj_is_int(rhs)) { size_t ivalue = mp_obj_get_int(rhs); if((ivalue > 0) && (ivalue < 256)) { CREATE_SINGLE_ITEM(RHS, uint8_t, NDARRAY_UINT8, ivalue); } else if((ivalue > 0) && (ivalue < 65535)) { CREATE_SINGLE_ITEM(RHS, uint16_t, NDARRAY_UINT16, ivalue); } else if((ivalue < 0) && (ivalue > -128)) { CREATE_SINGLE_ITEM(RHS, int8_t, NDARRAY_INT8, ivalue); } else if((ivalue < 0) && (ivalue > -32767)) { CREATE_SINGLE_ITEM(RHS, int16_t, NDARRAY_INT16, ivalue); } } else if(mp_obj_is_float(rhs)) { float fvalue = mp_obj_get_float(rhs); CREATE_SINGLE_ITEM(RHS, float, NDARRAY_FLOAT, fvalue); } else { RHS = rhs; rhs_is_scalar = false; } //else if(mp_obj_is_type(lhs, &ulab_ndarray_type) && mp_obj_is_type(RHS, &ulab_ndarray_type)) { // next, the ndarray stuff ndarray_obj_t *ol = MP_OBJ_TO_PTR(lhs); ndarray_obj_t *or = MP_OBJ_TO_PTR(RHS); if(!rhs_is_scalar && ((ol->m != or->m) || (ol->n != or->n))) { mp_raise_ValueError("operands could not be broadcast together"); } // At this point, the operands should have the same shape switch(op) { case MP_BINARY_OP_EQUAL: // Two arrays are equal, if their shape, typecode, and elements are equal if((ol->m != or->m) || (ol->n != or->n) || (ol->array->typecode != or->array->typecode)) { return mp_const_false; } else { size_t i = ol->bytes; uint8_t *l = (uint8_t *)ol->array->items; uint8_t *r = (uint8_t *)or->array->items; while(i) { // At this point, we can simply compare the bytes, the type is irrelevant if(*l++ != *r++) { return mp_const_false; } i--; } return mp_const_true; } break; case MP_BINARY_OP_LESS: case MP_BINARY_OP_LESS_EQUAL: case MP_BINARY_OP_MORE: case MP_BINARY_OP_MORE_EQUAL: case MP_BINARY_OP_ADD: case MP_BINARY_OP_SUBTRACT: case MP_BINARY_OP_TRUE_DIVIDE: case MP_BINARY_OP_MULTIPLY: // These are the upcasting rules // float always becomes float // operation on identical types preserves type // uint8 + int8 => int16 // uint8 + int16 => int16 // uint8 + uint16 => uint16 // int8 + int16 => int16 // int8 + uint16 => uint16 // uint16 + int16 => float // The parameters of RUN_BINARY_LOOP are // typecode of result, type_out, type_left, type_right, lhs operand, rhs operand, operator if(ol->array->typecode == NDARRAY_UINT8) { if(or->array->typecode == NDARRAY_UINT8) { RUN_BINARY_LOOP(NDARRAY_UINT8, uint8_t, uint8_t, uint8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT8) { RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_UINT16) { RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint8_t, uint16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT16) { RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_FLOAT) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint8_t, float, ol, or, op); } } else if(ol->array->typecode == NDARRAY_INT8) { if(or->array->typecode == NDARRAY_UINT8) { RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT8) { RUN_BINARY_LOOP(NDARRAY_INT8, int8_t, int8_t, int8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_UINT16) { RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT16) { RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, int16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_FLOAT) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, int8_t, float, ol, or, op); } } else if(ol->array->typecode == NDARRAY_UINT16) { if(or->array->typecode == NDARRAY_UINT8) { RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT8) { RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, int8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_UINT16) { RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT16) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint16_t, int16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_FLOAT) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint8_t, float, ol, or, op); } } else if(ol->array->typecode == NDARRAY_INT16) { if(or->array->typecode == NDARRAY_UINT8) { RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, uint8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT8) { RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_UINT16) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, int16_t, uint16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT16) { RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_FLOAT) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, uint16_t, float, ol, or, op); } } else if(ol->array->typecode == NDARRAY_FLOAT) { if(or->array->typecode == NDARRAY_UINT8) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, uint8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT8) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, int8_t, ol, or, op); } else if(or->array->typecode == NDARRAY_UINT16) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, uint16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_INT16) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, int16_t, ol, or, op); } else if(or->array->typecode == NDARRAY_FLOAT) { RUN_BINARY_LOOP(NDARRAY_FLOAT, float, float, float, ol, or, op); } } else { // this should never happen mp_raise_TypeError("wrong input type"); } // this instruction should never be reached, but we have to make the compiler happy return MP_OBJ_NULL; default: return MP_OBJ_NULL; // op not supported } } else { mp_raise_TypeError("wrong operand type on the right hand side"); } } mp_obj_t ndarray_unary_op(mp_unary_op_t op, mp_obj_t self_in) { ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in); ndarray_obj_t *ndarray = NULL; switch (op) { case MP_UNARY_OP_LEN: if(self->m > 1) { return mp_obj_new_int(self->m); } else { return mp_obj_new_int(self->len); } break; case MP_UNARY_OP_INVERT: if(self->array->typecode == NDARRAY_FLOAT) { mp_raise_ValueError("operation is not supported for given type"); } // we can invert the content byte by byte, there is no need to distinguish // between different typecodes ndarray = MP_OBJ_TO_PTR(ndarray_copy(self_in)); uint8_t *array = (uint8_t *)ndarray->array->items; for(size_t i=0; i < self->bytes; i++) array[i] = ~array[i]; return MP_OBJ_FROM_PTR(ndarray); break; case MP_UNARY_OP_NEGATIVE: ndarray = MP_OBJ_TO_PTR(ndarray_copy(self_in)); if(self->array->typecode == NDARRAY_UINT8) { uint8_t *array = (uint8_t *)ndarray->array->items; for(size_t i=0; i < self->array->len; i++) array[i] = -array[i]; } else if(self->array->typecode == NDARRAY_INT8) { int8_t *array = (int8_t *)ndarray->array->items; for(size_t i=0; i < self->array->len; i++) array[i] = -array[i]; } else if(self->array->typecode == NDARRAY_UINT16) { uint16_t *array = (uint16_t *)ndarray->array->items; for(size_t i=0; i < self->array->len; i++) array[i] = -array[i]; } else if(self->array->typecode == NDARRAY_INT16) { int16_t *array = (int16_t *)ndarray->array->items; for(size_t i=0; i < self->array->len; i++) array[i] = -array[i]; } else { float *array = (float *)ndarray->array->items; for(size_t i=0; i < self->array->len; i++) array[i] = -array[i]; } return MP_OBJ_FROM_PTR(ndarray); break; case MP_UNARY_OP_POSITIVE: return ndarray_copy(self_in); case MP_UNARY_OP_ABS: if((self->array->typecode == NDARRAY_UINT8) || (self->array->typecode == NDARRAY_UINT16)) { return ndarray_copy(self_in); } ndarray = MP_OBJ_TO_PTR(ndarray_copy(self_in)); if((self->array->typecode == NDARRAY_INT8)) { int8_t *array = (int8_t *)ndarray->array->items; for(size_t i=0; i < self->array->len; i++) { if(array[i] < 0) array[i] = -array[i]; } } else if((self->array->typecode == NDARRAY_INT16)) { int16_t *array = (int16_t *)ndarray->array->items; for(size_t i=0; i < self->array->len; i++) { if(array[i] < 0) array[i] = -array[i]; } } else { float *array = (float *)ndarray->array->items; for(size_t i=0; i < self->array->len; i++) { if(array[i] < 0) array[i] = -array[i]; } } return MP_OBJ_FROM_PTR(ndarray); break; default: return MP_OBJ_NULL; // operator not supported } } ipython3.. code :: a = array([1, 2, 3, 4]) a < a[2] .. parsed-literal:: array([ True, True, False, False]) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([1, -2, 3], dtype=ulab.int8) print(abs(a)) a = ulab.ndarray([1, 2, 3], dtype=ulab.uint8) print(~a) a = ulab.ndarray([1, 2, 3], dtype=ulab.int8) print(~a) .. parsed-literal:: ndarray([1, 2, 3], dtype=int8) ndarray([254, 253, 252], dtype=uint8) ndarray([-2, -3, -4], dtype=int8) ipython3.. code :: a = array([1, -2, 3], dtype=int8) print(-a, +a) a = array([1, 2, 3], dtype=uint8) print(-a, +a) a = array([1, 2, -3], dtype=float) print(-a, +a) .. parsed-literal:: [-1 2 -3] [ 1 -2 3] [255 254 253] [1 2 3] [-1. -2. 3.] [ 1. 2. -3.] ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([1, -2, 3], dtype=ulab.int8) print(-a, +a) a = ulab.ndarray([1, 2, 3], dtype=ulab.uint8) print(-a, +a) a = ulab.ndarray([1, 2, -3], dtype=ulab.float) print(-a, +a) .. parsed-literal:: ndarray([-1, 2, -3], dtype=int8) ndarray([1, -2, 3], dtype=int8) ndarray([255, 254, 253], dtype=uint8) ndarray([1, 2, 3], dtype=uint8) ndarray([-1.0, -2.0, 3.0], dtype=float) ndarray([1.0, 2.0, -3.0], dtype=float) ipython3.. code :: a = array([1, 2, 3], dtype=int8) print(~a) .. parsed-literal:: [-2 -3 -4] ipython3.. code :: a = array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) ipython3.. code :: a[0, 1] .. parsed-literal:: 2 Linear algebra ============== 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. 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. Examples -------- Transpose of one- and two-dimensional arrays, .transpose() ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ipython3.. code :: %%micropython -unix 1 from ulab import ndarray a = ndarray(range(10)) print('1D array: ', a) print('shape of a: ', a.shape()) a.transpose() print('\ntranspose of array: ', a) print('shape of a: ', a.shape()) a = ndarray([[1, 2, 3, 4], [5, 6, 7, 8]]) print('\n2D array: \n', a) print('shape of a: ', a.shape()) a.transpose() print('\ntranspose of array: \n', a) print('shape of a: ', a.shape()) .. parsed-literal:: 1D array: ndarray([0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0], dtype=float) shape of a: (10, 1) transpose of array: ndarray([0.0, 1.0, 2.0, ..., 7.0, 8.0, 9.0], dtype=float) shape of a: (1, 10) 2D array: ndarray([[1.0, 2.0, 3.0, 4.0], [5.0, 6.0, 7.0, 8.0]], dtype=float) shape of a: (2, 4) transpose of array: ndarray([[1.0, 5.0], [2.0, 6.0], [3.0, 7.0], [4.0, 8.0]], dtype=float) shape of a: (4, 2) .reshape() ~~~~~~~~~~ ipython3.. code :: %%micropython -unix 1 from ulab import ndarray a = ndarray(range(15)) print('1D array: ', a) print('shape of a: ', a.shape()) a.reshape((3, 5)) print('\n2D array: \n', a) print('shape of a: ', a.shape()) .. parsed-literal:: 1D array: ndarray([0.0, 1.0, 2.0, ..., 12.0, 13.0, 14.0], dtype=float) shape of a: (15, 1) 2D array: ndarray([[0.0, 1.0, 2.0, 3.0, 4.0], [5.0, 6.0, 7.0, 8.0, 9.0], [10.0, 11.0, 12.0, 13.0, 14.0]], dtype=float) shape of a: (3, 5) inverse of a matrix (inv) ~~~~~~~~~~~~~~~~~~~~~~~~~ ipython3.. code :: %%micropython -unix 1 from ulab import ndarray, inv a = ndarray([[1, 2], [3, 4]]) print('2D matrix (a): \n', a) b = inv(a) print('\ninverse of a: \n', b) .. parsed-literal:: 2D matrix (a): ndarray([[1.0, 2.0], [3.0, 4.0]], dtype=float) inverse of a: ndarray([[-2.0, 1.0], [1.5, -0.5]], dtype=float) matrix multiplication (dot) ~~~~~~~~~~~~~~~~~~~~~~~~~~~ With the ``dot`` function, we can now check, whether the inverse of the matrix was correct: ipython3.. code :: %%micropython -unix 1 from ulab import ndarray, inv, dot a = ndarray([[1, 2], [3, 4]]) print('2D matrix (a): \n', a) b = inv(a) print('\ninverse of a: \n', b) c = dot(a, b) print('\na multiplied by its inverse: \n', c) .. parsed-literal:: 2D matrix (a): ndarray([[1.0, 2.0], [3.0, 4.0]], dtype=float) inverse of a: ndarray([[-2.0, 1.0], [1.5, -0.5]], dtype=float) a multiplied by its inverse: ndarray([[1.0, 0.0], [0.0, 1.0]], dtype=float) zeros, ones, eye ~~~~~~~~~~~~~~~~ ipython3.. code :: %%micropython -unix 1 import ulab print(ulab.zeros(3, dtype=ulab.int16)) print(ulab.zeros((5, 3), dtype=ulab.float)) print("\n====================\n"); print(ulab.ones(3, dtype=ulab.int16)) print(ulab.ones((5, 3), dtype=ulab.float)) print("\n====================\n"); print(ulab.eye(5, dtype=ulab.int16)) print(ulab.eye(5, M=3, dtype=ulab.float)) print(ulab.eye(5, k=1, dtype=ulab.uint8)) print(ulab.eye(5, k=-3, dtype=ulab.uint8)) .. parsed-literal:: ndarray([0, 0, 0], dtype=int16) ndarray([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], dtype=float) ==================== ndarray([1, 1, 1], dtype=int16) ndarray([[1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 1.0, 1.0]], dtype=float) ==================== ndarray([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1]], dtype=int16) ndarray([[1.0, 0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0, 0.0]], dtype=float) ndarray([[0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1], [0, 0, 0, 0, 0]], dtype=uint8) ndarray([[0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [1, 0, 0, 0, 0], [0, 1, 0, 0, 0]], dtype=uint8) linalg.h -------- https://github.com/v923z/micropython-ulab/tree/master/code/linalg.h .. code:: cpp #ifndef _LINALG_ #define _LINALG_ #include "ndarray.h" #define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; } #define epsilon 1e-6 mp_obj_t linalg_transpose(mp_obj_t ); mp_obj_t linalg_reshape(mp_obj_t , mp_obj_t ); bool linalg_invert_matrix(float *, size_t ); mp_obj_t linalg_inv(mp_obj_t ); mp_obj_t linalg_dot(mp_obj_t , mp_obj_t ); mp_obj_t linalg_zeros(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t linalg_ones(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t linalg_eye(size_t , const mp_obj_t *, mp_map_t *); #endif linalg.c -------- https://github.com/v923z/micropython-ulab/tree/master/code/linalg.c .. code:: cpp #include #include #include "py/obj.h" #include "py/runtime.h" #include "py/misc.h" #include "linalg.h" mp_obj_t linalg_transpose(mp_obj_t self_in) { ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in); // the size of a single item in the array uint8_t _sizeof = mp_binary_get_size('@', self->array->typecode, NULL); // NOTE: In principle, we could simply specify the stride direction, and then we wouldn't // even have to shuffle the elements. The downside of that approach is that we would have // to implement two versions of the matrix multiplication and inversion functions // NOTE: // if the matrices are square, we can simply swap items, but // generic matrices can't be transposed in place, so we have to // declare a temporary variable // NOTE: // In the old matrix, the coordinate (m, n) is m*self->n + n // We have to assign this to the coordinate (n, m) in the new // matrix, i.e., to n*self->m + m // one-dimensional arrays can be transposed by simply swapping the dimensions if((self->m != 1) && (self->n != 1)) { uint8_t *c = (uint8_t *)self->array->items; // self->bytes is the size of the bytearray, irrespective of the typecode uint8_t *tmp = m_new(uint8_t, self->bytes); for(size_t m=0; m < self->m; m++) { for(size_t n=0; n < self->n; n++) { memcpy(tmp+_sizeof*(n*self->m + m), c+_sizeof*(m*self->n + n), _sizeof); } } memcpy(self->array->items, tmp, self->bytes); m_del(uint8_t, tmp, self->bytes); } SWAP(size_t, self->m, self->n); return mp_const_none; } mp_obj_t linalg_reshape(mp_obj_t self_in, mp_obj_t shape) { ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in); if(!MP_OBJ_IS_TYPE(shape, &mp_type_tuple) || (MP_OBJ_SMALL_INT_VALUE(mp_obj_len_maybe(shape)) != 2)) { mp_raise_ValueError("shape must be a 2-tuple"); } mp_obj_iter_buf_t iter_buf; mp_obj_t item, iterable = mp_getiter(shape, &iter_buf); uint16_t m, n; item = mp_iternext(iterable); m = mp_obj_get_int(item); item = mp_iternext(iterable); n = mp_obj_get_int(item); if(m*n != self->m*self->n) { // TODO: the proper error message would be "cannot reshape array of size %d into shape (%d, %d)" mp_raise_ValueError("cannot reshape array (incompatible input/output shape)"); } self->m = m; self->n = n; return MP_OBJ_FROM_PTR(self); } bool linalg_invert_matrix(float *data, size_t N) { // returns true, of the inversion was successful, // false, if the matrix is singular // initially, this is the unit matrix: the contents of this matrix is what // will be returned after all the transformations float *unit = m_new(float, N*N); float elem = 1.0; // initialise the unit matrix memset(unit, 0, sizeof(float)*N*N); for(size_t m=0; m < N; m++) { memcpy(&unit[m*(N+1)], &elem, sizeof(float)); } for(size_t m=0; m < N; m++){ // this could be faster with ((c < epsilon) && (c > -epsilon)) if(abs(data[m*(N+1)]) < epsilon) { m_del(float, unit, N*N); return false; } for(size_t n=0; n < N; n++){ if(m != n){ elem = data[N*n+m] / data[m*(N+1)]; for(size_t k=0; k < N; k++){ data[N*n+k] -= elem * data[N*m+k]; unit[N*n+k] -= elem * unit[N*m+k]; } } } } for(size_t m=0; m < N; m++){ elem = data[m*(N+1)]; for(size_t n=0; n < N; n++){ data[N*m+n] /= elem; unit[N*m+n] /= elem; } } memcpy(data, unit, sizeof(float)*N*N); m_del(float, unit, N*N); return true; } mp_obj_t linalg_inv(mp_obj_t o_in) { // since inv is not a class method, we have to inspect the input argument first if(!MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) { mp_raise_TypeError("only ndarrays can be inverted"); } ndarray_obj_t *o = MP_OBJ_TO_PTR(o_in); if(!MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) { mp_raise_TypeError("only ndarray objects can be inverted"); } if(o->m != o->n) { mp_raise_ValueError("only square matrices can be inverted"); } ndarray_obj_t *inverted = create_new_ndarray(o->m, o->n, NDARRAY_FLOAT); float *data = (float *)inverted->array->items; mp_obj_t elem; for(size_t m=0; m < o->m; m++) { // rows first for(size_t n=0; n < o->n; n++) { // columns next // this could, perhaps, be done in single line... // On the other hand, we probably spend little time here elem = mp_binary_get_val_array(o->array->typecode, o->array->items, m*o->n+n); data[m*o->n+n] = (float)mp_obj_get_float(elem); } } if(!linalg_invert_matrix(data, o->m)) { // TODO: I am not sure this is needed here. Otherwise, // how should we free up the unused RAM of inverted? m_del(float, inverted->array->items, o->n*o->n); mp_raise_ValueError("input matrix is singular"); } return MP_OBJ_FROM_PTR(inverted); } mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) { // TODO: should the results be upcast? ndarray_obj_t *m1 = MP_OBJ_TO_PTR(_m1); ndarray_obj_t *m2 = MP_OBJ_TO_PTR(_m2); if(m1->n != m2->m) { mp_raise_ValueError("matrix dimensions do not match"); } ndarray_obj_t *out = create_new_ndarray(m1->m, m2->n, NDARRAY_FLOAT); float *outdata = (float *)out->array->items; float sum, v1, v2; for(size_t i=0; i < m1->n; i++) { for(size_t j=0; j < m2->m; j++) { sum = 0.0; for(size_t k=0; k < m1->m; k++) { // (j, k) * (k, j) v1 = ndarray_get_float_value(m1->array->items, m1->array->typecode, i*m1->n+k); v2 = ndarray_get_float_value(m2->array->items, m2->array->typecode, k*m2->n+j); sum += v1 * v2; } outdata[i*m1->m+j] = sum; } } return MP_OBJ_FROM_PTR(out); } mp_obj_t linalg_zeros_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t kind) { static const mp_arg_t allowed_args[] = { { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_obj = MP_OBJ_NULL} } , { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} }, }; mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); uint8_t dtype = args[1].u_int; if(!mp_obj_is_int(args[0].u_obj) && !mp_obj_is_type(args[0].u_obj, &mp_type_tuple)) { mp_raise_TypeError("input argument must be an integer or a 2-tuple"); } ndarray_obj_t *ndarray = NULL; if(mp_obj_is_int(args[0].u_obj)) { size_t n = mp_obj_get_int(args[0].u_obj); ndarray = create_new_ndarray(1, n, dtype); } else if(mp_obj_is_type(args[0].u_obj, &mp_type_tuple)) { mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(args[0].u_obj); if(tuple->len != 2) { mp_raise_TypeError("input argument must be an integer or a 2-tuple"); } ndarray = create_new_ndarray(mp_obj_get_int(tuple->items[0]), mp_obj_get_int(tuple->items[1]), dtype); } if(kind == 1) { mp_obj_t one = mp_obj_new_int(1); for(size_t i=0; i < ndarray->array->len; i++) { mp_binary_set_val_array(dtype, ndarray->array->items, i, one); } } return MP_OBJ_FROM_PTR(ndarray); } mp_obj_t linalg_zeros(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { return linalg_zeros_ones(n_args, pos_args, kw_args, 0); } mp_obj_t linalg_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { return linalg_zeros_ones(n_args, pos_args, kw_args, 1); } mp_obj_t linalg_eye(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { static const mp_arg_t allowed_args[] = { { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_INT, {.u_int = 0} }, { MP_QSTR_M, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } }, { MP_QSTR_k, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} }, { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} }, }; mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); size_t n = args[0].u_int, m; int16_t k = args[2].u_int; uint8_t dtype = args[3].u_int; if(args[1].u_rom_obj == mp_const_none) { m = n; } else { m = mp_obj_get_int(args[1].u_rom_obj); } ndarray_obj_t *ndarray = create_new_ndarray(m, n, dtype); mp_obj_t one = mp_obj_new_int(1); size_t i = 0; if((k >= 0) && (k < n)) { while(k < n) { mp_binary_set_val_array(dtype, ndarray->array->items, i*n+k, one); k++; i++; } } else if((k < 0) && (-k < m)) { k = -k; i = 0; while(k < m) { mp_binary_set_val_array(dtype, ndarray->array->items, k*n+i, one); k++; i++; } } return MP_OBJ_FROM_PTR(ndarray); } Vectorising mathematical operations =================================== General comments ---------------- The following module implements the common mathematical functions for scalars, ndarrays (linear or matrix), and iterables. If the input argument is a scalar, a scalar is returned (i.e., for such arguments, these functions are identical to the functions in the ``math`` module), while for ndarrays, and iterables, the return value is an ndarray of type ``float``. Examples -------- ipython3.. code :: %%micropython -unix 1 import ulab # initialise an array a = ulab.ndarray([1, 2, 3, 4, 5]) print('1D array: ', a) print('\nexponent of an array (range(5)): ', ulab.exp(range(5))) print('\nexponent of a scalar (2.0): ', ulab.exp(2.0)) print('\n exponent of a 1D ndarray (a): ', ulab.exp(a)) # initialise a matrix b = ulab.ndarray([[1, 2, 3], [4, 5, 6]]) print('\n2D matrix: ', b) print('exponent of a 2D matrix (b): ', ulab.exp(b)) .. parsed-literal:: 1D array: ndarray([1.0, 2.0, 3.0, 4.0, 5.0], dtype=float) exponent of an array (range(5)): ndarray([1.0, 2.718281745910645, 7.389056205749512, 20.08553695678711, 54.59814834594727], dtype=float) exponent of a scalar (2.0): 7.38905609893065 exponent of a 1D ndarray (a): ndarray([2.718281745910645, 7.389056205749512, 20.08553695678711, 54.59814834594727, 148.4131622314453], dtype=float) 2D matrix: ndarray([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], dtype=float) exponent of a 2D matrix (b): ndarray([[2.718281745910645, 7.389056205749512, 20.08553695678711], [54.59814834594727, 148.4131622314453, 403.4288024902343]], dtype=float) Note that ndarrays are linear arrays in memory, even if the ``shape`` of the ndarray is a matrix. This means that we can treat both cases in a *single* loop. Since ``ndarray``\ s are iterable, we could treat ``ndarray``\ s, ``list``\ s, ``tuples``, and ``range``\ s on the same footing. However, that would mean extra trips to a lot of functions, therefore, reading out the values of the ``ndarray`` directly is probably significantly faster. vectorise.h ----------- https://github.com/v923z/micropython-ulab/tree/master/code/vectorise.h .. code:: cpp #ifndef _VECTORISE_ #define _VECTORISE_ #include "ndarray.h" mp_obj_t vectorise_acos(mp_obj_t ); mp_obj_t vectorise_acosh(mp_obj_t ); mp_obj_t vectorise_asin(mp_obj_t ); mp_obj_t vectorise_asinh(mp_obj_t ); mp_obj_t vectorise_atan(mp_obj_t ); mp_obj_t vectorise_atanh(mp_obj_t ); mp_obj_t vectorise_ceil(mp_obj_t ); mp_obj_t vectorise_cos(mp_obj_t ); mp_obj_t vectorise_erf(mp_obj_t ); mp_obj_t vectorise_erfc(mp_obj_t ); mp_obj_t vectorise_exp(mp_obj_t ); mp_obj_t vectorise_expm1(mp_obj_t ); mp_obj_t vectorise_floor(mp_obj_t ); mp_obj_t vectorise_gamma(mp_obj_t ); mp_obj_t vectorise_lgamma(mp_obj_t ); mp_obj_t vectorise_log(mp_obj_t ); mp_obj_t vectorise_log10(mp_obj_t ); mp_obj_t vectorise_log2(mp_obj_t ); mp_obj_t vectorise_sin(mp_obj_t ); mp_obj_t vectorise_sinh(mp_obj_t ); mp_obj_t vectorise_sqrt(mp_obj_t ); mp_obj_t vectorise_tan(mp_obj_t ); mp_obj_t vectorise_tanh(mp_obj_t ); #endif vectorise.c ----------- https://github.com/v923z/micropython-ulab/tree/master/code/vectorise.c .. code:: cpp #include #include #include #include "py/runtime.h" #include "py/binary.h" #include "py/obj.h" #include "py/objarray.h" #include "vectorise.h" #ifndef MP_PI #define MP_PI MICROPY_FLOAT_CONST(3.14159265358979323846) #endif mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) { // Return a single value, if o_in is not iterable if(mp_obj_is_float(o_in) || mp_obj_is_integer(o_in)) { return mp_obj_new_float(f(mp_obj_get_float(o_in))); } mp_float_t x; if(MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) { ndarray_obj_t *o = MP_OBJ_TO_PTR(o_in); ndarray_obj_t *out = create_new_ndarray(o->m, o->n, NDARRAY_FLOAT); float *datain = (float *)o->array->items; float *dataout = (float *)out->array->items; for(size_t i=0; i < o->array->len; i++) { dataout[i] = f(datain[i]); } return MP_OBJ_FROM_PTR(out); } else if(MP_OBJ_IS_TYPE(o_in, &mp_type_tuple) || MP_OBJ_IS_TYPE(o_in, &mp_type_list) || MP_OBJ_IS_TYPE(o_in, &mp_type_range)) { mp_obj_array_t *o = MP_OBJ_TO_PTR(o_in); ndarray_obj_t *out = create_new_ndarray(1, o->len, NDARRAY_FLOAT); float *dataout = (float *)out->array->items; mp_obj_iter_buf_t iter_buf; mp_obj_t item, iterable = mp_getiter(o_in, &iter_buf); size_t i=0; while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) { x = mp_obj_get_float(item); dataout[i++] = f(x); } return MP_OBJ_FROM_PTR(out); } return mp_const_none; } #define MATH_FUN_1(py_name, c_name) \ mp_obj_t vectorise_ ## py_name(mp_obj_t x_obj) { \ return vectorise_generic_vector(x_obj, MICROPY_FLOAT_C_FUN(c_name)); \ } // _degrees won't compile for the unix port /* mp_float_t _degreesf(mp_float_t x) { return(180*x/MP_PI); } MATH_FUN_1(degrees, _degrees); // _radians won't compile for the unix port mp_float_t _radiansf(mp_float_t x) { return(MP_PI*x/180.0); } MATH_FUN_1(radians, _radians); STATIC mp_float_t _fabsf(mp_float_t x) { return fabsf(x); } MATH_FUN_1(fabs, _fabs); */ MATH_FUN_1(acos, acos); MATH_FUN_1(acosh, acosh); MATH_FUN_1(asin, asin); MATH_FUN_1(asinh, asinh); MATH_FUN_1(atan, atan); MATH_FUN_1(atanh, atanh); MATH_FUN_1(ceil, ceil); MATH_FUN_1(cos, cos); MATH_FUN_1(erf, erf); MATH_FUN_1(erfc, erfc); MATH_FUN_1(exp, exp); MATH_FUN_1(expm1, expm1); MATH_FUN_1(floor, floor); MATH_FUN_1(gamma, tgamma); MATH_FUN_1(lgamma, lgamma); MATH_FUN_1(log, log); MATH_FUN_1(log10, log10); MATH_FUN_1(log2, log2); MATH_FUN_1(sin, sin); MATH_FUN_1(sinh, sinh); MATH_FUN_1(sqrt, sqrt); MATH_FUN_1(tan, tan); MATH_FUN_1(tanh, tanh); Polynomials =========== This module has two functions, ``polyval``, and ``polyfit``. The background for ``polyfit`` can be found under https://en.wikipedia.org/wiki/Polynomial_regression, and one can take the matrix inversion function from ``linalg``. Background ---------- An estimate, :math:`\beta`, for the coefficients of a polynomial fit can be gotten from :raw-latex:`\begin{equation} \vec{\beta} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T \vec{y} \end{equation}` where :math:`\vec{y}` are the dependent values, and the matrix :math:`X` is constructed from the independent values as :raw-latex:`\begin{equation} X = \begin{pmatrix} 1 & x_1^2 & x_1^2 & ... & x_1^m \\ 1 & x_2^2 & x_2^2 & ... & x_2^m \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n^2 & x_n^2 & ... & x_n^m \end{pmatrix} \end{equation}` Note that the calculation of :math:`X^T` is trivial, and we need :math:`X` only once, namely in the product :math:`X^TX`. We will save RAM by storing only :math:`X^T`, and expressing :math:`X` from :math:`X^T`, when we need it. The routine calculates the coefficients in increasing order, therefore, before returning, we have to reverse the array. Examples -------- polyval ~~~~~~~ ipython3.. code :: %%micropython -unix 1 import ulab p = [1, 1, 1, 0] x = [0, 1, 2, 3, 4] print('coefficients: ', p) print('independent values: ', x) print('\nvalues of p(x): ', ulab.polyval(p, x)) # the same works with ndarrays a = ulab.ndarray(x) print('\nndarray (a): ', a) print('value of p(a): ', ulab.polyval(p, a)) .. parsed-literal:: coefficients: [1, 1, 1, 0] independent values: [0, 1, 2, 3, 4] values of p(x): ndarray([0.0, 3.0, 14.0, 39.0, 84.0], dtype=float) ndarray (a): ndarray([0.0, 1.0, 2.0, 3.0, 4.0], dtype=float) value of p(a): ndarray([0.0, 3.0, 14.0, 39.0, 84.0], dtype=float) polyfit ~~~~~~~ First a perfect parabola with zero shift, and leading coefficient of 1. ipython3.. code :: %%micropython -unix 1 import ulab x = ulab.ndarray([-3, -2, -1, 0, 1, 2, 3]) y = ulab.ndarray([9, 4, 1, 0, 1, 4, 9]) print('independent values: ', x) print('dependent values: ', y) print('fit values', ulab.polyfit(x, y, 2)) .. parsed-literal:: independent values: ndarray([-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0], dtype=float) dependent values: ndarray([9.0, 4.0, 1.0, 0.0, 1.0, 4.0, 9.0], dtype=float) fit values ndarray([1.00000011920929, 0.0, 0.0], dtype=float) We can now take a more meaningful example: the data points scatter here: ipython3.. code :: %%micropython -unix 1 import ulab x = ulab.ndarray([-3, -2, -1, 0, 1, 2, 3]) y = ulab.ndarray([10, 5, 1, 0, 1, 4.2, 9.1]) print('independent values: ', x) print('dependent values: ', y) print('fit values', ulab.polyfit(x, y, 2)) .. parsed-literal:: independent values: ndarray([-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0], dtype=float) dependent values: ndarray([10.0, 5.0, 1.0, 0.0, 1.0, 4.199999809265137, 9.100000381469727], dtype=float) fit values ndarray([1.065476179122925, -0.1535714119672775, 0.06666660308837891], dtype=float) Finally, let us see, what this looks like in numpy: ipython3.. code :: x = array([-3, -2, -1, 0, 1, 2, 3]) y = array([10, 5, 1, 0, 1, 4.2, 9.1]) print('independent values: ', x) print('dependent values: ', y) print('fit values: ', polyfit(x, y, 2)) .. parsed-literal:: independent values: [-3 -2 -1 0 1 2 3] dependent values: [10. 5. 1. 0. 1. 4.2 9.1] fit values: [ 1.06547619 -0.15357143 0.06666667] Look at that! poly.h ------ https://github.com/v923z/micropython-ulab/tree/master/code/poly.h .. code:: cpp #ifndef _POLY_ #define _POLY_ mp_obj_t poly_polyval(mp_obj_t , mp_obj_t ); mp_obj_t poly_polyfit(size_t , const mp_obj_t *); #endif poly.c ------ https://github.com/v923z/micropython-ulab/tree/master/code/poly.c .. code:: cpp #include "py/obj.h" #include "py/runtime.h" #include "py/objarray.h" #include "ndarray.h" #include "linalg.h" #include "poly.h" bool object_is_nditerable(mp_obj_t o_in) { if(mp_obj_is_type(o_in, &ulab_ndarray_type) || mp_obj_is_type(o_in, &mp_type_tuple) || mp_obj_is_type(o_in, &mp_type_list) || mp_obj_is_type(o_in, &mp_type_range)) { return true; } return false; } size_t get_nditerable_len(mp_obj_t o_in) { if(mp_obj_is_type(o_in, &ulab_ndarray_type)) { ndarray_obj_t *in = MP_OBJ_TO_PTR(o_in); return in->array->len; } else { return (size_t)mp_obj_get_int(mp_obj_len_maybe(o_in)); } } mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) { // TODO: return immediately, if o_p is not an iterable size_t m, n; if(MP_OBJ_IS_TYPE(o_x, &ulab_ndarray_type)) { ndarray_obj_t *ndx = MP_OBJ_TO_PTR(o_x); m = ndx->m; n = ndx->n; } else { mp_obj_array_t *ix = MP_OBJ_TO_PTR(o_x); m = 1; n = ix->len; } // polynomials are going to be of type float, except, when both // the coefficients and the independent variable are integers ndarray_obj_t *out = create_new_ndarray(m, n, NDARRAY_FLOAT); mp_obj_iter_buf_t x_buf; mp_obj_t x_item, x_iterable = mp_getiter(o_x, &x_buf); mp_obj_iter_buf_t p_buf; mp_obj_t p_item, p_iterable; mp_float_t x, y; float *outf = (float *)out->array->items; uint8_t plen = mp_obj_get_int(mp_obj_len_maybe(o_p)); float *p = m_new(float, plen); p_iterable = mp_getiter(o_p, &p_buf); uint16_t i = 0; while((p_item = mp_iternext(p_iterable)) != MP_OBJ_STOP_ITERATION) { p[i] = (float)mp_obj_get_float(p_item); i++; } i = 0; while ((x_item = mp_iternext(x_iterable)) != MP_OBJ_STOP_ITERATION) { x = mp_obj_get_float(x_item); y = p[0]; for(uint8_t j=0; j < plen-1; j++) { y *= x; y += p[j+1]; } outf[i++] = y; } m_del(float, p, plen); return MP_OBJ_FROM_PTR(out); } mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) { if((n_args != 2) && (n_args != 3)) { mp_raise_ValueError("number of arguments must be 2, or 3"); } if(!object_is_nditerable(args[0])) { mp_raise_ValueError("input data must be an iterable"); } uint16_t lenx, leny; uint8_t deg; float *x, *XT, *y, *prod; if(n_args == 2) { // only the y values are supplied // TODO: this is actually not enough: the first argument can very well be a matrix, // in which case we are between the rock and a hard place leny = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0])); deg = (uint8_t)mp_obj_get_int(args[1]); if(leny < deg) { mp_raise_ValueError("more degrees of freedom than data points"); } lenx = leny; x = m_new(float, lenx); // assume uniformly spaced data points for(size_t i=0; i < lenx; i++) { x[i] = i; } y = m_new(float, leny); fill_array_iterable(y, args[0]); } else if(n_args == 3) { lenx = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0])); leny = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0])); if(lenx != leny) { mp_raise_ValueError("input vectors must be of equal length"); } deg = (uint8_t)mp_obj_get_int(args[2]); if(leny < deg) { mp_raise_ValueError("more degrees of freedom than data points"); } x = m_new(float, lenx); fill_array_iterable(x, args[0]); y = m_new(float, leny); fill_array_iterable(y, args[1]); } // one could probably express X as a function of XT, // and thereby save RAM, because X is used only in the product XT = m_new(float, (deg+1)*leny); // XT is a matrix of shape (deg+1, len) (rows, columns) for(uint8_t i=0; i < leny; i++) { // column index XT[i+0*lenx] = 1.0; // top row for(uint8_t j=1; j < deg+1; j++) { // row index XT[i+j*leny] = XT[i+(j-1)*leny]*x[i]; } } prod = m_new(float, (deg+1)*(deg+1)); // the product matrix is of shape (deg+1, deg+1) float sum; for(uint16_t i=0; i < deg+1; i++) { // column index for(uint16_t j=0; j < deg+1; j++) { // row index sum = 0.0; for(size_t k=0; k < lenx; k++) { // (j, k) * (k, i) // Note that the second matrix is simply the transpose of the first: // X(k, i) = XT(i, k) = XT[k*lenx+i] sum += XT[j*lenx+k]*XT[i*lenx+k]; // X[k*(deg+1)+i]; } prod[j*(deg+1)+i] = sum; } } if(!linalg_invert_matrix(prod, deg+1)) { // Although X was a Vandermonde matrix, whose inverse is guaranteed to exist, // we bail out here, if prod couldn't be inverted: if the values in x are not all // distinct, prod is singular m_del(float, XT, (deg+1)*lenx); m_del(float, x, lenx); m_del(float, y, lenx); m_del(float, prod, (deg+1)*(deg+1)); mp_raise_ValueError("could not invert Vandermonde matrix"); } // at this point, we have the inverse of X^T * X // y is a column vector; x is free now, we can use it for storing intermediate values for(uint16_t i=0; i < deg+1; i++) { // row index sum = 0.0; for(uint16_t j=0; j < lenx; j++) { // column index sum += XT[i*lenx+j]*y[j]; } x[i] = sum; } // XT is no longer needed m_del(float, XT, (deg+1)*leny); ndarray_obj_t *beta = create_new_ndarray(deg+1, 1, NDARRAY_FLOAT); float *betav = (float *)beta->array->items; // x[0..(deg+1)] contains now the product X^T * y; we can get rid of y m_del(float, y, leny); // now, we calculate beta, i.e., we apply prod = (X^T * X)^(-1) on x = X^T * y; x is a column vector now for(uint16_t i=0; i < deg+1; i++) { sum = 0.0; for(uint16_t j=0; j < deg+1; j++) { sum += prod[i*(deg+1)+j]*x[j]; } betav[i] = sum; } m_del(float, x, lenx); m_del(float, prod, (deg+1)*(deg+1)); for(uint8_t i=0; i < (deg+1)/2; i++) { // We have to reverse the array, for the leading coefficient comes first. SWAP(float, betav[i], betav[deg-i]); } return MP_OBJ_FROM_PTR(beta); } Fast Fourier transform ====================== The original idea of the implementation of the fast Fourier transform is taken from Numerical recipes. The main modification is that the present FFT kernel requires two input vectors of float type: one for the real part, and one for the imaginary part, while in Numerical recipes, the real and imaginary parts occupy alternating positions in the same array. However, since ``ndarray`` cannot hold complex types, it makes sense to starts with two separate vectors. This is especially true for our particular case, since the data are most probably real, coming from an ADC or similar. By separating the real and imaginary parts at the very beginning, we can process *real* data by not providing the imaginary part. If only one argument is supplied, it is assumed to be real, and the imaginary part is automatically filled in. Now, the implementation computes the transform in place. This means that RAM space could be saved, if the old data are not required anymore. The problem, however, is that the results are of type float, irrespective of the input type. If one can somehow guarantee that the input type is also float, then the old data can be overwritten. This is what happens in the ``spectrum`` function that overwrites the input array. Examples -------- Full FFT ~~~~~~~~ ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([0, 1, 2, 3, 0, 1, 2, 3]) re, im = ulab.fft(a) print('real part: ', re) print('imag part: ', im) .. parsed-literal:: real part: ndarray([12.0, 0.0, -3.999999761581421, 0.0, -4.0, 0.0, -4.0, 0.0], dtype=float) imag part: ndarray([0.0, 0.0, 3.999999523162842, 0.0, 0.0, 0.0, -3.999999523162842, 0.0], dtype=float) The same Fourier transform on numpy: ipython3.. code :: fft.fft([0, 1, 2, 3, 0, 1, 2, 3]) .. parsed-literal:: array([12.+0.j, 0.+0.j, -4.+4.j, 0.+0.j, -4.+0.j, 0.+0.j, -4.-4.j, 0.+0.j]) Spectrum ~~~~~~~~ ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([0, 1, 2, 3, 0, 1, 2, 3]) ulab.spectrum(a) print(a) .. parsed-literal:: ndarray([12.0, 0.0, 5.656853675842285, 0.0, 4.0, 0.0, 5.656853675842285, 0.0], dtype=float) And watch this: if you need the spectrum, but do not want to overwrite your data, you can do the following ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([0, 1, 2, 3, 0, 1, 2, 3]) re, im = ulab.fft(a) print('spectrum: ', ulab.sqrt(re*re+im*im)) .. parsed-literal:: spectrum: ndarray([12.0, 0.0, 5.656853675842285, 0.0, 4.0, 0.0, 5.656853675842285, 0.0], dtype=float) fft.h ----- https://github.com/v923z/micropython-ulab/tree/master/code/fft.h .. code:: cpp #ifndef _FFT_ #define _FFT_ #ifndef MP_PI #define MP_PI MICROPY_FLOAT_CONST(3.14159265358979323846) #endif #define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; } mp_obj_t fft_fft(size_t , const mp_obj_t *); mp_obj_t fft_spectrum(mp_obj_t ); #endif fft.c ----- https://github.com/v923z/micropython-ulab/tree/master/code/fft.c .. code:: cpp #include #include #include #include #include "py/runtime.h" #include "py/binary.h" #include "py/obj.h" #include "py/objarray.h" #include "ndarray.h" #include "fft.h" void fft_kernel(float *real, float *imag, int n, int isign) { // This is basically a modification of four1 from Numerical Recipes // The main difference is that this function takes two arrays, one // for the real, and one for the imaginary parts. int j, m, mmax, istep; float tempr, tempi; float wtemp, wr, wpr, wpi, wi, theta; j = 0; for(int i = 0; i < n; i++) { if (j > i) { SWAP(float, real[i], real[j]); SWAP(float, imag[i], imag[j]); } m = n >> 1; while (j >= m && m > 0) { j -= m; m >>= 1; } j += m; } mmax = 1; while (n > mmax) { istep = mmax << 1; theta = -1.0*isign*6.28318530717959/istep; wtemp = sinf(0.5 * theta); wpr = -2.0 * wtemp * wtemp; wpi = sinf(theta); wr = 1.0; wi = 0.0; for(m = 0; m < mmax; m++) { for(int i = m; i < n; i += istep) { j = i + mmax; tempr = wr * real[j] - wi * imag[j]; tempi = wr * imag[j] + wi * real[j]; real[j] = real[i] - tempr; imag[j] = imag[i] - tempi; real[i] += tempr; imag[i] += tempi; } wtemp = wr; wr = wr*wpr - wi*wpi + wr; wi = wi*wpr + wtemp*wpi + wi; } mmax = istep; } } mp_obj_t fft_fft(size_t n_args, const mp_obj_t *args) { // TODO: return the absolute value, if keyword argument is specified // TODO: transform the data in place, if keyword argument is specified if(!MP_OBJ_IS_TYPE(args[0], &ulab_ndarray_type)) { mp_raise_NotImplementedError("FFT is defined for ndarrays only"); } if(n_args == 2) { if(!MP_OBJ_IS_TYPE(args[1], &ulab_ndarray_type)) { mp_raise_NotImplementedError("FFT is defined for ndarrays only"); } } // Check if input is of length of power of 2 ndarray_obj_t *re = MP_OBJ_TO_PTR(args[0]); uint16_t len = re->array->len; if((len & (len-1)) != 0) { mp_raise_ValueError("input array length must be power of 2"); } ndarray_obj_t *out_re = create_new_ndarray(1, len, NDARRAY_FLOAT); float *data_re = (float *)out_re->array->items; if(re->array->typecode == NDARRAY_FLOAT) { memcpy((float *)out_re->array->items, (float *)re->array->items, re->bytes); } else { for(size_t i=0; i < len; i++) { data_re[i] = ndarray_get_float_value(re->array->items, re->array->typecode, i); } } ndarray_obj_t *out_im = create_new_ndarray(1, len, NDARRAY_FLOAT); float *data_im = (float *)out_im->array->items; if(n_args == 2) { ndarray_obj_t *im = MP_OBJ_TO_PTR(args[1]); if (re->array->len != im->array->len) { mp_raise_ValueError("real and imaginary parts must be of equal length"); } if(im->array->typecode == NDARRAY_FLOAT) { memcpy((float *)out_im->array->items, (float *)im->array->items, im->bytes); } else { for(size_t i=0; i < len; i++) { data_im[i] = ndarray_get_float_value(im->array->items, im->array->typecode, i); } } } fft_kernel(data_re, data_im, len, 1); mp_obj_t tuple[2]; tuple[0] = out_re; tuple[1] = out_im; return mp_obj_new_tuple(2, tuple); } mp_obj_t fft_spectrum(mp_obj_t oin) { // calculates the the spectrum of a single real ndarray in place if(!MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) { mp_raise_NotImplementedError("FFT is defined for ndarrays only"); } ndarray_obj_t *re = MP_OBJ_TO_PTR(oin); uint16_t len = re->array->len; if((re->m > 1) && (re->n > 1)) { mp_raise_ValueError("input data must be an array"); } if((len & (len-1)) != 0) { mp_raise_ValueError("input array length must be power of 2"); } if(re->array->typecode != NDARRAY_FLOAT) { mp_raise_TypeError("input array must be of type float"); } float *data_re = (float *)re->array->items; float *data_im = m_new(float, len); fft_kernel(data_re, data_im, len, 1); for(size_t i=0; i < len; i++) { data_re[i] = sqrtf(data_re[i]*data_re[i] + data_im[i]*data_im[i]); } m_del(float, data_im, len); return mp_const_none; } Numerical ========= General comments ---------------- 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. 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``. roll ~~~~ 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. Examples -------- ipython3.. code :: %%micropython -unix 1 import ulab print(ulab.linspace(0, 10, 11)) print(ulab.sum([1, 2, 3])) a = ulab.ndarray([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [55, 66, 77, 88, 99]], dtype=ulab.int8) print(a) ulab.roll(a, -1, axis=0) print(a) ulab.roll(a, 1, axis=1) print(a) .. parsed-literal:: ndarray([0.0, 1.0, 2.0, ..., 8.0, 9.0, 10.0], dtype=float) 6.0 ndarray([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [55, 66, 77, 88, 99]], dtype=int8) ndarray([[55, 66, 77, 88, 99], [1, 2, 3, 4, 5], [6, 7, 8, 9, 10]], dtype=int8) ndarray([[66, 77, 88, 99, 55], [2, 3, 4, 5, 1], [7, 8, 9, 10, 6]], dtype=int8) ipython3.. code :: linspace(0, 10, 11, endpoint=False, dtype=int8, retstep=True) .. parsed-literal:: (array([0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=int8), 0.9090909090909091) ipython3.. code :: a = array([0, 1, 2, 3], dtype=float) abs(a), -a .. parsed-literal:: (array([0., 1., 2., 3.]), array([-0., -1., -2., -3.])) ipython3.. code :: %%micropython -unix 1 import ulab print(ulab.linspace(0, 10, num=11, endpoint=True, retstep=True, dtype=ulab.int8)) print(ulab.linspace(0, 10, num=11, endpoint=False, retstep=True, dtype=ulab.int16)) .. parsed-literal:: (ndarray([0, 1, 2, ..., 8, 9, 10], dtype=int8), 1.0) (ndarray([0, 0, 1, ..., 7, 8, 9], dtype=int16), 0.9090909361839294) ipython3.. code :: %%micropython -unix 1 import ulab a = ulab.ndarray([0, 1, 2, -3], dtype=ulab.float) print(abs(a)) .. parsed-literal:: ndarray([0.0, 1.0, 2.0, 3.0], dtype=float) numerical.h ----------- https://github.com/v923z/micropython-ulab/tree/master/code/numerical.h .. code:: cpp #ifndef _NUMERICAL_ #define _NUMERICAL_ #include "ndarray.h" mp_obj_t numerical_linspace(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t numerical_sum(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t numerical_mean(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t numerical_std(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t numerical_min(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t numerical_max(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t numerical_argmin(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t numerical_argmax(size_t , const mp_obj_t *, mp_map_t *); mp_obj_t numerical_roll(size_t , const mp_obj_t *, mp_map_t *); #endif numerical.c ----------- Parsing of arguments ~~~~~~~~~~~~~~~~~~~~ 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. https://github.com/v923z/micropython-ulab/tree/master/code/numerical.c .. code:: cpp #include #include #include #include "py/obj.h" #include "py/runtime.h" #include "py/builtin.h" #include "py/misc.h" #include "numerical.h" enum NUMERICAL_FUNCTION_TYPE { NUMERICAL_MIN, NUMERICAL_MAX, NUMERICAL_ARGMIN, NUMERICAL_ARGMAX, NUMERICAL_SUM, NUMERICAL_MEAN, NUMERICAL_STD, }; mp_obj_t numerical_linspace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { static const mp_arg_t allowed_args[] = { { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } }, { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } }, { MP_QSTR_num, MP_ARG_INT, {.u_int = 50} }, { MP_QSTR_endpoint, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_true_obj)} }, { MP_QSTR_retstep, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_false_obj)} }, { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} }, }; mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; mp_arg_parse_all(2, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); uint16_t len = args[2].u_int; if(len < 2) { mp_raise_ValueError("number of points must be at least 2"); } float value, step; value = mp_obj_get_float(args[0].u_obj); uint8_t typecode = args[5].u_int; if(args[3].u_obj == mp_const_true) step = (mp_obj_get_float(args[1].u_obj)-value)/(len-1); else step = (mp_obj_get_float(args[1].u_obj)-value)/len; ndarray_obj_t *ndarray = create_new_ndarray(1, len, typecode); if(typecode == NDARRAY_UINT8) { uint8_t *array = (uint8_t *)ndarray->array->items; for(size_t i=0; i < len; i++, value += step) array[i] = (uint8_t)value; } else if(typecode == NDARRAY_INT8) { int8_t *array = (int8_t *)ndarray->array->items; for(size_t i=0; i < len; i++, value += step) array[i] = (int8_t)value; } else if(typecode == NDARRAY_UINT16) { uint16_t *array = (uint16_t *)ndarray->array->items; for(size_t i=0; i < len; i++, value += step) array[i] = (uint16_t)value; } else if(typecode == NDARRAY_INT16) { int16_t *array = (int16_t *)ndarray->array->items; for(size_t i=0; i < len; i++, value += step) array[i] = (int16_t)value; } else { float *array = (float *)ndarray->array->items; for(size_t i=0; i < len; i++, value += step) array[i] = value; } if(args[4].u_obj == mp_const_false) { return MP_OBJ_FROM_PTR(ndarray); } else { mp_obj_t tuple[2]; tuple[0] = ndarray; tuple[1] = mp_obj_new_float(step); return mp_obj_new_tuple(2, tuple); } } mp_obj_t numerical_sum_mean_std_array(mp_obj_t oin, uint8_t optype) { mp_float_t value, sum = 0.0, sq_sum = 0.0; mp_obj_iter_buf_t iter_buf; mp_obj_t item, iterable = mp_getiter(oin, &iter_buf); mp_int_t len = mp_obj_get_int(mp_obj_len(oin)); while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) { value = mp_obj_get_float(item); sum += value; if(optype == NUMERICAL_STD) { sq_sum += value*value; } } if(optype == NUMERICAL_SUM) { return mp_obj_new_float(sum); } else if(optype == NUMERICAL_MEAN) { return mp_obj_new_float(sum/len); } else { sum /= len; // this is now the mean! return mp_obj_new_float(sqrtf((sq_sum/len-sum*sum))); } } STATIC mp_obj_t numerical_argmin_argmax_array(mp_obj_t o_in, mp_uint_t op, uint8_t type) { size_t idx = 0, best_idx = 0; mp_obj_iter_buf_t iter_buf; mp_obj_t iterable = mp_getiter(o_in, &iter_buf); mp_obj_t best_obj = MP_OBJ_NULL; mp_obj_t item; while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) { if ((best_obj == MP_OBJ_NULL) || (mp_binary_op(op, item, best_obj) == mp_const_true)) { best_obj = item; best_idx = idx; } idx++; } if((type == NUMERICAL_ARGMIN) || (type == NUMERICAL_ARGMAX)) { return MP_OBJ_NEW_SMALL_INT(best_idx); } else { return best_obj; } } STATIC size_t numerical_argmin_argmax_single_line(void *data, size_t start, size_t stop, size_t stride, uint8_t typecode, uint8_t optype) { size_t best_idx = start; mp_float_t value, best_value = ndarray_get_float_value(data, typecode, start); for(size_t i=start; i < stop; i+=stride) { value = ndarray_get_float_value(data, typecode, i); if((optype == NUMERICAL_MIN) || (optype == NUMERICAL_ARGMIN)) { if(best_value > value) { best_value = value; best_idx = i; } } else if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_ARGMAX)) { if(best_value < value) { best_value = value; best_idx = i; } } } return best_idx; } STATIC mp_obj_t numerical_argmin_argmax_matrix(mp_obj_t oin, mp_obj_t axis, uint8_t optype) { ndarray_obj_t *in = MP_OBJ_TO_PTR(oin); size_t best_idx; if((axis == mp_const_none) || (in->m == 1) || (in->n == 1)) { // return the value for the flattened array best_idx = numerical_argmin_argmax_single_line(in->array->items, 0, in->array->len, 1, in->array->typecode, optype); if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) { return MP_OBJ_NEW_SMALL_INT(best_idx); } else { // TODO: do we have to do type conversion here, depending on the type of the input array? return mp_obj_new_float(ndarray_get_float_value(in->array->items, in->array->typecode, best_idx)); } } else { uint8_t _axis = mp_obj_get_int(axis); size_t m = (_axis == 0) ? 1 : in->m; size_t n = (_axis == 0) ? in->n : 1; size_t len = in->array->len; // TODO: pass in->array->typcode to create_new_ndarray ndarray_obj_t *out = create_new_ndarray(m, n, NDARRAY_FLOAT); // TODO: these two cases could probably be combined in a more elegant fashion... if(_axis == 0) { // vertical for(size_t i=0; i < n; i++) { best_idx = numerical_argmin_argmax_single_line(in->array->items, i, len, n, in->array->typecode, optype); if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) { ((float_t *)out->array->items)[i] = (float)best_idx; } else { ((float_t *)out->array->items)[i] = ndarray_get_float_value(in->array->items, in->array->typecode, best_idx); } } } else { // horizontal for(size_t i=0; i < m; i++) { best_idx = numerical_argmin_argmax_single_line(in->array->items, i*in->n, (i+1)*in->n, 1, in->array->typecode, optype); if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) { ((float_t *)out->array->items)[i] = (float)best_idx; } else { ((float_t *)out->array->items)[i] = ndarray_get_float_value(in->array->items, in->array->typecode, best_idx); } } } return MP_OBJ_FROM_PTR(out); } return mp_const_none; } STATIC mp_float_t numerical_sum_mean_std_single_line(void *data, size_t start, size_t stop, size_t stride, uint8_t typecode, uint8_t optype) { mp_float_t sum = 0.0, sq_sum = 0.0, value; size_t len = 0; for(size_t i=start; i < stop; i+=stride, len++) { value = ndarray_get_float_value(data, typecode, i); sum += value; if(optype == NUMERICAL_STD) { sq_sum += value*value; } } if(len == 0) { mp_raise_ValueError("data length is 0!"); } if(optype == NUMERICAL_SUM) { return sum; } else if(optype == NUMERICAL_MEAN) { return sum/len; } else { sum /= len; // this is now the mean! return sqrtf((sq_sum/len-sum*sum)); } } STATIC mp_obj_t numerical_sum_mean_std_matrix(mp_obj_t oin, mp_obj_t axis, uint8_t optype) { ndarray_obj_t *in = MP_OBJ_TO_PTR(oin); if((axis == mp_const_none) || (in->m == 1) || (in->n == 1)) { // return the value for the flattened array return mp_obj_new_float(numerical_sum_mean_std_single_line(in->array->items, 0, in->array->len, 1, in->array->typecode, optype)); } else { uint8_t _axis = mp_obj_get_int(axis); size_t m = (_axis == 0) ? 1 : in->m; size_t n = (_axis == 0) ? in->n : 1; size_t len = in->array->len; mp_float_t sms; // TODO: pass in->array->typcode to create_new_ndarray ndarray_obj_t *out = create_new_ndarray(m, n, NDARRAY_FLOAT); // TODO: these two cases could probably be combined in a more elegant fashion... if(_axis == 0) { // vertical for(size_t i=0; i < n; i++) { sms = numerical_sum_mean_std_single_line(in->array->items, i, len, n, in->array->typecode, optype); ((float_t *)out->array->items)[i] = sms; } } else { // horizontal for(size_t i=0; i < m; i++) { sms = numerical_sum_mean_std_single_line(in->array->items, i*in->n, (i+1)*in->n, 1, in->array->typecode, optype); ((float_t *)out->array->items)[i] = sms; } } return MP_OBJ_FROM_PTR(out); } } STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t type) { static const mp_arg_t allowed_args[] = { { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} } , { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} }, }; mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); mp_obj_t oin = args[0].u_obj; mp_obj_t axis = args[1].u_obj; if((axis != mp_const_none) && (mp_obj_get_int(axis) != 0) && (mp_obj_get_int(axis) != 1)) { // this seems to pass with False, and True... mp_raise_ValueError("axis must be None, 0, or 1"); } if(MP_OBJ_IS_TYPE(oin, &mp_type_tuple) || MP_OBJ_IS_TYPE(oin, &mp_type_list) || MP_OBJ_IS_TYPE(oin, &mp_type_range)) { switch(type) { case NUMERICAL_MIN: case NUMERICAL_ARGMIN: return numerical_argmin_argmax_array(oin, MP_BINARY_OP_LESS, type); case NUMERICAL_MAX: case NUMERICAL_ARGMAX: return numerical_argmin_argmax_array(oin, MP_BINARY_OP_MORE, type); case NUMERICAL_SUM: case NUMERICAL_MEAN: case NUMERICAL_STD: return numerical_sum_mean_std_array(oin, type); default: // we should never reach this point, but whatever return mp_const_none; } } else if(MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) { switch(type) { case NUMERICAL_MIN: case NUMERICAL_MAX: case NUMERICAL_ARGMIN: case NUMERICAL_ARGMAX: return numerical_argmin_argmax_matrix(oin, axis, type); case NUMERICAL_SUM: case NUMERICAL_MEAN: case NUMERICAL_STD: return numerical_sum_mean_std_matrix(oin, axis, type); default: mp_raise_NotImplementedError("operation is not implemented on ndarrays"); } } else { mp_raise_TypeError("input must be tuple, list, range, or ndarray"); } return mp_const_none; } mp_obj_t numerical_min(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MIN); } mp_obj_t numerical_max(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MAX); } mp_obj_t numerical_argmin(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { return numerical_function(n_args, pos_args, kw_args, NUMERICAL_ARGMIN); } mp_obj_t numerical_argmax(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { return numerical_function(n_args, pos_args, kw_args, NUMERICAL_ARGMAX); } mp_obj_t numerical_sum(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { return numerical_function(n_args, pos_args, kw_args, NUMERICAL_SUM); } mp_obj_t numerical_mean(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MEAN); } mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { return numerical_function(n_args, pos_args, kw_args, NUMERICAL_STD); } mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { // TODO: replace memcpy by memmove static const mp_arg_t allowed_args[] = { { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj)} }, { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&mp_const_none_obj) } }, { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} }, }; mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; mp_arg_parse_all(2, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); mp_obj_t oin = args[0].u_obj; int16_t shift = mp_obj_get_int(args[1].u_obj); int8_t axis = args[2].u_int; if((axis != 0) && (axis != 1)) { mp_raise_ValueError("axis must be None, 0, or 1"); } ndarray_obj_t *in = MP_OBJ_TO_PTR(oin); uint8_t _sizeof = mp_binary_get_size('@', in->array->typecode, NULL); size_t len; int16_t _shift; uint8_t *data = (uint8_t *)in->array->items; // TODO: transpose the matrix, if axis == 0 if(shift < 0) { _shift = -shift; } else { _shift = shift; } if(axis == 0) { len = in->m; // temporary buffer uint8_t *_data = m_new(uint8_t, _sizeof*len); _shift = _shift % len; if(shift < 0) _shift = len - _shift; _shift *= _sizeof; uint8_t *tmp = m_new(uint8_t, _shift); for(size_t n=0; n < in->n; n++) { for(size_t m=0; m < len; m++) { // this loop should fill up the temporary buffer memcpy(&_data[m*_sizeof], &data[(m*in->n+n)*_sizeof], _sizeof); } // now, the actual shift memcpy(tmp, _data, _shift); memcpy(_data, &_data[_shift], len*_sizeof-_shift); memcpy(&_data[len*_sizeof-_shift], tmp, _shift); for(size_t m=0; m < len; m++) { // this loop should dump the content of the temporary buffer into data memcpy(&data[(m*in->n+n)*_sizeof], &_data[m*_sizeof], _sizeof); } } m_del(uint8_t, tmp, _shift); m_del(uint8_t, _data, _sizeof*len); return mp_const_none; } len = in->n; if((in->m == 1) || (in->n == 1)) { len = in->array->len; } _shift = _shift % len; if(shift < 0) _shift = len - _shift; // TODO: if(shift > len/2), we should move in the opposite direction. That would save RAM _shift *= _sizeof; uint8_t *tmp = m_new(uint8_t, _shift); for(size_t m=0; m < in->m; m++) { memcpy(tmp, &data[m*len*_sizeof], _shift); memcpy(&data[m*len*_sizeof], &data[m*len*_sizeof+_shift], len*_sizeof-_shift); memcpy(&data[(m+1)*len*_sizeof-_shift], tmp, _shift); } m_del(uint8_t, tmp, _shift); return mp_const_none; } ulab module =========== This module simply brings all components together, and does not contain new function definitions. ulab.c ------ https://github.com/v923z/micropython-ulab/tree/master/code/ulab.c .. code:: cpp #include #include #include #include #include "py/runtime.h" #include "py/binary.h" #include "py/obj.h" #include "py/objarray.h" // this can in the future be dropped #include "ndarray.h" #include "linalg.h" #include "vectorise.h" #include "poly.h" #include "fft.h" #include "numerical.h" #define ULAB_VERSION 0.15 typedef struct _mp_obj_float_t { mp_obj_base_t base; mp_float_t value; } mp_obj_float_t; mp_obj_float_t ulab_version = {{&mp_type_float}, ULAB_VERSION}; MP_DEFINE_CONST_FUN_OBJ_1(ndarray_shape_obj, ndarray_shape); MP_DEFINE_CONST_FUN_OBJ_2(ndarray_size_obj, ndarray_size); MP_DEFINE_CONST_FUN_OBJ_1(ndarray_rawsize_obj, ndarray_rawsize); MP_DEFINE_CONST_FUN_OBJ_KW(ndarray_flatten_obj, 1, ndarray_flatten); MP_DEFINE_CONST_FUN_OBJ_1(linalg_transpose_obj, linalg_transpose); MP_DEFINE_CONST_FUN_OBJ_2(linalg_reshape_obj, linalg_reshape); MP_DEFINE_CONST_FUN_OBJ_1(linalg_inv_obj, linalg_inv); MP_DEFINE_CONST_FUN_OBJ_2(linalg_dot_obj, linalg_dot); MP_DEFINE_CONST_FUN_OBJ_KW(linalg_zeros_obj, 0, linalg_zeros); MP_DEFINE_CONST_FUN_OBJ_KW(linalg_ones_obj, 0, linalg_ones); MP_DEFINE_CONST_FUN_OBJ_KW(linalg_eye_obj, 0, linalg_eye); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acos_obj, vectorise_acos); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_acosh_obj, vectorise_acosh); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_asin_obj, vectorise_asin); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_asinh_obj, vectorise_asinh); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_atan_obj, vectorise_atan); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_atanh_obj, vectorise_atanh); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_ceil_obj, vectorise_ceil); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_cos_obj, vectorise_cos); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_erf_obj, vectorise_erf); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_erfc_obj, vectorise_erfc); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_exp_obj, vectorise_exp); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_expm1_obj, vectorise_expm1); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_floor_obj, vectorise_floor); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_gamma_obj, vectorise_gamma); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_lgamma_obj, vectorise_lgamma); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log_obj, vectorise_log); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log10_obj, vectorise_log10); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_log2_obj, vectorise_log2); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sin_obj, vectorise_sin); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sinh_obj, vectorise_sinh); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_sqrt_obj, vectorise_sqrt); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tan_obj, vectorise_tan); MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tanh_obj, vectorise_tanh); STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_linspace_obj, 2, numerical_linspace); STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sum_obj, 1, numerical_sum); STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_mean_obj, 1, numerical_mean); STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_std_obj, 1, numerical_std); STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_min_obj, 1, numerical_min); STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_max_obj, 1, numerical_max); STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argmin_obj, 1, numerical_argmin); STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argmax_obj, 1, numerical_argmax); STATIC MP_DEFINE_CONST_FUN_OBJ_KW(numerical_roll_obj, 2, numerical_roll); STATIC MP_DEFINE_CONST_FUN_OBJ_2(poly_polyval_obj, poly_polyval); STATIC MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(poly_polyfit_obj, 2, 3, poly_polyfit); STATIC MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj, 1, 2, fft_fft); STATIC MP_DEFINE_CONST_FUN_OBJ_1(fft_spectrum_obj, fft_spectrum); STATIC const mp_rom_map_elem_t ulab_ndarray_locals_dict_table[] = { { MP_ROM_QSTR(MP_QSTR_shape), MP_ROM_PTR(&ndarray_shape_obj) }, { MP_ROM_QSTR(MP_QSTR_size), MP_ROM_PTR(&ndarray_size_obj) }, { MP_ROM_QSTR(MP_QSTR_rawsize), MP_ROM_PTR(&ndarray_rawsize_obj) }, { MP_ROM_QSTR(MP_QSTR_flatten), MP_ROM_PTR(&ndarray_flatten_obj) }, { MP_ROM_QSTR(MP_QSTR_transpose), MP_ROM_PTR(&linalg_transpose_obj) }, { MP_ROM_QSTR(MP_QSTR_reshape), MP_ROM_PTR(&linalg_reshape_obj) }, }; STATIC MP_DEFINE_CONST_DICT(ulab_ndarray_locals_dict, ulab_ndarray_locals_dict_table); const mp_obj_type_t ulab_ndarray_type = { { &mp_type_type }, .name = MP_QSTR_ndarray, .print = ndarray_print, .make_new = ndarray_make_new, .subscr = ndarray_subscr, .getiter = ndarray_getiter, .unary_op = ndarray_unary_op, .binary_op = ndarray_binary_op, .locals_dict = (mp_obj_dict_t*)&ulab_ndarray_locals_dict, }; STATIC const mp_map_elem_t ulab_globals_table[] = { { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_ulab) }, { MP_ROM_QSTR(MP_QSTR___version__), MP_ROM_PTR(&ulab_version) }, { MP_OBJ_NEW_QSTR(MP_QSTR_ndarray), (mp_obj_t)&ulab_ndarray_type }, { MP_OBJ_NEW_QSTR(MP_QSTR_inv), (mp_obj_t)&linalg_inv_obj }, { MP_ROM_QSTR(MP_QSTR_dot), (mp_obj_t)&linalg_dot_obj }, { MP_ROM_QSTR(MP_QSTR_zeros), (mp_obj_t)&linalg_zeros_obj }, { MP_ROM_QSTR(MP_QSTR_ones), (mp_obj_t)&linalg_ones_obj }, { MP_ROM_QSTR(MP_QSTR_eye), (mp_obj_t)&linalg_eye_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_acos), (mp_obj_t)&vectorise_acos_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_acosh), (mp_obj_t)&vectorise_acosh_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_asin), (mp_obj_t)&vectorise_asin_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_asinh), (mp_obj_t)&vectorise_asinh_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_atan), (mp_obj_t)&vectorise_atan_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_atanh), (mp_obj_t)&vectorise_atanh_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_ceil), (mp_obj_t)&vectorise_ceil_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_cos), (mp_obj_t)&vectorise_cos_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_erf), (mp_obj_t)&vectorise_erf_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_erfc), (mp_obj_t)&vectorise_erfc_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_exp), (mp_obj_t)&vectorise_exp_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_expm1), (mp_obj_t)&vectorise_expm1_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_floor), (mp_obj_t)&vectorise_floor_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_gamma), (mp_obj_t)&vectorise_gamma_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_lgamma), (mp_obj_t)&vectorise_lgamma_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_log), (mp_obj_t)&vectorise_log_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_log10), (mp_obj_t)&vectorise_log10_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_log2), (mp_obj_t)&vectorise_log2_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_sin), (mp_obj_t)&vectorise_sin_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_sinh), (mp_obj_t)&vectorise_sinh_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_sqrt), (mp_obj_t)&vectorise_sqrt_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_tan), (mp_obj_t)&vectorise_tan_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_tanh), (mp_obj_t)&vectorise_tanh_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_linspace), (mp_obj_t)&numerical_linspace_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_sum), (mp_obj_t)&numerical_sum_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_mean), (mp_obj_t)&numerical_mean_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_std), (mp_obj_t)&numerical_std_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_min), (mp_obj_t)&numerical_min_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_max), (mp_obj_t)&numerical_max_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_argmin), (mp_obj_t)&numerical_argmin_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_argmax), (mp_obj_t)&numerical_argmax_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_roll), (mp_obj_t)&numerical_roll_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_polyval), (mp_obj_t)&poly_polyval_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_polyfit), (mp_obj_t)&poly_polyfit_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_fft), (mp_obj_t)&fft_fft_obj }, { MP_OBJ_NEW_QSTR(MP_QSTR_spectrum), (mp_obj_t)&fft_spectrum_obj }, // class constants { MP_ROM_QSTR(MP_QSTR_uint8), MP_ROM_INT(NDARRAY_UINT8) }, { MP_ROM_QSTR(MP_QSTR_int8), MP_ROM_INT(NDARRAY_INT8) }, { MP_ROM_QSTR(MP_QSTR_uint16), MP_ROM_INT(NDARRAY_UINT16) }, { MP_ROM_QSTR(MP_QSTR_int16), MP_ROM_INT(NDARRAY_INT16) }, { MP_ROM_QSTR(MP_QSTR_float), MP_ROM_INT(NDARRAY_FLOAT) }, }; STATIC MP_DEFINE_CONST_DICT ( mp_module_ulab_globals, ulab_globals_table ); const mp_obj_module_t ulab_user_cmodule = { .base = { &mp_type_module }, .globals = (mp_obj_dict_t*)&mp_module_ulab_globals, }; MP_REGISTER_MODULE(MP_QSTR_ulab, ulab_user_cmodule, MODULE_ULAB_ENABLED); makefile -------- ipython3.. code :: %%writefile ../../../ulab/code/micropython.mk USERMODULES_DIR := $(USERMOD_DIR) # Add all C files to SRC_USERMOD. SRC_USERMOD += $(USERMODULES_DIR)/ndarray.c SRC_USERMOD += $(USERMODULES_DIR)/linalg.c SRC_USERMOD += $(USERMODULES_DIR)/vectorise.c SRC_USERMOD += $(USERMODULES_DIR)/poly.c SRC_USERMOD += $(USERMODULES_DIR)/fft.c SRC_USERMOD += $(USERMODULES_DIR)/numerical.c SRC_USERMOD += $(USERMODULES_DIR)/ulab.c # We can add our module folder to include paths if needed # This is not actually needed in this example. CFLAGS_USERMOD += -I$(USERMODULES_DIR) .. parsed-literal:: Overwriting ../../../ulab/code/micropython.mk make ---- unix port ~~~~~~~~~ ipython3.. code :: %cd ../../../micropython/ports/unix/ .. parsed-literal:: /home/v923z/sandbox/micropython/v1.11/micropython/ports/unix .. code:: bash !make USER_C_MODULES=../../../ulab all .. parsed-literal:: Use make V=1 or set BUILD_VERBOSE in your environment to increase build verbosity. Including User C Module from ../../../ulab/code GEN build/genhdr/moduledefs.h GEN build/genhdr/qstr.i.last GEN build/genhdr/qstr.split GEN build/genhdr/qstrdefs.collected.h QSTR not updated CC ../../py/objmodule.c CC ../../../ulab/code/ndarray.c CC ../../../ulab/code/ulab.c LINK micropython text data bss dec hex filename 2085 6862 0 8947 22f3 build/build/frozen_mpy.o 2 0 0 2 2 build/build/frozen.o 477228 57856 2104 537188 83264 micropython stm32 port ~~~~~~~~~~ ipython3.. code :: %cd ../../../micropython/ports/stm32/ .. parsed-literal:: /home/v923z/sandbox/micropython/v1.11/micropython/ports/stm32 .. code:: bash !make BOARD=PYBV11 CROSS_COMPILE=../../../../compiler/bin/arm-none-eabi- USER_C_MODULES=../../../ulab all Change log ========== ipython3.. code :: %%writefile ../../../ulab/docs/ulab-change-log.md Tue, 8 Oct 2019 version 0.15 added inv, neg, pos, and abs unary operators to ndarray.c Mon, 7 Oct 2019 version 0.14 made the internal binary_op function tighter, and added keyword arguments to linspace Sat, 4 Oct 2019 version 0.13 added the <, <=, >, >= binary operators to ndarray Fri, 4 Oct 2019 version 0.12 added .flatten to ndarray, ones, zeros, and eye to linalg Thu, 3 Oct 2019 version 0.11 binary operators are now based on macros .. parsed-literal:: Overwriting ../../../ulab/docs/ulab-change-log.md ipython3.. code ::