{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Effective computation in Biomechanics\n", "\n", "Romain Martinez " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [Python fundamentals](01.01-python-base.ipynb) | [Contents](index.ipynb) | [Biomechanical analysis with Pyomeca](01.03-intro-to-pyomeca.ipynb) >" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Scientific computing with Numpy\n", "\n", "\n", "\n", "[Numpy](https://numpy.org/) is probably the most fundamental package for scientific computing in Python.\n", "\n", "It provides a highly efficient interface to create and interact with multidimensional arrays.\n", "\n", "Nearly every other scientific package uses or integrates with NumPy in some way.\n", "\n", "The performance-sensitive parts of NumPy are all written in the C language, so they are very fast." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numpy arrays" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fundamental object of NumPy is its `numpy.array`, an $n$-dimensional array that is also present in some form in array-oriented languages such as MATLAB." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np # we import numpy as \"np\" to avoid typing \"numpy\" each time\n", "\n", "%load_ext lab_black\n", "\n", "vector = np.array([10, 11, 12])" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "print(type(vector))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(3,)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vector.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let’s start things off by forming a 3-dimensional array with 36 consecutive numbers" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[[ 0, 1, 2],\n", " [ 3, 4, 5],\n", " [ 6, 7, 8],\n", " [ 9, 10, 11]],\n", "\n", " [[12, 13, 14],\n", " [15, 16, 17],\n", " [18, 19, 20],\n", " [21, 22, 23]],\n", "\n", " [[24, 25, 26],\n", " [27, 28, 29],\n", " [30, 31, 32],\n", " [33, 34, 35]]])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "matrix = np.arange(36).reshape(3, 4, 3)\n", "matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Picturing high-dimensional arrays in two dimensions can be difficult.\n", "\n", "One intuitive way to think about an array’s shape is to simply say:\n", "\n", "> matrix is a 3 by 4 by 3 array" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(3, 4, 3)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "matrix.shape" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "matrix.ndim=3\n", "matrix.shape=(3, 4, 3)\n", "matrix.size=36\n", "matrix.dtype=dtype('int64')\n" ] } ], "source": [ "print(f\"{matrix.ndim=}\")\n", "print(f\"{matrix.shape=}\")\n", "print(f\"{matrix.size=}\")\n", "print(f\"{matrix.dtype=}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visually, `matrix` could be thought of as a container of three 4x3 grids (or a rectangular prism) and would look like this:\n", "\n", "\n", "\n", "Data with greater than two dimensions are tough to picture, but they are __everywhere__ in biomechanics, so better get used to it!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating a 4x3 array" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0, 1, 2],\n", " [ 3, 4, 5],\n", " [ 6, 7, 8],\n", " [ 9, 10, 11]])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10, 11]]) # from a list" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0, 1, 2],\n", " [ 3, 4, 5],\n", " [ 6, 7, 8],\n", " [ 9, 10, 11]])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.arange(12).reshape(4, 3) # with consecutive numbers" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0., 0., 0.],\n", " [0., 0., 0.],\n", " [0., 0., 0.],\n", " [0., 0., 0.]])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.zeros((4, 3)) # with zeros" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1., 1., 1.],\n", " [1., 1., 1.],\n", " [1., 1., 1.],\n", " [1., 1., 1.]])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.ones((4, 3)) # with ones" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1., 1., 1.],\n", " [1., 1., 1.],\n", " [1., 1., 1.],\n", " [1., 1., 1.]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.empty((4, 3)) # with last value used" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.55456578, 0.29574095, 0.00245861],\n", " [0.40599595, 0.56963864, 0.10437213],\n", " [0.83959092, 0.32007413, 0.0839273 ],\n", " [0.46797033, 0.07439805, 0.56752482]])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.random.rand(4, 3) # with random numbers" ] }, { "cell_type": "markdown", "metadata": { "toc-hr-collapsed": false }, "source": [ "### Array indexing\n", "\n", "More details in the [indexing](https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html) section of the documentation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Integer indexing" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0, 1, 2],\n", " [ 3, 4, 5],\n", " [ 6, 7, 8],\n", " [ 9, 10, 11]])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "array = np.arange(12).reshape(4, 3)\n", "array" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 5])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# first row second column & second row third column\n", "array[[0, 1], [1, 2]]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Slicing" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0, 1, 2],\n", " [3, 4, 5]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# first two rows\n", "array[:2, :]" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([10, 11])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# last row and columns 1 & 2\n", "array[-1, 1:3]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Warning**: a slice of an array is a view into the same data, so modifying it will **modify the original array**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Boolean indexing" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[False, False, False],\n", " [False, False, False],\n", " [ True, True, True],\n", " [ True, True, True]])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "array > 5" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 6, 7, 8, 9, 10, 11])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "array[array > 5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Array math\n", "\n", "More details in the [mathematical functions](https://docs.scipy.org/doc/numpy/reference/routines.math.html) section of the official documentation." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "x = np.array([[1, 2], [3, 4]])\n", "y = np.array([[5, 6], [7, 8]])" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 6, 8],\n", " [10, 12]])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x + y\n", "# or np.add(x, y)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-4, -4],\n", " [-4, -4]])" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x - y\n", "# or np.substract(x, y)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 5, 12],\n", " [21, 32]])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x * y\n", "# or np.multiply(x, y)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.2 , 0.33333333],\n", " [0.42857143, 0.5 ]])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x / y\n", "# or np.divide(x, y)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1. , 1.41421356],\n", " [1.73205081, 2. ]])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.sqrt(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Warning**: unlike MATLAB, `*` is elementwise multiplication, not matrix multiplication.\n", "\n", "We instead use the dot function to compute inner products of vectors, to multiply a vector by a matrix, and to multiply matrices." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[19, 22],\n", " [43, 50]])" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x @ y\n", "# or np.dot(x, y)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([4, 6])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.sum(x, axis=0)\n", "# or x.sum(axis=0)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([3, 7])" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.sum(x, axis=1)\n", "# or x.sum(axis=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reshaping" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 2],\n", " [3, 4]])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 3],\n", " [2, 4]])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.T" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 2, 3, 4])" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.reshape(-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Your turn\n", "\n", "Here are some EMG data from the anterior deltoid during a box lifting task:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "" ], "text/plain": [ "alt.LayerChart(...)" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "import altair as alt\n", "\n", "emg_dataframe = pd.read_csv(\"../data/emgs.csv\").query('filename == \"12kg_H2_3\"')\n", "\n", "base = alt.Chart(emg_dataframe)\n", "line = base.mark_line().encode(x=\"index\", y=\"delt_ant\")\n", "rule = base.mark_rule(color=\"firebrick\", size=2).encode(y=\"mean(delt_ant)\")\n", "\n", "above_base = alt.Chart(emg_dataframe.query(\"delt_ant > delt_ant.mean()\"))\n", "points = above_base.mark_point(color=\"red\").encode(x=\"index\", y=\"delt_ant\")\n", "rect = above_base.mark_errorband(extent=\"stdev\").encode(y=\"delt_ant\")\n", "\n", "line + points + rule + rect" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "emg = emg_dataframe[\"delt_ant\"].values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the `emg` Numpy array:\n", "1. Selects only points that are above the average of the EMG signal (red points).\n", "2. From these points, calculate how many points are within plus or minus one standard deviation of the mean (red points outside the blue area)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [Python fundamentals](01.01-python-base.ipynb) | [Contents](index.ipynb) | [Biomechanical analysis with Pyomeca](01.03-intro-to-pyomeca.ipynb) >" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "conda-env-tutorials-py" }, "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.1" } }, "nbformat": 4, "nbformat_minor": 4 }