path: root/circuitpython/extmod/ulab/docs/ulab-approx.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-08T12:50:51.417613Z",
     "start_time": "2021-01-08T12:50:51.208257Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook magic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-08T12:50:52.581876Z",
     "start_time": "2021-01-08T12:50:52.567901Z"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.magic import Magics, magics_class, line_cell_magic\n",
    "from IPython.core.magic import cell_magic, register_cell_magic, register_line_magic\n",
    "from IPython.core.magic_arguments import argument, magic_arguments, parse_argstring\n",
    "import subprocess\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-08T12:50:53.516712Z",
     "start_time": "2021-01-08T12:50:53.454984Z"
    }
   },
   "outputs": [],
   "source": [
    "@magics_class\n",
    "class PyboardMagic(Magics):\n",
    "    @cell_magic\n",
    "    @magic_arguments()\n",
    "    @argument('-skip')\n",
    "    @argument('-unix')\n",
    "    @argument('-pyboard')\n",
    "    @argument('-file')\n",
    "    @argument('-data')\n",
    "    @argument('-time')\n",
    "    @argument('-memory')\n",
    "    def micropython(self, line='', cell=None):\n",
    "        args = parse_argstring(self.micropython, line)\n",
    "        if args.skip: # doesn't care about the cell's content\n",
    "            print('skipped execution')\n",
    "            return None # do not parse the rest\n",
    "        if args.unix: # tests the code on the unix port. Note that this works on unix only\n",
    "            with open('/dev/shm/micropython.py', 'w') as fout:\n",
    "                fout.write(cell)\n",
    "            proc = subprocess.Popen([\"../../micropython/ports/unix/micropython\", \"/dev/shm/micropython.py\"], \n",
    "                                    stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "            print(proc.stdout.read().decode(\"utf-8\"))\n",
    "            print(proc.stderr.read().decode(\"utf-8\"))\n",
    "            return None\n",
    "        if args.file: # can be used to copy the cell content onto the pyboard's flash\n",
    "            spaces = \"    \"\n",
    "            try:\n",
    "                with open(args.file, 'w') as fout:\n",
    "                    fout.write(cell.replace('\\t', spaces))\n",
    "                    printf('written cell to {}'.format(args.file))\n",
    "            except:\n",
    "                print('Failed to write to disc!')\n",
    "            return None # do not parse the rest\n",
    "        if args.data: # can be used to load data from the pyboard directly into kernel space\n",
    "            message = pyb.exec(cell)\n",
    "            if len(message) == 0:\n",
    "                print('pyboard >>>')\n",
    "            else:\n",
    "                print(message.decode('utf-8'))\n",
    "                # register new variable in user namespace\n",
    "                self.shell.user_ns[args.data] = string_to_matrix(message.decode(\"utf-8\"))\n",
    "        \n",
    "        if args.time: # measures the time of executions\n",
    "            pyb.exec('import utime')\n",
    "            message = pyb.exec('t = utime.ticks_us()\\n' + cell + '\\ndelta = utime.ticks_diff(utime.ticks_us(), t)' + \n",
    "                               \"\\nprint('execution time: {:d} us'.format(delta))\")\n",
    "            print(message.decode('utf-8'))\n",
    "        \n",
    "        if args.memory: # prints out memory information \n",
    "            message = pyb.exec('from micropython import mem_info\\nprint(mem_info())\\n')\n",
    "            print(\"memory before execution:\\n========================\\n\", message.decode('utf-8'))\n",
    "            message = pyb.exec(cell)\n",
    "            print(\">>> \", message.decode('utf-8'))\n",
    "            message = pyb.exec('print(mem_info())')\n",
    "            print(\"memory after execution:\\n========================\\n\", message.decode('utf-8'))\n",
    "\n",
    "        if args.pyboard:\n",
    "            message = pyb.exec(cell)\n",
    "            print(message.decode('utf-8'))\n",
    "\n",
    "ip = get_ipython()\n",
    "ip.register_magics(PyboardMagic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## pyboard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-07T07:35:35.126401Z",
     "start_time": "2020-05-07T07:35:35.105824Z"
    }
   },
   "outputs": [],
   "source": [
    "import pyboard\n",
    "pyb = pyboard.Pyboard('/dev/ttyACM0')\n",
    "pyb.enter_raw_repl()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-19T19:11:18.145548Z",
     "start_time": "2020-05-19T19:11:18.137468Z"
    }
   },
   "outputs": [],
   "source": [
    "pyb.exit_raw_repl()\n",
    "pyb.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-07T07:35:38.725924Z",
     "start_time": "2020-05-07T07:35:38.645488Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "%%micropython -pyboard 1\n",
    "\n",
    "import utime\n",
    "import ulab as np\n",
    "\n",
    "def timeit(n=1000):\n",
    "    def wrapper(f, *args, **kwargs):\n",
    "        func_name = str(f).split(' ')[1]\n",
    "        def new_func(*args, **kwargs):\n",
    "            run_times = np.zeros(n, dtype=np.uint16)\n",
    "            for i in range(n):\n",
    "                t = utime.ticks_us()\n",
    "                result = f(*args, **kwargs)\n",
    "                run_times[i] = utime.ticks_diff(utime.ticks_us(), t)\n",
    "            print('{}() execution times based on {} cycles'.format(func_name, n, (delta2-delta1)/n))\n",
    "            print('\\tbest: %d us'%np.min(run_times))\n",
    "            print('\\tworst: %d us'%np.max(run_times))\n",
    "            print('\\taverage: %d us'%np.mean(run_times))\n",
    "            print('\\tdeviation: +/-%.3f us'%np.std(run_times))            \n",
    "            return result\n",
    "        return new_func\n",
    "    return wrapper\n",
    "\n",
    "def timeit(f, *args, **kwargs):\n",
    "    func_name = str(f).split(' ')[1]\n",
    "    def new_func(*args, **kwargs):\n",
    "        t = utime.ticks_us()\n",
    "        result = f(*args, **kwargs)\n",
    "        print('execution time: ', utime.ticks_diff(utime.ticks_us(), t), ' us')\n",
    "        return result\n",
    "    return new_func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__END_OF_DEFS__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximation methods\n",
    "\n",
    "`ulab` implements five functions that can be used for interpolating, root finding, and minimising arbitrary `python` functions in one dimension. Two of these functions, namely, `interp`, and `trapz` are defined in `numpy`, while the other three are parts of `scipy`'s `optimize` module. \n",
    "\n",
    "Note that routines that work with user-defined functions still have to call the underlying `python` code, and therefore, gains in speed are not as significant as with other vectorised operations. As a rule of thumb, a factor of two can be expected, when compared to an optimised `python` implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## interp\n",
    "\n",
    "`numpy`: https://docs.scipy.org/doc/numpy/numpy.interp\n",
    "\n",
    "The `interp` function returns the linearly interpolated values of a one-dimensional numerical array. It requires three positional arguments,`x`, at which the interpolated values are evaluated, `xp`, the array\n",
    "of the independent data variable, and `fp`, the array of the dependent values of the data. `xp` must be a monotonically increasing sequence of numbers.\n",
    "\n",
    "Two keyword arguments, `left`, and `right` can also be supplied; these determine the return values, if `x < xp[0]`, and `x > xp[-1]`, respectively. If these arguments are not supplied, `left`, and `right` default to `fp[0]`, and `fp[-1]`, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-08T12:54:58.895801Z",
     "start_time": "2021-01-08T12:54:58.869338Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([0.8, 1.8, 2.8, 3.8, 4.8], dtype=float64)\n",
      "array([1.0, 1.8, 2.8, 4.6, 5.0], dtype=float64)\n",
      "array([0.0, 1.8, 2.8, 4.6, 5.0], dtype=float64)\n",
      "array([1.0, 1.8, 2.8, 4.6, 10.0], dtype=float64)\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%micropython -unix 1\n",
    "\n",
    "from ulab import numpy as np\n",
    "\n",
    "x = np.array([1, 2, 3, 4, 5]) - 0.2\n",
    "xp = np.array([1, 2, 3, 4])\n",
    "fp = np.array([1, 2, 3, 5])\n",
    "\n",
    "print(x)\n",
    "print(np.interp(x, xp, fp))\n",
    "print(np.interp(x, xp, fp, left=0.0))\n",
    "print(np.interp(x, xp, fp, right=10.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## newton\n",
    "\n",
    "`scipy`:https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton.html\n",
    "\n",
    "`newton` finds a zero of a real, user-defined function using the Newton-Raphson (or secant or Halley’s) method. The routine requires two positional arguments, the function, and the initial value. Three keyword\n",
    "arguments can be supplied to control the iteration. These are the absolute and relative tolerances `tol`, and `rtol`, respectively, and the number of iterations before stopping, `maxiter`. The function retuns a single scalar, the position of the root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-08T12:56:35.139958Z",
     "start_time": "2021-01-08T12:56:35.119712Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.260135727246117\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%micropython -unix 1\n",
    "\n",
    "from ulab import scipy as spy\n",
    "    \n",
    "def f(x):\n",
    "    return x*x*x - 2.0\n",
    "\n",
    "print(spy.optimize.newton(f, 3., tol=0.001, rtol=0.01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## bisect \n",
    "\n",
    "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.bisect.html\n",
    "\n",
    "`bisect` finds the root of a function of one variable using a simple bisection routine. It takes three positional arguments, the function itself, and two starting points. The function must have opposite signs\n",
    "at the starting points. Returned is the position of the root.\n",
    "\n",
    "Two keyword arguments, `xtol`, and `maxiter` can be supplied to control the accuracy, and the number of bisections, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-08T12:58:28.444300Z",
     "start_time": "2021-01-08T12:58:28.421989Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999997615814209\n",
      "only 8 bisections:  0.984375\n",
      "with 0.1 accuracy:  0.9375\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%micropython -unix 1\n",
    "\n",
    "from ulab import scipy as spy\n",
    "    \n",
    "def f(x):\n",
    "    return x*x - 1\n",
    "\n",
    "print(spy.optimize.bisect(f, 0, 4))\n",
    "\n",
    "print('only 8 bisections: ',  spy.optimize.bisect(f, 0, 4, maxiter=8))\n",
    "\n",
    "print('with 0.1 accuracy: ',  spy.optimize.bisect(f, 0, 4, xtol=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performance\n",
    "\n",
    "Since the `bisect` routine calls user-defined `python` functions, the speed gain is only about a factor of two, if compared to a purely `python` implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-19T19:08:24.750562Z",
     "start_time": "2020-05-19T19:08:24.682959Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bisect running in python\r\n",
      "execution time:  1270  us\r\n",
      "bisect running in C\r\n",
      "execution time:  642  us\r\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%micropython -pyboard 1\n",
    "\n",
    "from ulab import scipy as spy\n",
    "\n",
    "def f(x):\n",
    "    return (x-1)*(x-1) - 2.0\n",
    "\n",
    "def bisect(f, a, b, xtol=2.4e-7, maxiter=100):\n",
    "    if f(a) * f(b) > 0:\n",
    "        raise ValueError\n",
    "\n",
    "    rtb = a if f(a) < 0.0 else b\n",
    "    dx = b - a if f(a) < 0.0 else a - b\n",
    "    for i in range(maxiter):\n",
    "        dx *= 0.5\n",
    "        x_mid = rtb + dx\n",
    "        mid_value = f(x_mid)\n",
    "        if mid_value < 0:\n",
    "            rtb = x_mid\n",
    "        if abs(dx) < xtol:\n",
    "            break\n",
    "\n",
    "    return rtb\n",
    "\n",
    "@timeit\n",
    "def bisect_scipy(f, a, b):\n",
    "    return spy.optimize.bisect(f, a, b)\n",
    "\n",
    "@timeit\n",
    "def bisect_timed(f, a, b):\n",
    "    return bisect(f, a, b)\n",
    "\n",
    "print('bisect running in python')\n",
    "bisect_timed(f, 3, 2)\n",
    "\n",
    "print('bisect running in C')\n",
    "bisect_scipy(f, 3, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fmin\n",
    "\n",
    "`scipy`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin.html\n",
    "\n",
    "The `fmin` function finds the position of the minimum of a user-defined function by using the downhill simplex method. Requires two positional arguments, the function, and the initial value. Three keyword arguments, `xatol`, `fatol`, and `maxiter` stipulate conditions for stopping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-08T13:00:26.729947Z",
     "start_time": "2021-01-08T13:00:26.702748Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9996093749999952\n",
      "1.199999999999996\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%micropython -unix 1\n",
    "\n",
    "from ulab import scipy as spy\n",
    "\n",
    "def f(x):\n",
    "    return (x-1)**2 - 1\n",
    "\n",
    "print(spy.optimize.fmin(f, 3.0))\n",
    "print(spy.optimize.fmin(f, 3.0, xatol=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## trapz\n",
    "\n",
    "`numpy`: https://numpy.org/doc/stable/reference/generated/numpy.trapz.html\n",
    "\n",
    "The function takes one or two one-dimensional `ndarray`s, and integrates the dependent values (`y`) using the trapezoidal rule. If the independent variable (`x`) is given, that is taken as the sample points corresponding to `y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-08T13:01:29.515166Z",
     "start_time": "2021-01-08T13:01:29.494285Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x:  array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0], dtype=float64)\n",
      "y:  array([0.0, 1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0], dtype=float64)\n",
      "============================\n",
      "integral of y:  244.5\n",
      "integral of y at x:  244.5\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%micropython -unix 1\n",
    "\n",
    "from ulab import numpy as np\n",
    "\n",
    "x = np.linspace(0, 9, num=10)\n",
    "y = x*x\n",
    "\n",
    "print('x: ',  x)\n",
    "print('y: ',  y)\n",
    "print('============================')\n",
    "print('integral of y: ', np.trapz(y))\n",
    "print('integral of y at x: ', np.trapz(y, x=x))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "382.797px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}