diff --git a/src/test/data/A Reaction-Diffusion Equation Solver in Python with Numpy.ipynb b/src/test/data/A Reaction-Diffusion Equation Solver in Python with Numpy.ipynb
new file mode 100644
index 00000000..d91fd725
--- /dev/null
+++ b/src/test/data/A Reaction-Diffusion Equation Solver in Python with Numpy.ipynb
@@ -0,0 +1,785 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This notebook demonstrates how IPython notebooks can be used to discuss the theory and implementation of numerical algorithms on one page.\n",
+ "\n",
+ "With `ipython nbconvert --to markdown name.ipynb` a notebook like this one can be made into a \n",
+ "[blog post](http://georg.io/2013/12/Crank_Nicolson) in one easy step. To display the graphics in your resultant blog post use,\n",
+ "for instance, your [Dropbox Public folder](https://www.dropbox.com/help/16/en) that you can \n",
+ "where $N$ and $J$ are the number of discrete time and space points in our grid respectively.\n",
+ "$\\Delta t$ and $\\Delta x$ are the time step and space step respectively and defined as follows:\n",
+ "\n",
+ "$$\\Delta t = T / N,$$\n",
+ "\n",
+ "$$\\Delta x = L / J,$$\n",
+ "\n",
+ "where $T$ is the point in time up to which we will integrate $u$ numerically.\n",
+ "\n",
+ "Our ultimate goal is to construct a numerical method that allows us to approximate the unknonwn analytic solution $u(x,t)$\n",
+ "reasonably well in these discrete grid points.\n",
+ "\n",
+ "That is we want construct a method that computes values $U(j \\Delta x, n \\Delta t)$ (note: capital $U$) so that\n",
+ "\n",
+ "$$U(j \\Delta x, n \\Delta t) \\approx u(j \\Delta x, n \\Delta t)$$\n",
+ "\n",
+ "As a shorthand we will write $U_j^n = U(j \\Delta x, n \\Delta t)$ and $(j,n)$ to refer to grid point $(j \\Delta x, n \\Delta t)$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## The Crank-Nicolson Stencil"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Based on the two-dimensional grid we construct we then approximate the operators of our reaction-diffusion system.\n",
+ "\n",
+ "For instance, to approximate the time derivative on the left-hand side in grid point $(j,n)$ we use the values of $U$ in two specific grid points:\n",
+ "\n",
+ "$$\\frac{\\partial u}{\\partial t}\\Bigg|_{x = j \\Delta x, t = n \\Delta t} \\approx \\frac{U_j^{n+1} - U_j^n}{\\Delta t}.$$\n",
+ "\n",
+ "We can think of this scheme as a stencil that we superimpose on our $(x,t)$-grid and this particular stencil is\n",
+ "commonly referred to as [forward difference](http://en.wikipedia.org/wiki/Finite_difference#Forward.2C_backward.2C_and_central_differences).\n",
+ "\n",
+ "The spatial part of the [Crank-Nicolson stencil](http://journals.cambridge.org/abstract_S0305004100023197)\n",
+ "(or see [Table 3.2.2 of Trefethen](http://people.maths.ox.ac.uk/trefethen/3all.pdf))\n",
+ "for the heat equation ($u_t = u_{xx}$) approximates the \n",
+ "[Laplace operator](http://en.wikipedia.org/wiki/Laplace_operator) of our equation and takes the following form\n",
+ "To approximate $f(u(j \\Delta x, n \\Delta t))$ we write simply $f(U_j^n)$.\n",
+ "\n",
+ "These approximations define the stencil for our numerical method as pictured on [Wikipedia](http://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method).\n",
+ "Note that since our numerical integration starts with a well-defined initial condition at $n=0$, $\\mathbf{U}^0$, the\n",
+ "vector $\\mathbf{U}^{n+1}$ on the left-hand side is the only unknown in this system of linear equations.\n",
+ "\n",
+ "Thus, to integrate numerically our reaction-diffusion system from time point $n$ to $n+1$ we need to solve numerically for vector $\\mathbf{U}^{n+1}$.\n",
+ "\n",
+ "Let us call the matrix on the left-hand side $A$, the one on the right-hand side $B$,\n",
+ "and the vector on the right-hand side $\\mathbf{f}^n$.\n",
+ "Using this notation we can write the above system as\n",
+ "\n",
+ "$$A \\mathbf{U}^{n+1} = B \\mathbf{U}^n + f^n.$$\n",
+ "\n",
+ "In this linear equation, matrices $A$ and $B$ are defined by our problem: we need to specify these matrices once for our\n",
+ "problem and incorporate our boundary conditions in them.\n",
+ "Vector $\\mathbf{f}^n$ is a function of $\\mathbf{U}^n$ and so needs to be reevaluated in every time point $n$.\n",
+ "We also need to carry out one matrix-vector multiplication every time point, $B \\mathbf{U}^n$, and\n",
+ "Our one-dimensional domain has unit length and we define `J = 100` equally spaced\n",
+ "grid points in this domain.\n",
+ "This divides our domain into `J-1` subintervals, each of length `dx`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "L = 1.\n",
+ "J = 100\n",
+ "dx = float(L)/float(J-1)\n",
+ "x_grid = numpy.array([j*dx for j in range(J)])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Equally, we define `N = 1000` equally spaced grid points on our time domain of length `T = 200` thus dividing our time domain into `N-1` intervals of length `dt`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "T = 200\n",
+ "N = 1000\n",
+ "dt = float(T)/float(N-1)\n",
+ "t_grid = numpy.array([n*dt for n in range(N)])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Specify System Parameters and the Reaction Term"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We choose our parameter values based on the work by\n",
+ "A_u = numpy.diagflat([-sigma_u for i in range(J-1)], -1) +\\\n",
+ " numpy.diagflat([1.+sigma_u]+[1.+2.*sigma_u for i in range(J-2)]+[1.+sigma_u]) +\\\n",
+ " numpy.diagflat([-sigma_u for i in range(J-1)], 1)\n",
+ " \n",
+ "B_u = numpy.diagflat([sigma_u for i in range(J-1)], -1) +\\\n",
+ " numpy.diagflat([1.-sigma_u]+[1.-2.*sigma_u for i in range(J-2)]+[1.-sigma_u]) +\\\n",
+ " numpy.diagflat([sigma_u for i in range(J-1)], 1)\n",
+ " \n",
+ "A_v = numpy.diagflat([-sigma_v for i in range(J-1)], -1) +\\\n",
+ " numpy.diagflat([1.+sigma_v]+[1.+2.*sigma_v for i in range(J-2)]+[1.+sigma_v]) +\\\n",
+ " numpy.diagflat([-sigma_v for i in range(J-1)], 1)\n",
+ " \n",
+ "B_v = numpy.diagflat([sigma_v for i in range(J-1)], -1) +\\\n",
+ " numpy.diagflat([1.-sigma_v]+[1.-2.*sigma_v for i in range(J-2)]+[1.-sigma_v]) +\\\n",
+ " numpy.diagflat([sigma_v for i in range(J-1)], 1)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To confirm, this is what `A_u` looks like:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[ 1.981 -0.981 0. ... 0. 0. 0. ]\n",
+ " [-0.981 2.962 -0.981 ... 0. 0. 0. ]\n",
+ " [ 0. -0.981 2.962 ... 0. 0. 0. ]\n",
+ " ...\n",
+ " [ 0. 0. 0. ... 2.962 -0.981 0. ]\n",
+ " [ 0. 0. 0. ... -0.981 2.962 -0.981]\n",
+ " [ 0. 0. 0. ... 0. -0.981 1.981]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "print A_u"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Solve the System Iteratively"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To advance our system by one time step, we need to do one matrix-vector multiplication followed by one vector-vector addition on the right hand side.\n",
+ "\n",
+ "To facilitate this, we rewrite our reaction term so that it accepts concentration vectors $\\mathbf{U}^n$ and $\\mathbf{V}^n$ as arguments\n",
+ "and returns vector $\\mathbf{f}^n$.\n",
+ "\n",
+ "As a reminder, this is our non-vectorial definition of $f$\n",
+ "\n",
+ " f = lambda u, v: v*(k0 + float(u*u)/float(1. + u*u)) - u"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "f_vec = lambda U, V: numpy.multiply(dt, numpy.subtract(numpy.multiply(V, \n",
+ "- using the online [IPython notebook slide viewer](https://slideviewer.herokuapp.com/) (some slides of the notebook might not be properly rendered).\n",
+ "\n",
+ "This and other related IPython notebooks can be found at the course github repository:\n",
+ "* [Random variable](http://en.wikipedia.org/wiki/Random_variable): a variable whose value is subject to variations due to __chance__. A random variable can take on a set of possible different values, each with an associated probability, in contrast to other mathematical variables.\n",
+ "\n",
+ "* [Probability distribution](http://en.wikipedia.org/wiki/Probability_distribution): mathematical function describing the possible values of a random variable and their associated probabilities.\n",
+ "\n",
+ "* [Probability density function (pdf)](http://en.wikipedia.org/wiki/Probability_density_function) of a __continuous random variable__ is a function that describes the relative likelihood for this random variable to take on a given value. \n",
+ " * The probability of the random variable falling within a particular range of values is given by the integral of this variable’s density over that range.\n",
+ " * The probability density function is nonnegative everywhere, and its integral over the entire space is equal to one.\n",
+ "The probability distribution of a random variable is often characterised by a small number of parameters, which also have a practical interpretation.\n",
+ "\n",
+ "* [Mean](http://en.wikipedia.org/wiki/Mean) (a.k.a expected value) refers to one measure of the central tendency either of a probability distribution or of the random variable characterized by that distribution.\n",
+ " * population mean: $\\mu = \\operatorname{E}[X]$.\n",
+ " * estimation of sample mean: $\\bar{x}$.\n",
+ "* [Standard deviation](http://en.wikipedia.org/wiki/Standard_deviation) measures the amount of variation or dispersion from the mean.\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
+ " warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
+ "In applied mathematics, [test functions](http://en.wikipedia.org/wiki/Test_functions_for_optimization), also known as artificial landscapes, are useful to evaluate characteristics of optimization algorithms, such as:\n",
+ "\n",
+ "* Velocity of convergence.\n",
+ "* Precision.\n",
+ "* Robustness.\n",
+ "* General performance.\n",
+ "\n",
+ "DEAP has a number of test problems already implemented. See http://deap.readthedocs.org/en/latest/api/benchmarks.html"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " s)\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ "* In an evolution strategy, new candidate solutions are sampled according to a multivariate normal distribution in the $\\mathbb{R}^n$. \n",
+ "* Recombination amounts to selecting a new mean value for the distribution. \n",
+ "* Mutation amounts to adding a random vector, a perturbation with zero mean. \n",
+ "* Pairwise dependencies between the variables in the distribution are represented by a covariance matrix. \n",
+ "\n",
+ "### The covariance matrix adaptation (CMA) is a method to update the covariance matrix of this distribution. \n",
+ "\n",
+ "> This is particularly useful, if the objective function $f()$ is ill-conditioned."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "internals": {
+ "slide_helper": "subslide_end",
+ "slide_type": "subslide"
+ },
+ "slide_helper": "slide_end",
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "### CMA-ES features\n",
+ "\n",
+ "* Adaptation of the covariance matrix amounts to learning a second order model of the underlying objective function.\n",
+ "* This is similar to the approximation of the inverse Hessian matrix in the Quasi-Newton method in classical optimization. \n",
+ "* In contrast to most classical methods, fewer assumptions on the nature of the underlying objective function are made. \n",
+ "* *Only the ranking between candidate solutions is exploited* for learning the sample distribution and neither derivatives nor even the function values themselves are required by the method."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 46,
+ "metadata": {
+ "internals": {
+ "slide_type": "subslide"
+ },
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from deap import cma"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "internals": {},
+ "slideshow": {
+ "slide_type": "-"
+ }
+ },
+ "source": [
+ "A similar setup to the previous one."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 47,
+ "metadata": {
+ "internals": {
+ "slide_helper": "subslide_end"
+ },
+ "slide_helper": "slide_end",
+ "slideshow": {
+ "slide_type": "-"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/home/mmmm1998/.local/lib/python3.6/site-packages/deap/creator.py:141: RuntimeWarning: A class named 'Individual' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.\n",
+ "### OK, but wouldn't it be nice to have an animated plot of how CMA-ES progressed? \n",
+ "\n",
+ "* We need to do some coding to make this animation work.\n",
+ "* We are going to create a class named `PlotableStrategy` that inherits from `deap.cma.Strategy`. This class logs the features we need to make the plots as evolution takes place. That is, for every iteration we store:\n",
+ " * Current centroid and covariance ellipsoid.\n",
+ " * Updated centroid and covariance.\n",
+ " * Sampled individuals.\n",
+ " * Evolution path.\n",
+ " \n",
+ "_Note_: I think that DEAP's implementation of CMA-ES has the drawback of storing information that should be stored as part of \"individuals\". I leave this for an afternoon hack."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 50,
+ "metadata": {
+ "internals": {
+ "slide_helper": "subslide_end",
+ "slide_type": "subslide"
+ },
+ "slide_helper": "slide_end",
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from math import sqrt, log, exp\n",
+ "class PlotableStrategy(cma.Strategy):\n",
+ " \"\"\"This is a modification of deap.cma.Strategy class.\n",
+ " We store the execution data in order to plot it.\n",
+ " **Note:** This class should not be used for other uses than\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n",
+ " This work is licensed under a [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License](http://creativecommons.org/licenses/by-nc-sa/4.0/).\n",
diff --git a/src/test/data/Automata and Computability using Jupyter.ipynb b/src/test/data/Automata and Computability using Jupyter.ipynb
new file mode 100644
index 00000000..145a82bd
--- /dev/null
+++ b/src/test/data/Automata and Computability using Jupyter.ipynb
@@ -0,0 +1,2095 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ }
+ },
+ "source": [
+ "# Jove helps teach models of computation using Jupyter \n",
+ "\n",
+ "Included are modules on:\n",
+ "\n",
+ "* Sets, strings and languages\n",
+ "* Language operations\n",
+ "* Construction of and operations on DFA and NFA\n",
+ "* Regular expression parsing and automata inter-conversion\n",
+ "* Derivate-based parsing\n",
+ "* Pushdown automata\n",
+ "* The construction of parsers using context-free productions, including\n",
+ " a full lexer/parser for Jove's own markdown syntax\n",
+ "* Studies of parsing: ambiguity, associativity, precedence\n",
+ "* Turing machines (including one for the Collatz problem)\n",
+ "\n",
+ "For a complete Jove top-level reference, kindly refer to https://github.com/ganeshutah/Jove from where you can download and obtain Jove. You can also visit this Github link now and poke around (the NBViewer will display the contents).\n",
+ "\n",
+ "Once you are in the top-level Gallery link we provide, feel free to explore the hierarchy of modules found there.\n",
+ "\n",
+ "These notebooks should give you an idea of the contents.\n",
+ "\n",
+ "* [DFA Illustrations (has a Youtube)](http://nbviewer.jupyter.org/github/ganeshutah/Jove/blob/master/notebooks/tutorial/DFAUnit2.ipynb)\n",
+ "<IPython.lib.display.YouTubeVideo at 0x7fa7a1ee4c50>"
+ ]
+ },
+ "execution_count": 1,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "from IPython.display import YouTubeVideo\n",
+ "YouTubeVideo('dGcLHtYLgDU')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "You may use any of these help commands:\n",
+ "help(ResetStNum)\n",
+ "help(NxtStateStr)\n",
+ "\n",
+ "You may use any of these help commands:\n",
+ "help(mkp_dfa)\n",
+ "help(mk_dfa)\n",
+ "help(totalize_dfa)\n",
+ "help(addtosigma_delta)\n",
+ "help(step_dfa)\n",
+ "help(run_dfa)\n",
+ "help(accepts_dfa)\n",
+ "help(comp_dfa)\n",
+ "help(union_dfa)\n",
+ "help(intersect_dfa)\n",
+ "help(pruneUnreach)\n",
+ "help(iso_dfa)\n",
+ "help(langeq_dfa)\n",
+ "help(same_status)\n",
+ "help(h_langeq_dfa)\n",
+ "help(fixptDist)\n",
+ "help(min_dfa)\n",
+ "help(pairFR)\n",
+ "help(state_combos)\n",
+ "help(sepFinNonFin)\n",
+ "help(bash_eql_classes)\n",
+ "help(listminus)\n",
+ "help(bash_1)\n",
+ "help(mk_rep_eqc)\n",
+ "help(F_of)\n",
+ "help(rep_of_s)\n",
+ "help(q0_of)\n",
+ "help(Delta_of)\n",
+ "help(mk_state_eqc_name)\n",
+ "\n",
+ "You may use any of these help commands:\n",
+ "help(mk_nfa)\n",
+ "help(totalize_nfa)\n",
+ "help(step_nfa)\n",
+ "help(run_nfa)\n",
+ "help(ec_step_nfa)\n",
+ "help(Eclosure)\n",
+ "help(Echelp)\n",
+ "help(accepts_nfa)\n",
+ "help(nfa2dfa)\n",
+ "help(n2d)\n",
+ "help(inSets)\n",
+ "help(rev_dfa)\n",
+ "help(min_dfa_brz)\n",
+ "\n",
+ "You may use any of these help commands:\n",
+ "help(re2nfa)\n",
+ "\n",
+ "You may use any of these help commands:\n",
+ "help(RE2Str)\n",
+ "help(mk_gnfa)\n",
+ "help(mk_gnfa_from_D)\n",
+ "help(dfa2nfa)\n",
+ "help(del_gnfa_states)\n",
+ "help(gnfa_w_REStr)\n",
+ "help(del_one_gnfa_state)\n",
+ "help(Edges_Exist_Via)\n",
+ "help(choose_state_to_del)\n",
+ "help(form_alt_RE)\n",
+ "help(form_concat_RE)\n",
+ "help(form_kleene_RE)\n",
+ "\n",
+ "You may use any of these help commands:\n",
+ "help(md2mc)\n",
+ ".. and if you want to dig more, then ..\n",
+ "help(default_line_attr)\n",
+ "help(length_ok_input_items)\n",
+ "help(union_line_attr_list_fld)\n",
+ "help(extend_rsltdict)\n",
+ "help(form_delta)\n",
+ "help(get_machine_components)\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "import sys\n",
+ "sys.path[0:0] = ['/home/mmmm1998/Документы/Репозитории/Jove','/home/mmmm1998/Документы/Репозитории/Jove/3rdparty'] # Put these at the head of the search path\n",
+ "from jove.DotBashers import *\n",
+ "from jove.Def_DFA import *\n",
+ "from jove.Def_NFA import *\n",
+ "from jove.Def_RE2NFA import *\n",
+ "from jove.Def_NFA2RE import *\n",
+ "from jove.Def_md2mc import *"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ }
+ },
+ "source": [
+ " # Jove allows you to set problems in markdown and have students solve"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ }
+ },
+ "source": [
+ "1) LOdd1Three0 : Set of strings over {0,1} with an odd # of 1s OR exactly three 0s. \n",
+ "\n",
+ "* Hint on how to arrive at the language:\n",
+ "\n",
+ " - develop NFAs for the two cases and perform their union. Obtain DFA\n",
+ "\n",
+ " - develop REs for the two cases and perform the union. \n",
+ "\n",
+ " - Testing the creations:\n",
+ "\n",
+ " . Come up with language for even # of 1s and separately for \"other than three 0s\". \n",
+ " \n",
+ " . Do two intersections. \n",
+ " \n",
+ " . Is the language empty?\n",
+ "\n",
+ "\n",
+ "2) Language of strings over {0,1} with exactly two occurrences of 0101 in it.\n",
+ "\n",
+ " * Come up with it directly (take overlaps into account, i.e. 010101 has two occurrences in it\n",
+ "\n",
+ " * Come up in another way\n",
+ "\n",
+ "Notes:\n",
+ "\n",
+ "* Most of the problem students will have in this course is interpreting English (technical English)\n",
+ "\n",
+ "* So again, read the writeup at the beginning of Module6 (should be ready soon today) and work on using the tool.\n",
+ "\n",
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ }
+ },
+ "source": [
+ "__Solutions__\n",
+ "\n",
+ "1) LOdd1Three0 : Set of strings over {0,1} with an odd # of 1s OR exactly three 0s. \n",
+ "\n",
+ "* Hint on how to arrive at the language:\n",
+ "\n",
+ " - develop NFAs for the two cases and perform their union. Obtain DFA\n",
+ "\n",
+ " - develop REs for the two cases and perform the union. \n",
+ "\n",
+ " - Testing the creations:\n",
+ "\n",
+ " . Come up with language for even # of 1s and separately for \"other than three 0s\". \n",
+ " \n",
+ " . Do two intersections. \n",
+ " \n",
+ " . Is the language empty?\n",
+ "\n",
+ "\n",
+ "2) Language of strings over {0,1} with exactly two occurrences of 0101 in it.\n",
+ "\n",
+ " * Come up with it directly (take overlaps into account, i.e. 010101 has two occurrences in it\n",
+ "\n",
+ " * Come up in another way\n",
+ "\n",
+ "Notes:\n",
+ "\n",
+ "* Most of the problem students will have in this course is interpreting English (technical English)\n",
+ "\n",
+ "* So again, read the writeup at the beginning of Module6 (should be ready soon today) and work on using the tool.\n",
+ "<p>We live in a complex environment and must constantly integrate sensory information to interact with the world around us. Inputs from different modalities might not always be congruent with each other, but dissociating the true nature of the stimulus may be a matter of life or death for an organism.</p>\n",
+ "<p>You hear and see evidence of a rattlesnake in tall grass near you. You get an auditory and a visual cue of the snake's location $s$. Both cues are associated with a likelihood function indicating the probability of that cue for all possible locations of the snake. The likelihood function associated with the visual cue, $p(c_\\mathrm{V}|s)$, has high uncertainty, because of the tall grass. The auditory cue is easier to localize, so its associated likelihood function, $p(c_\\mathrm{A}|s)$, is sharper. In accordance with Bayes' Rule, and assuming a flat prior over the snake's location, an optimal estimate of the location of the snake can be computed by multiplying the two likelihoods. This joint likelihood will be between the two cues but closer to the less uncertain cue, and will have less uncertainty than both unimodal likelihood functions.</p>"
+ "Behavioral experiments have demonstrated that humans perform near-optimal Bayesian inference on ambiguous sensory information (van Beers *et al.*, 1999; Ernst & Banks, 2002; Kording & Wolpert, 2004; Stocker & Simoncelli, 2006). This has been demonstrated in cue combination experiments in which subjects report a near-optimal estimate of the stimulus given two noisy measurements of that stimulus. However, the neural basis for how humans might perform these computations is unclear. \n",
+ "\n",
+ "Ma *et. al.* (2006) propose that variance in cortical activity, rather than impairing sensory systems, is an adaptive mechanism to encode uncertainty in sensory measurements. They provide theory showing how the brain might use probabilistic population codes to perform near-optimal cue combination. We will re-derive the theory in here, and demonstrate it by simulating and decoding neural populations.\n",
+ "\n",
+ "## Cues can be represented by neural populations\n",
+ "\n",
+ "To return to our deadly rattlesnake, let's now assume that $c_\\mathrm{V}$ and $c_\\mathrm{A}$ are represented by populations of neurons $\\mathbf{r}_\\mathrm{V}$ and $\\mathbf{r}_\\mathrm{A}$, respectively. For our math and simulations, we assume that $\\mathbf{r}_\\mathrm{V}$ and $\\mathbf{r}_\\mathrm{A}$ are each composed of $N$ neurons that:\n",
+ "\n",
+ "* have independent Poisson variability\n",
+ "* have regularly spaced Gaussian tuning curves that are identical in mean and variance for neurons with the same index in both populations\n",
+ "\n",
+ "The populations may have different gains, $g_\\mathrm{V}$ and $g_\\mathrm{A}$.\n",
+ "\n",
+ "These are the tuning curves for the neurons in $\\mathbf{r}_\\mathrm{V}$ (purple) and $\\mathbf{r}_\\mathrm{A}$ (red). Each curve represents the mean firing rate of a single neuron given a location $s$. Each neuron thus has a preferred location, which is where its tuning curve peaks."
+ "The tuning curves are dense enough that we can also assume that $\\sum_{i=0}^N f_i(s) = k$ (*i.e.*, the sum of the tuning curves in a population is constant.)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First, we will show how the brain can decode a likelihood over stimulus from neural activity. Then we will ask how the brain can compute joint likelihoods.\n",
+ "### How can the brain decode $p(\\mathbf{r_\\mathrm{V}}|s)$?\n",
+ "## How can the brain compute the joint likelihood $p(\\mathbf{r}_\\mathrm{V}|s)\\ p(\\mathbf{r}_\\mathrm{A}|s)$?\n",
+ "The fact that we see neurons from $\\mathbf{r}_\\mathrm{V}$ and $\\mathbf{r}_\\mathrm{A}$ being added on a neuron-by-neuron basis in the exponent above suggests that we could construct a third population vector, $\\mathbf{r}_\\mathrm{V}+\\mathbf{r}_\\mathrm{A}$, and decode that.\n",
+ "\n",
+ "### First, we must prove that the sum of two Poisson-distributed random variables $X+Y$ is again Poisson-distributed.\n",
+ "Since equations $(9)$ and $(27)$ are proportional, we have shown that optimal cue combination can be executed by decoding linear sums of populations."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$x = 2$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Simulation\n",
+ "Here are the spike counts (during 1 s) from the two populations on one trial. Depicted in blue is a third population vector that is the sum of $\\mathbf{r}_\\mathrm{V}$ and $\\mathbf{r}_\\mathrm{A}$."
+ "Here are the decoded likelihoods for each population alone $(6)$, the joint likelihood $(9)$, and the likelihood for the summed population $(27)$. Note that the joint likelihood (gold) is less uncertain than either unimodal likelihood. Also note that it is identical to the likelihood for the summed population (blue)."
+ "Here, we break the assumption that the two populations have the same tuning curve width. Note that the joint likelihood (gold) is no longer identical to the likelihood for the summed population (blue)."
+ "Now you can play interactively with the parameters of the simulation using these sliders, and watch the decoded likelihoods shift around. Every time you change a parameter, new sets of spikes are generated and used to infer $s$.\n",
+ "\n",
+ "For the simulation to be interactive, you'll have to download this notebook."
+ "i = ipywidgets.interactive(spikes_and_inference,\n",
+ " true_stimulus = (-40, 40, .1),\n",
+ " number_of_neurons = (2, 200, 1),\n",
+ " r_V_gain = (0, 100, 1),\n",
+ " r_A_gain = (0, 100, 1),\n",
+ " r_V_tuning_curve_sigma = (0.1, 50, .1),\n",
+ " r_A_tuning_curve_sigma = (0.1, 50, .1),\n",
+ " tuning_curve_baseline = (0, 20, .1));\n",
+ "display(ipywidgets.VBox(i.children[2:-1]))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Conclusion\n",
+ "\n",
+ "It has been shown behaviorally that humans perform near-optimal Bayesian inference on ambiguous sensory information. As suggested by Ma *et. al.* (2006) and shown here, it is possible that the brain does this operation by simply performing linear combinations of populations of Poisson neurons receiving various sensory input. Cortical neurons may be particularly well suited for this task because they have Poisson-like firing rates, displaying reliable variability from trial to trial (Tolhurst, Movshon & Dean, 1982; Softky & Koch, 1993).\n",
+ "\n",
+ "High levels of noise in these populations might at first be difficult to reconcile considering highly precise behavioral data. However, variability in neural populations might be direcly representative of uncertainty in environmental stimuli. Variability in cortical populations would then be critical for precise neural coding.\n",
+ "\n",
+ "## References\n",
+ "\n",
+ "* Ernst MO, Banks MS. (2002). Humans integrate visual and haptic information in a statistically optimal fashion. *Nature.*\n",
+ "* Ma WJ, Beck JM, Latham PE, Pouget A. (2006). Bayesian inference with probabilistic population codes. *Nature Neuroscience.*\n",
+ "* Softky WR, Koch C. (1993). The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. *Journal of Neuroscience.*\n",
+ "* Stocker AA, Simoncelli EP. (2006). Noise characteristics and prior expectations in human visual speed perception. *Nature Neuroscience.*\n",
+ "* Tolhurst, DJ, Movshon JA, Dean AF. (1983). The statistical reliability of signals in single neurons in cat and monkey visual cortex. *Vision Research.*\n",
+ "* van Beers RJ, Sittig AC, Gon JJ. (1999). Integration of proprioceptive and visual position-information: An experimentally supported model. *Journal of Neurophysiology.*"
+ "Author: J. R. Johansson (robert@riken.jp), http://dml.riken.jp/~rob/\n",
+ "\n",
+ "The latest version of this [IPython notebook](http://ipython.org/ipython-doc/dev/interactive/htmlnotebook.html) lecture is available at [http://github.com/jrjohansson/qutip-lectures](http://github.com/jrjohansson/qutip-lectures).\n",
+ "\n",
+ "The other notebooks in this lecture series are indexed at [http://jrjohansson.github.com](http://jrjohansson.github.com)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# setup the matplotlib graphics library and configure it to show \n",
+ "# figures inline in the notebook\n",
+ "%matplotlib inline\n",
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# make qutip available in the rest of the notebook\n",
+ "from qutip import *\n",
+ "\n",
+ "from IPython.display import Image"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Introduction and model\n",
+ "\n",
+ "Consider a single atom coupled to a single cavity mode, as illustrated in the figure below. If there atom excitation rate $\\Gamma$ exceeds the relaxation rate, a population inversion can occur in the atom, and if coupled to the cavity the atom can then act as a photon pump on the cavity."
+ "### Create a list of collapse operators that describe the dissipation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# collapse operators\n",
+ "c_ops = []\n",
+ "\n",
+ "rate = kappa * (1 + n_th_a)\n",
+ "if rate > 0.0:\n",
+ " c_ops.append(sqrt(rate) * a)\n",
+ "\n",
+ "rate = kappa * n_th_a\n",
+ "if rate > 0.0:\n",
+ " c_ops.append(sqrt(rate) * a.dag())\n",
+ "\n",
+ "rate = gamma\n",
+ "if rate > 0.0:\n",
+ " c_ops.append(sqrt(rate) * sm)\n",
+ "\n",
+ "rate = Gamma\n",
+ "if rate > 0.0:\n",
+ " c_ops.append(sqrt(rate) * sm.dag())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Evolve the system\n",
+ "\n",
+ "Here we evolve the system with the Lindblad master equation solver, and we request that the expectation values of the operators $a^\\dagger a$ and $\\sigma_+\\sigma_-$ are returned by the solver by passing the list `[a.dag()*a, sm.dag()*sm]` as the fifth argument to the solver."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "opt = Odeoptions(nsteps=2000) # allow extra time-steps \n",
+ "In order to study population genetics, we first need a model of a population. And even before that, we need to define what we mean by *population*. Populations can be defined on many levels and with many diffferent criteria. For our purposes, we will simply say that a population is a set of individuals sharing a common environment. And because this is population *genetics* we can think of individuals as entities comprising of specific genes or chromosomes. \n",
+ "\n",
+ "So where do we get a population from? As you may have discussed in previous workshops, there are very large datasets containing sequencing information from different populations. So we could download one of these datasets and perform some analysis on it. But I find this can be dry and tedious. So why download data when we can simply create our own?\n",
+ "\n",
+ "In this workshop we're going to be creating and studying our own \"artificial\" populations to illustrate some important population genetics concepts and methodologies. Not only will this help you learn population genetics, but you will get a lot more programming practice than if we were to simply parse data files and go from there. \n",
+ "\n",
+ "More specifically, we're going to build our own RNA world.\n",
+ "\n",
+ "As you may know, RNA is widely thought to be the first self replicating life-form to arise x billion years ago. One of the strongest arguments for this theory is that RNA is able to carry information in its nucleotides like DNA, and like protein, it is able to adopt higher order structures to catalyze reactions, such as self replication. So it is likely, and there is growing evidence that this is the case, that the first form of replicating life was RNA. And because of this dual property of RNA as an information vessel as well as a structural/functional element we can use RNA molecules to build very nice population models. \n",
+ "\n",
+ "So in this notebook, I'll be walking you through building genetic populations, simulating their evolution, and using statistics and other mathematical tools for understanding key properties of populations.\n",
+ "\n",
+ "### Building an RNA population\n",
+ "\n",
+ "As we saw earlier, RNA has the nice property of posessing a strong mapping between information carrying (sequence) and function (structure). This is analogous to what is known in evolutionary terms as a genotype and a phenotype. With these properties, we have everything we need to model a population, and simulate its evolution.\n",
+ "\n",
+ "#### RNA sequence-structure\n",
+ "\n",
+ "We can think of the genotype as a sequence $s$ consisting of letters/nucleotides from the alphabet $\\{U,A,C,G\\}$. The corresponding phenotype $\\omega$ is the secondary structure of $s$ which can be thought of as a pairing between nucleotides in the primary sequence that give rise to a 2D architecture. Because it has been shown that the function of many biomolecules, including RNA, is driven by structure this gives us a good proxy for phenotype. \n",
+ "\n",
+ "Below is an example of what an RNA secondary structure, or pairing, looks like."
+ "As you can see, unparied positions are forming loop-like structures, and paired positions are forming stem-like structures. It is this spatial arrangement of nucleotides that drives RNA's function. Therefore, another sequence that adopts a similar shape, is likely to behave in a similar manner. Another thing to notice is that, although in reality this is often not the case, in general we only allow pairs between $\\{C,G\\}$ and $\\{A, U\\}$ nucleotides, most modern approaches allow for non-canonical pairings and you will find some examples of this in the above structure.\n",
+ "\n",
+ "*How do we go from a sequence to a structure?*\n",
+ "\n",
+ "So a secondary structure is just a list of pairings between positions. How do we get the optimal pairing?\n",
+ "\n",
+ "The algorithm we're going to be using in our simulations is known as the Nussinov Algorithm. The Nussinov algorithm is one of the first and simplest attempts at predicting RNA structure. Because bonds tend to stabilize RNA, the algorithm tries to maximize the number of pairs in the structure and return that as its solution. Current approaches achieve more accurate solutions by using energy models based one experimental values to then obtain a structure that minimizes free energy. But since we're not really concerned with the accuracy of our predictions, Nussinov is a good entry point. Furthermore, the main algorithmic concepts are the same between Nussinov and state of the art RNA structure prediction algorithms. I implemented the algorithm in a separate file called `fold.py` that we can import and use its functions. I'm not going to go into detail here on how the algorithm works because it is beyond the scope of this workshop but there is a bonus exercise at the end if you're curious.\n",
+ "\n",
+ "You can predict a secondary structure by calling `nussinov()` with a sequence string and it will return a tuple in the form `(structure, pairs)`."
+ "You will see a funny dot-bracket string in the output. This is a representation of the structure of an RNA. Quite simply, a matching parir of parentheses (open and close) correspond to the nucleotides at those positions being paried. Whereas, a dot means that that position is unpaired in the structure. Feel free to play around with the input sequence to get a better understanding of the notation.\n",
+ "\n",
+ "So that's enough about RNA structure prediction. Let's move on to building our populations.\n",
+ "\n",
+ "### Fitness of a sequence: Target Structure\n",
+ "\n",
+ "Now that we have a good way of getting a phenotype (secondary structure), we need a way to evaluate the fitness of that phenotype. If we think in real life terms, fitness is the ability of a genotype to replicate into the next generation. If you have a gene carrying a mutation that causes some kind of disease, your fitness is decreased and you have a lower chance of contributing offspring to the next generation. On a molecular level the same concept applies. A molecule needs to accomplish a certain function, i.e. bind to some other molecule or send some kind of signal. And as we've seen before, the most important factor that determines how well it can carry out this function is its structure. So we can imagine that a certain structure, we can call this a 'target' structure, is required in order to accomplish a certain function. So a sequence that folds correctly to a target structure is seen as having a greater fitness than one that does not. Since we've encoded structures as simple dot-bracket strings, we can easily compare structures and thus evaluate the fitness between a given structure and the target, or 'correct' structure. \n",
+ "\n",
+ "There are many ways to compare structures $w_{1}$ and $w_{2}$, but we're going to use one of the simplest ways, which is base-pair distance. This is just the number of pairs in $w_{1}$ that are not in $w_{2}$. Again, this is beyond the scope of this workshop so I'll just give you the code for it and if you would like to know more you can ask me."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "('((....))', [(0, 7), (1, 6)])\n",
+ "('.(....).', [(1, 6)])\n",
+ "1\n"
+ ]
+ }
+ ],
+ "source": [
+ "### 3\n",
+ "\n",
+ "#ss_to_bp() and bp_distance() by Vladimir Reinharz.\n",
+ "def ss_to_bp(ss):\n",
+ " bps = set()\n",
+ " l = []\n",
+ " for i, x in enumerate(ss):\n",
+ " if x == '(':\n",
+ " l.append(i)\n",
+ " elif x == ')':\n",
+ " bps.add((l.pop(), i))\n",
+ " return bps\n",
+ "\n",
+ "def bp_distance(w1, w2):\n",
+ " \"\"\"\n",
+ " return base pair distance between structures w1 and w1. \n",
+ " w1 and w1 are lists of tuples representing pairing indices.\n",
+ "#give the list of pairs to bp_distance and see what the distance is.\n",
+ "print(bp_distance(w1[-1], w2[-1]))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Defining a cell: a little bit of Object Oriented Programming (OOP)\n",
+ "\n",
+ "Since we're going to be playing aroudn with sequences and structures and fitness values a lot, it's best to package it all nicely into an object. As you'll have seen with Vlad, objects are just a nice way of grouping data into an easily accessible form. \n",
+ "\n",
+ "We're trying to simulate evolution on a very simple kind of organism, or cell. It contains two copies of a RNA gene, each with a corresponding structure. "
+ "Now we've defined a 'Cell'. Since a population is a collection of individuals our populations will naturally consist of **lists** of 'Cell' objects, each with their own sequences. Here we initialize all the Cells with random sequences and add them to the 'population' list."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "### 5\n",
+ "import random\n",
+ "\n",
+ "def populate(target, pop_size=100):\n",
+ " \n",
+ " population = []\n",
+ "\n",
+ " for i in range(pop_size):\n",
+ " #get a random sequence to start with\n",
+ " sequence = \"\".join([random.choice(\"AUCG\") for _ in range(len(target))])\n",
+ " #use nussinov to get the secondary structure for the sequence\n",
+ " structure = nussinov(sequence)\n",
+ " #add a new Cell object to the population list\n",
+ "Now that we can store populatoins of cells, we need a way to evaluate the fitness of a given Cell. Recall that a Cell is simply an object that contains two RNA sequences (think of it as two copies of a gene on each chromosome). \n",
+ "\n",
+ "So we simply need to loop through each Cell in a population and compute base pair distance to the target structure. However, simply using base-pair distance is not a very good measure of fitness. There are two reasons for this: \n",
+ "\n",
+ "1. We want fitness to represent a *probability* that a cell will reproduce, and base pair distance is an integer.\n",
+ "2. We want this probability to be a *relative* measure. That is, we want to be the fitness to be proportional to how good a cell is with respect to all others in the population. This touches on an important principle in evolution where we only need to be 'better' than the competition and not good in some absolute measure. For example, if you and I are being chased by a bear. In order to survive, I only need to be faster than you, and not necessarily some absolute level of fitness.\n",
+ "\n",
+ "In order to get a probability (number between 0 and 1) we use the following equation to define the fitness of a structure $\\omega$ on a target structure $T$:\n",
+ "plt.plot([math.exp(beta * (bp_dist / float(target_length))) for bp_dist in range(target_length)])\n",
+ "plt.xlabel(\"Base pair distance to target structure\")\n",
+ "plt.ylabel(\"P(w, T)\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As you can see, it's a very simple function that evaluates to 1 (highest fitness) if the base pair distance is 0, and decreases as the structures get further and further away from the target. I didn't include the $N$ in the plotting as it will be a bit more annoying to compute, but it is simply a scaling factor so the shape and main idea won't be different.\n",
+ "\n",
+ "Now we can use this function to get a fitness value for each Cell in our population."
+ "Evolution would go nowhere without random mutations. While mutations are technically just random errors in the copying of genetic material, they are essential in the process of evolution. This is because they introduce novel diversity to populatons, which with a low frequency can be beneficial. And when a beneficial mutation arises (i.e. a mutation that increases fitness, or replication probability) it quickly takes over the population and the populatioin as a whole has a higher fitness.\n",
+ "\n",
+ "Implementing mutations in our model will be quite straightforward. Since mutations happen at the genotype/sequence level, we simply have to iterate through our strings of nucleotides (sequences) and randomly introduce changes."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "AAAAGGAGUGUGUAUGU\n",
+ "('AcAAGgAuUGUuaAaGa', True)\n"
+ ]
+ }
+ ],
+ "source": [
+ "def mutate(sequence, mutation_rate=0.001):\n",
+ " \"\"\"Takes a sequence and mutates bases with probability mutation_rate\"\"\"\n",
+ " \n",
+ " #start an empty string to store the mutated sequence\n",
+ " new_sequence = \"\"\n",
+ " #boolean storing whether or not the sequence got mutated\n",
+ " mutated = False\n",
+ " #go through every bp in the sequence\n",
+ " for bp in sequence:\n",
+ " #generate a random number between 0 and 1\n",
+ " r = random.random()\n",
+ " #if r is below mutation rate, introduce a mutation\n",
+ " if r < mutation_rate:\n",
+ " #add a randomly sampled nucleotide to the new sequence\n",
+ "The final process in this evolution model is selection. Once you have populations with a diverse range of fitnesses, we need to select the fittest individuals and let them replicate and contribute offspring to the next generation. In real populations this is just the process of reproduction. If you're fit enough you will be likely to reproduce more than another individual who is not as well suited to the environment.\n",
+ "\n",
+ "In order to represent this process in our model, we will use the fitness values that we assigned to each Cell earlier and use that to select replicating Cells. This is equivalent to sampling from a population with the sampling being weighted by the fitness of each Cell. Thankfully, `numpy.random.choice` comes to the rescue here. Once we have sampled enough Cells to build our next generation, we introduce mutations and compute the fitness values of the new generation."
+ "## Gathering information on our populations\n",
+ "\n",
+ "Here we simply store some statistics (in a dictionary) on the population at each generation such as the average base pair distance and the average fitness of the populations. No coding to do here, it's not a very interesting function but feel free to give it a look."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def record_stats(pop, population_stats):\n",
+ " \"\"\"\n",
+ " Takes a population list and a dictionary and updates it with stats on the population.\n",
+ " \"\"\"\n",
+ " generation_bp_distance_1 = [rna.bp_distance_1 for rna in pop]\n",
+ " generation_bp_distance_2 = [rna.bp_distance_2 for rna in pop]\n",
+ "Let's see if it actually worked by plotting the average base pair distance as a function of generations for both genes in each cell. We should expect a gradual decrease as the populations get closer to the target structure."
+ " Plot base pair distance for each chromosome over generations.\n",
+ " \"\"\"\n",
+ " for m in ['mean_bp_distance_1', 'mean_bp_distance_2']:\n",
+ " plt.plot(pops_stats[m], label=m)\n",
+ " plt.legend()\n",
+ " plt.xlabel(\"Generations\")\n",
+ " plt.ylabel(\"Mean Base Pair Distance\")\n",
+ " \n",
+ "evo_plot(pops_stats)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "You should see a nice drop in base pair distance! Another way of visualizing this is by plotting a histogram of the base pair distance of all Cells in the initial population versus the final population."
+ " Plots histograms of base pair distance in initial and final populations.\n",
+ " \"\"\"\n",
+ " #plot bp_distance_1 for rnas in first population\n",
+ " g = sns.distplot([rna.bp_distance_1 for rna in pops[0]], label='initial population')\n",
+ " #plot bp_distance_1 for rnas in first population\n",
+ " g = sns.distplot([rna.bp_distance_1 for rna in pops[-1]], label='final population')\n",
+ " g.set(xlabel='Mean Base Pair Distance')\n",
+ " g.legend()\n",
+ "bp_distance_distributions(pops)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Studying our evolved sequences with some Population Genetics tools\n",
+ "\n",
+ "Now that we've generated some sequences, we can analyze them!\n",
+ "\n",
+ "So after several rounds of selection, what do we get? We have a bunch of different sequences. We would like a way to characterize this diversity. One important tool for doing this is by making what is known as phylogenetic trees. \n",
+ "\n",
+ "Phylogenetic trees tell us about which groups of similar sequences are present and how they are likely related in evolutionary time. \n",
+ "\n",
+ "There are several ways of building phylogenetic trees using BioPython. Here we will go over one type and I'll leave another one as an exercise.\n",
+ "\n",
+ "### UPGMA (Unweighted Pair Group Method with Arithmetic Means)\n",
+ "\n",
+ "This is basically a clustering method based on the distance (or number of differences) between every pair of sequences. It assumes that sequences that are more similar are more likely to be related than the other way around. \n",
+ "\n",
+ "For $N$ sequences, the algorithm builds an $NxN$ matrix that stores the distance between each sequence to every other sequence. The algorithm goes through this matrix and finds the pair of sequences that is most similar and merges it into a 'cluster' or in tree terms, connects them to a common node. This process is repeated until all the sequences have been assigned to a group. Refer to the wikipedia article on [UPGMA](https://en.wikipedia.org/wiki/UPGMA) for a more detailed explanation. \n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {
+ "slideshow": {
+ "slide_type": "-"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from Bio import SeqIO\n",
+ "from Bio.Seq import Seq\n",
+ "from Bio.SeqRecord import SeqRecord\n",
+ "from Bio import AlignIO"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sequences = []\n",
+ "#let's take the first 10 sequences of our population to keep things simple\n",
+ "for seq in pops[-1][:10]:\n",
+ " #store each sequence in the sequences list as a SeqRecord object\n",
+ "The UPGMA algorithm requires a `MultipleSeqAlignment` object to build the distance matrix. So now that we have the `seq.fasta` file, we can give it to an online multiple sequence alignment tool. We can do this through BioPython but it requires some installation and setup so we will skip that for now. Go to the [MUSCLE Web Server](http://www.ebi.ac.uk/Tools/msa/muscle/) and give it the `seq.fasta` file. It will take a few seconds and it will give you an alignment and click *Download Alignment File*, copy paste the whole thing to a new file called `aln.clustal`. This is the alignment we will use to build our tree."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#open the alignmnent file\n",
+ "with open(\"aln.clustal\", \"r\") as aln:\n",
+ " #use AlignIO to read the alignment file in 'clustal' format\n",
+ "The populations we generated evolved asexually. This means that individuals do not mate or exchange genetic information. So to make our simulation a bit more interesting let's let the Cells mate. This is going to require a few small changes in the `selection()` function. Previously, when we selected sequences to go into the next generation we just let them provide one offspring which was a copy of itself and introduced mutations. Now instead of choosing one Cell at a time, we will randomly choose two 'parents' that will mate. When they mate, each parent will contribute one of its chromosomes to the child. We'll repeat this process until we have filled the next generation."
+ "When we are presented with data from a population we don't know much about. It is often useful to try to learn whether there are any evolutionary or behavioural influences that are shaping population dynamics. This could be in the form of selective pressure, mating preference, genetic drift, mutations, gene flow, etc. So in order to detect if something like this is happening we need to develop a test. This is where Hardy Weinberg comes in. \n",
+ "\n",
+ "The Hardy Weinberg equilibrium states that \"allele and genotype frequencies remain constant in the absence of other evolutionary influences. (such as the ones we mentioned above)\" - Wikipedia.\n",
+ "\n",
+ "So if we can measure allele/genotype frequencies (which we can do because we have sequences), we can see whether the HW principle holds true. If it does not, then we can do more digging to see what could be happening to shift populations away from equilibrium.\n",
+ "\n",
+ "In order to do this we need to define an 'allele'. An allele (for our purproses) will be a locus (position in a sequence) that can take one of two states, a *reference* state or an *alternate* state. For example, we can look at locus number **5** (position 5 in our RNA sequences) and call reference **C**, and alternate **G**. If we are in HW we can predict the frequency of each allele in our population.\n",
+ "\n",
+ "To simplify our notation we will call the alternate allele *A* and the reference allele *a*. We can write the probability of each allele as $p_{A} + p_{a} = 1$. Since we are dealing with diploid populations, each individual will have two copies of each locus so it can be $p_{AA}, p{Aa}, p{aA}, p{aa}$. By simple probability laws we can get an expression for the probability of each genotype based on the probabilities of the single loci $p_{a}$ and $p_{A}$.\n",
+ "\n",
+ "$$p_{aa}\\simeq p_{a}^2$$\n",
+ "\n",
+ "$$p_{AA}\\simeq p_{A}^2$$\n",
+ "\n",
+ "$$p_{Aa,~aA} \\simeq 2 p_{a} p_{A}.$$\n",
+ "\n",
+ "Since it is hard to know what the true probability of observing either $p_{a}$ and $p_{A}$ we can estimate this probability from our data as follows:\n",
+ "Now we have expected values, and observed values. We need a test to determine whether we have a significant departure from the hypothesis of Hardy Weinberg equilibrium. The statistical test that is commonly used is known as the $\\chi^{2}$ test. If you take a look at the equation you'll see that the statistic simply takes the squared difference between our observed value and the expected value (divided by expected) and sums this for each possible genotype. The reason we take the squared difference is because we want to deal only with positive values, hence the name $\\chi^{2}$.\n",
+ "The first thing we need to do is get alleles from our sequence data. This boils down to going through each sequence at the position of interest and counting the number of $AA$, $Aa$, $aa$ we get.\n",
+ "\n",
+ "\n",
+ "\\** the sections on Hardy Weinberg and F-statistics are adapted from Simon Gravel's HGEN 661 Notes"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def allele_finder(pop, locus, ref, alt):\n",
+ " genotypes = []\n",
+ " for p in pop:\n",
+ " #get the nucleotide at the locus from the first chromosome \n",
+ " locus_1 = p.sequence_1[locus].upper()\n",
+ " #same for the second\n",
+ " locus_2 = p.sequence_2[locus].upper()\n",
+ " \n",
+ " #check that it is either ref or alt, we don't care about other alleles for now.\n",
+ " if locus_1 in (ref, alt) and locus_2 in (ref, alt):\n",
+ " #if the alelle is ref, store a value of 1 in allele_1, and 0 otherwise\n",
+ " allele_1 = int(locus_1 == ref)\n",
+ " #same for the second allele\n",
+ " allele_2 = int(locus_2 == ref)\n",
+ " \n",
+ " #add allele to our list of alleles as a tuple. \n",
+ "Now that we have alleles represented in the right form, we can see if our population is at Hardy Weinberg equilibrium using the $\\chi_{2}$ test and the equations above."
+ "Can we say that our population is at equilibrium? Can you find parameters for `evolution_with_mating()` that will give us populations outside of the HW equilibrium?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## A brief interlude on the p-value\n",
+ "\n",
+ "Let's take a minute to understand what the p-value means. The p-value is a probability. Specifically, it is the probability of observing a value equal to or more extreme than that our statistic given the test distribution. So in our case, it is the probability of observing a $X^2$ greater than or equal to the one the test gives us under a $\\chi^2$ distribution. When this value is very small, it suggests that it is unlikely that we are sampling from our assumed 'null' distribution and that some other alternate distribution is the true distribution. So a low p-value here would be evidence against the neutral Hardy Weinberg model and would suggest that our population is experiencing some influences such as mating preference, selection, mutation etc.\n",
+ "\n",
+ "A lot of research bases its conclusions solely on p-value and it is important to be very wary of this bad practice. It has become a bad convention that people say a p-value lower than some arbitrary threshold means one's findings are significant. However, very often the p-value does not give us the whole story and we need to know about things like sample size, size of impact, reproducibility, power of the test, etc. (check this out [American Statistical Association statement on p-values](http://www.nature.com/news/statisticians-issue-warning-over-misuse-of-p-values-1.19503), [p-hacking](http://fivethirtyeight.com/features/science-isnt-broken/#part1), and [this](http://allendowney.blogspot.ca/2016/06/there-is-still-only-one-test.html))\n",
+ "\n",
+ "Let's just visualize this very quickly using the $\\chi^{2}_{1}$ distribution. You will see that the p-value corresponds to the shaded red area under the curve. That area is the probability of observing a value as extreme or more than the one we found. When that is a very small area, we can be more confident that our assumption of HW is false."
+ "The last topic we'll cover is F-statistics. \n",
+ "\n",
+ "Once we find that our population strays from the HW condition we can begin to ask why that is the case. Often this deviation from the expected allele frequencies under HW is due to mating preference. Hardy Weinberg assumes that all individuals in a population have an equal probability of mating with any other individual (random mating). However, when certain individuals prefer to mate with specific others (in real populations this can be due to culture, race, geographic barriers, etc.), you get what is known as population structure. Population structure means that we begin to see *sub-populations* within our total population where individuals prefer to mate within their sub-population. This biased mating will result in a higher number of homozygotes than we would expect under Hardy-Weinberg equilibrium. Simply because mating preferences will tend to drive populations toward similar genotypes. So if this is the case, and no other factors are biasing allele dynamics, within sub-populations we should have Hardy-Weinberg like conditions. \n",
+ "\n",
+ "For example, if Raptors fans prefer to mate with other Raptors fans, then when we consider only Raptors fans, we should observe random mating. Simply because if the mating preference criterion is 'being a Raptor's fan' then any Raptor's fan will be equally likely to mate with any other Raptor's fan so we have Hardy Weinberg again.\n",
+ "\n",
+ "Let's express this in quantities we can measure.\n",
+ "\n",
+ "From before we calculated the observed and expected number of heterozygotes in a population. Let's call these $\\hat H$ and $H_{IT}$ respectively. $\\hat H$ is just the count of heterozygotes, and $H_{IT}$ is the same as the expected number of heterozygotes we calculated earlier.\n",
+ "\n",
+ "We define a quantity $e_{IT}$ as a measure of the 'excess heterozygosity' in the population when we consider all individuals $I$ in the total population $T$. $e_{IT} > 1$ when we have more heterozygotes than we expect under HW. And $0 < e_{IT} < 1$ if we have less heterozygotes than we would expect under HW.\n",
+ "\n",
+ "\n",
+ "$$e_{IT}=\\frac{\\mbox{observed proportion of hets}}{\\mbox{expected proportion of hets}}=\\frac{ H_{obs}}{H_{IT}}$$\n",
+ "\n",
+ "We use $e_{IT}$ to define the statistic $F_{IT}$\n",
+ "\n",
+ "$$F_{IT}=1-e_{IT}$$\n",
+ "\n",
+ "So $F_{IT} > 0$ when we have a lack of heterozygotes and $F_{IT} < 0$ when we have an excess of heterozygotes. $F_{IT} = 0$ under random mating.\n",
+ "\n",
+ "When we have a subpropulation $S$ we can calculate the equivalent quantity but instead of considering heterozygosity in the whole population we only take a sub-population into account.\n",
+ "\n",
+ "$$e_{IS} = \\frac{H_{obs}}{H_{IS}}$$\n",
+ "\n",
+ "And lastly, we have $F_{ST}$. This one is not as intuitive to derive so I'm not including the derivation here. But basically it measure the excess heterozygosity in the total population due to the presence of two subpopulations with allele frequencies $p_{1}$ and $p_{2}$.\n",
+ "Let's get some data for our F-tests. First we need to evolve two populations indepenently of each other, to simulate isolated mating. Then to simulate the total population we combine the two sub-populations. We then use our `allele_finder()` function to get all the alleles, and the `hardy_weinberg_chi_2_test()` function to get our expected and observed counts. Finally we plug those into the `f_statistics()` function."
+ "Try playing with different evolution parameters and see the effect on the different F-statistics. This workshop is a work in progress so there may be some biases in our simulation scheme that can make for come confusing F-statistics. If you come up with anything interesting I would love to know about it."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Exercises / Extra Material\n",
+ "\n",
+ "### Programming Exercises\n",
+ "\n",
+ "i. *Heatmap of mutation rates vs. population sizes.* (short) \n",
+ "\n",
+ "Make a heatmap that plots the base pair distance of the average base pair distance of the population at generation `-1` for mutation rates $\\mu = \\{0, 0.001, 0.01, 0.1, 0.5\\}$ and population sizes $N=\\{10, 100, 1000, 10000\\}$. The resulting heatmap will be `5x4` dimensions. You may choose how many generations to evolve your populations, just plot the last one in the heatmap."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "([<matplotlib.axis.YTick at 0x7fa02ba18ac8>,\n",
+ " <matplotlib.axis.YTick at 0x7fa02ba183c8>,\n",
+ " <matplotlib.axis.YTick at 0x7fa02ba59a58>,\n",
+ " <matplotlib.axis.YTick at 0x7fa02acbe978>,\n",
+ " <matplotlib.axis.YTick at 0x7fa02acc5048>],\n",
+ "#xticks/yticks takes a list of numbers that specify the position of the ticks and a list with the tick labels\n",
+ "plt.xticks([i + .5 for i in range(len(population_sizes))], population_sizes)\n",
+ "plt.ylabel(\"Mutation Rate\")\n",
+ "plt.yticks([i + .5 for i in range(len(mutation_rates))], mutation_rates)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "ii. *Introduce mating preferences within a population.* (medium length) \n",
+ "\n",
+ "Modify the `selection_with_mating()` function to allow for mating preferences within a population. In our example above we were just running two independent simulations to study barriers to gene flow. But now you will implement mating preferences within a single simulation. Your function will assign each Cell a new attribute called `self.preference` which will take a string value denoting the mating type the current cell prefers to mate with. For example we can have a population with three mating types: $\\{A, B, C\\}$. Your function will randomly assign preferences to each cell in the initial population. We will define a preference between types $A$ and $B$ as the probability that two cells of those given types will mate if selected. \n",
+ "\n",
+ "$$\n",
+ "preferences(A,B,C) = \n",
+ "\\begin{bmatrix}\n",
+ " 0.7 & 0.1 & 0.2 \\\\\n",
+ " 0.1 & 0.9 & 0 \\\\\n",
+ " 0.2 & 0 & 0.8 \\\\\n",
+ "\\end{bmatrix}\n",
+ "$$\n",
+ "\n",
+ "Once you selected two potential parents for mating (as we did earlier) you will use the matrix to evaluate whether or not the two parents will mate and contribute an offspring to the next generation. "
+ "### Population Genetics / Bioinformatics Exercises\n",
+ "\n",
+ "*Exercise 1. Make a tree using maximum parsimony.*\n",
+ "\n",
+ "We saw how to make trees using a distance score. Another popular method is known as the maximum parsimony approach. I won't go into too much detail on this since we are short on time, but I will give a quick intro and we'll look at how ot make a tree using maximum parsimony.\n",
+ "\n",
+ "This approach is based on the principle of parsimony, which states that the simplest explanation for our data is the most likely to be true. So given an alignment, we assume that the best tree is the one that minimizes the number of changes, or mutations. This is often a reasonable assumption to make since mutation rates in real populations are generally low, and things like back-mutations (e.g. A --> C --> A) are unlikely. Computing the tree that that maximizes parsimony directly is a difficult task, but evaluating the parsimony score of a tree given the tree is easy. So this approach basically generates many random trees for the data and scores them based on parsimony keeping the most parsimonious tree. Take a look at [the biopython manual to work through this example](http://biopython.org/wiki/Phylo), and [this one](http://biopython.org/DIST/docs/api/Bio.Phylo.TreeConstruction.ParsimonyTreeConstructor-class.html).\n",
+ "\n",
+ "Since we already have an alignment (`aln.clustal`) we will just re-use it and make a maximum parsimony tree instead. "
+ "We just saw two methods of growing phylogenetic trees given an alignment. However, as we saw with the maximum parsimony approach, there can be many different trees for a single data set. How do we know our tree is a good representation of the data? By 'good' here we will instead use the word 'robust'. Is the tree we use too sensitive to the particularities of the data we gave it? If we make a small change in the sequence will we get a very different tree? Normally these problems would be addressed by re-sampling and seeing if we obtain similar results. But we can't really re-sample evolution. It happened once and we can't make it happen again. So we use something called *bootstrapping* which is a technique often used in statistics where instead of generating new data, you re-sample from your present data.\n",
+ "\n",
+ "So we have a multiple sequence alignment with $M$ sequences (rows) each with sequences of length $N$ nucleotides (columns). For each row, we can randomly sample $N$ nucleotides with replacement to make a new 'bootstrapped' sequence also of length $N$. Think of it as a kind of shuffling of the data. This gives us a whole new alignment that we can again use to make a new tree.\n",
+ "\n",
+ "This process is repeated many times to obtain many trees. The differences in topology (shape/structure) of the trees we obtained are assessed. If after this shuffling/perturbations we still get similar enough looking trees we can say that our final tree is robust to small changes in the data. ([some more reading on this](http://projecteuclid.org/download/pdf_1/euclid.ss/1063994979))\n",
+ "\n",
+ "Let's run a small example of this using the bootstrapping functions in `BioPython`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 39,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CLUSTAL X (1.81) multiple sequence alignment\n",
+ "\n",
+ "\n",
+ "8 GAAGAGACAC\n",
+ "7 GGGGGAGCGC\n",
+ "0 GCCACAGCGC\n",
+ "3 ACCACAGUGU\n",
+ "6 CCCACAGUGU\n",
+ "9 CCCACUAGAG\n",
+ "1 CCCUCUCGCG\n",
+ "2 CCCUCUCGCG\n",
+ "4 CUUUUUCGCG\n",
+ "5 CUUUUUCGCG\n",
+ " \n",
+ "\n",
+ "\n",
+ "\n",
+ "CLUSTAL X (1.81) multiple sequence alignment\n",
+ "\n",
+ "\n",
+ "8 GGAUUUCUUA\n",
+ "7 GAGCCCCCCG\n",
+ "0 AAGGGGCGGG\n",
+ "3 AAGGGGUGGG\n",
+ "6 AAGGGGUGGG\n",
+ "9 AUAUUUGUUA\n",
+ "1 UUCUUUGUUC\n",
+ "2 UUCUUUGUUC\n",
+ "4 UUCUUUGUUC\n",
+ "5 UUCUUUGUUC\n",
+ " \n",
+ "\n",
+ "\n",
+ "\n",
+ "CLUSTAL X (1.81) multiple sequence alignment\n",
+ "\n",
+ "\n",
+ "8 UCUUGAGAUC\n",
+ "7 CCCGGGGGGC\n",
+ "0 GCGAGGGCAC\n",
+ "3 GUGCAGACCU\n",
+ "6 GUGCCGCCCU\n",
+ "9 UGUACACCAG\n",
+ "1 UGUGCCCCGG\n",
+ "2 UGUGCCCCGG\n",
+ "4 UGUGCCCUGG\n",
+ "5 UGUGCCCUGG\n",
+ " \n",
+ "\n",
+ "\n",
+ "\n",
+ "CLUSTAL X (1.81) multiple sequence alignment\n",
+ "\n",
+ "\n",
+ "8 GGUGCGCAUG\n",
+ "7 GACGCAUGGG\n",
+ "0 AAGGUAAGAA\n",
+ "3 AAGAUAUGCA\n",
+ "6 AAGCUAUGCA\n",
+ "9 AUUCUUAAAA\n",
+ "1 UUUCAUACGU\n",
+ "2 UUUCAUACGU\n",
+ "4 UUUCAUACGU\n",
+ "5 UUUCAUACGU\n",
+ " \n",
+ "\n",
+ "\n",
+ "\n",
+ "CLUSTAL X (1.81) multiple sequence alignment\n",
+ "\n",
+ "\n",
+ "8 GCUAGGAACC\n",
+ "7 GCGGGGGGCU\n",
+ "0 AUACGGGGCA\n",
+ "3 AUCCAAGGUU\n",
+ "6 AUCCCCGGUU\n",
+ "9 AUACCCAAGA\n",
+ "1 UAGCCCCCGA\n",
+ "2 UAGCCCCCGA\n",
+ "4 UAGUCCCCGA\n",
+ "5 UAGUCCCCGA\n",
+ " \n",
+ "\n",
+ "\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "from Bio.Phylo.Consensus import *\n",
+ "\n",
+ "#open our alignment file.\n",
+ "with open('aln.clustal', 'r') as align:\n",
+ " aln = AlignIO.read(align, 'clustal')\n",
+ "\n",
+ "#take 5 bootstrap samples from our alignment\n",
+ "bootstraps = bootstrap(aln,5)\n",
+ "\n",
+ "#let's print each new alignment in clustal format. you should see 5 different alignments.\n",
+ "Similarly to the $\\chi^{2}$ test we saw for testing deviations from HW equilibrium, we can use a T-test to compare differences in means between two independent samples. We can use this to revisit a the first programming question in the exercsies section. Does mutation rate and population size have an effect on the fitness of populations? We can translate this question to, is there a difference in the mean base pair distance between populations under different mutation and population size regimes?\n",
+ "\n",
+ "Scipy has a very useful function that implements the T-test called `scipy.stats.ttest_ind`. Run two independent simulations (with different mutation rates) and compute the difference in mean bp distance between the two at their final generation. Store the populations in two different variables. Give a list of `bp_distance_1` values for each memeber of the population to `ttest_ind()`. \n",
+ "\n",
+ "Make sure to read teh `ttest_ind()` documentation, particularly about the argumetn `equal_var`. What should we set it to?"
+ "1. *Nussinov Algorithm (Only try this if you are feeling brave and are done with the other exercises or are interested in getting a taste of Computer Science. It is beyond the scope of this workshop.)*\n",
+ "\n",
+ "There are several approaches for solving this problem, we will look at the simplest one here which is known as the Nussinov Algorithm. This algorithm is a popular example of a class of algorithms know as dynamic programming algorithms. The main idea behind these algorithms is that we can break down the problem into many subproblems which are easier to compute than the full problem. Once we have obtained the solution for the subproblems, we can retrieve the solution to the full problem by doing something called a backtrace (more on the backtrace later). \n",
+ "\n",
+ "Here, the problem is obtaining the optimal pairing on a string of nucleotides. In order to know how good our structure is, we assign a score to it. One possible scoring scheme could be adding 1 to the score per paired set of nucleotides, and 0 otherwise. So in other words, we want a pairing that will give us the highest possible score. We can write this quantity as $OPT(i, j)$ where $i$ and $j$ are the indices of the sequence between which we obtain the pairing score. Our algorithm is therefore going to compute a folding score for all substrings bound by $i$ and $j$ and store the value in what is known as a dynamic programming table. Our dynamic programming table will be a $N$ x $N$ array where $N$ is the length of our sequence. So now that we have a way of measuring how good a structure is, we need a way to evaluate scores given a subsequence. To do this, we set some rules on the structure of an RNA sequence:\n",
+ "\n",
+ "\n",
+ "If $i$ and $j$ form a pair:\n",
+ "1. The pair $i$ and $j$ must form a valid watson-crick pair.\n",
+ "2. $i < j-4$. This ensures that bonding is not happening between positions that are too close to each other, which would produce steric clashes.\n",
+ "3. If pair $(i,j)$ and $(k, l)$ are in the structure, then $i < k < j < l$. This ensures that there is no crossing over of pairs which would result in pseudoknots.\n",
+ "4. No base appears in more than one pair.\n",
+ "\n",
+ "Using these rules we can begin to build our algorithm. The first part of our algorithm needs to take as input indices $i$ and $j$ and return the value $OPT(i,j)$ which is the optimal score of a structure between $i$ and $j$. We start by thinking about values of $i$ and $j$ for which we can immediately know the solution, this is known as a 'base case'. This is a case where the solution is known and no further recursion is required. Once the algorithm reaches the base case, it can return a solution and propagate it upward to the first recursive call. So once we have reached $i$ and $j$ that are too close to form a structure (rule number 2), we know that the score is 0. \n",
+ "\n",
+ "Otherwise, we must weigh the possibility of forming a pair or not forming a pair. If $i$ and $j$ are unpaired, then $OPT(i,j)$ is just $OPT(i, j-1)$ since the score will not increase for unpaired indices. \n",
+ "\n",
+ "The other case is that $i$ is paired to some index $t$ on the interval $[i,j]$. We then add 1 to the score and consider the structure formed before and after the pairing between $i$ and $t$. We can write these two cases as $OPT(i, t-1)$ and $OPT(t+1, j)$. But how do we know which $t$ to pair $i$ with? Well we simply try all possible values of $t$ within the allowed range and choose the best one. \n",
+ "\n",
+ "All of this can be summed up as follows:\n",
+ "\n",
+ "$$ OPT(i,j) = max\\begin{cases}\n",
+ " OPT(i, j-1) \\quad \\text{If $i$ and $j$ are not paired with each other.}\\\\\n",
+ " max(1 + OPT(i, t-1) + OPT(t+1, j)) \\quad \\text{Where we try all values of $t$ < j - 4}\n",
+ " \\end{cases}$$\n",
+ "\n",
+ "\n",
+ "We can now use this recursion to fill our dynamic programming table. Once we have filled the table with scores, we can retrieve the optimal folding by a process called backtracking. We won't go into detail on how this works, but the main idea is that we can start by looking at the entry containing the score for the full sequence $OPT[0][N]$. We can then look at adjacent entries and deduce which case (pairing or not pairing) resulted in the current value. We can continue like this for the full table until we have retrieved the full structure."
+ "# Rigid-body transformations in a plane (2D)\n",
+ "\n",
+ "> Marcos Duarte \n",
+ "> Laboratory of Biomechanics and Motor Control ([http://demotu.org/](http://demotu.org/)) \n",
+ "> Federal University of ABC, Brazil"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The kinematics of a rigid body is completely described by its pose, i.e., its position and orientation in space (and the corresponding changes are translation and rotation). The translation and rotation of a rigid body are also known as rigid-body transformations (or simply, rigid transformations).\n",
+ "\n",
+ "Remember that in physics, a [rigid body](https://en.wikipedia.org/wiki/Rigid_body) is a model (an idealization) for a body in which deformation is neglected, i.e., the distance between every pair of points in the body is considered constant. Consequently, the position and orientation of a rigid body can be completely described by a corresponding coordinate system attached to it. For instance, two (or more) coordinate systems can be used to represent the same rigid body at two (or more) instants or two (or more) rigid bodies in space.\n",
+ "\n",
+ "Rigid-body transformations are used in motion analysis (e.g., of the human body) to describe the position and orientation of each segment (using a local (anatomical) coordinate system defined for each segment) in relation to a global coordinate system fixed at the laboratory. Furthermore, one can define an additional coordinate system called technical coordinate system also fixed at the rigid body but not based on anatomical landmarks. In this case, the position of the technical markers is first described in the laboratory coordinate system, and then the technical coordinate system is calculated to recreate the anatomical landmarks position in order to finally calculate the original anatomical coordinate system (and obtain its unknown position and orientation through time).\n",
+ "\n",
+ "In what follows, we will study rigid-body transformations by looking at the transformations between two coordinate systems. For simplicity, let's first analyze planar (two-dimensional) rigid-body transformations and later we will extend these concepts to three dimensions (where the study of rotations are more complicated)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Affine transformations\n",
+ "\n",
+ "Translation and rotation are two examples of [affine transformations](https://en.wikipedia.org/wiki/Affine_transformation). Affine transformations preserve straight lines, but not necessarily the distance between points. Other examples of affine transformations are scaling, shear, and reflection. The figure below illustrates different affine transformations in a plane. Note that a 3x3 matrix is shown on top of each transformation; these matrices are known as the transformation matrices and are the mathematical representation of the physical transformations. Next, we will study how to use this approach to describe the translation and rotation of a rigid-body. \n",
+ "<br>\n",
+ "<figure><img src='https://upload.wikimedia.org/wikipedia/commons/thumb/2/2c/2D_affine_transformation_matrix.svg/360px-2D_affine_transformation_matrix.svg.png' alt='Affine transformations'/> <figcaption><center><i>Figure. Examples of affine transformations in a plane applied to a square (with the letter <b>F</b> in it) and the corresponding transformation matrices (<a href=\"https://en.wikipedia.org/wiki/Affine_transformation\">image from Wikipedia</a>).</i></center></figcaption> </figure>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Translation\n",
+ "\n",
+ "In a two-dimensional space, two coordinates and one angle are sufficient to describe the pose of the rigid body, totalizing three degrees of freedom for a rigid body. Let's see first the transformation for translation, then for rotation, and combine them at last.\n",
+ "\n",
+ "A pure two-dimensional translation of a coordinate system in relation to other coordinate system and the representation of a point in these two coordinate systems are illustrated in the figure below (remember that this is equivalent to describing a translation between two rigid bodies). \n",
+ "<br>\n",
+ "<figure><img src='./../images/translation2D.png' alt='translation 2D'/> <figcaption><center><i>Figure. A point in two-dimensional space represented in two coordinate systems (Global and local), with one system translated.</i></center></figcaption> </figure>\n",
+ "\n",
+ "The position of point $\\mathbf{P}$ originally described in the local coordinate system but now described in the Global coordinate system in vector form is:\n",
+ "Because position and translation can be treated as vectors, the inverse operation, to describe the position at the local coordinate system in terms of the Global coordinate system, is simply:\n",
+ "From classical mechanics, this transformation is an example of [Galilean transformation](http://en.wikipedia.org/wiki/Galilean_transformation). \n",
+ "\n",
+ "For example, if the local coordinate system is translated by $\\mathbf{L_G}=[2, 3]$ in relation to the Global coordinate system, a point with coordinates $\\mathbf{P_l}=[4, 5]$ at the local coordinate system will have the position $\\mathbf{P_G}=[6, 8]$ at the Global coordinate system:"
+ "This operation also works if we have more than one data point (NumPy knows how to handle vectors with different dimensions):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 6, 8],\n",
+ " [ 8, 10],\n",
+ " [10, 12]])"
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "Pl = np.array([[4, 5], [6, 7], [8, 9]]) # 2D array with 3 rows and two columns\n",
+ "PG = LG + Pl\n",
+ "PG"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Rotation\n",
+ "\n",
+ "A pure two-dimensional rotation of a coordinate system in relation to other coordinate system and the representation of a point in these two coordinate systems are illustrated in the figure below (remember that this is equivalent to describing a rotation between two rigid bodies). The rotation is around an axis orthogonal to this page, not shown in the figure (for a three-dimensional coordinate system the rotation would be around the $\\mathbf{Z}$ axis). \n",
+ "<br>\n",
+ "<figure><img src='./../images/rotation2D.png' alt='rotation 2D'/> <figcaption><center><i>Figure. A point in the two-dimensional space represented in two coordinate systems (Global and local), with one system rotated in relation to the other around an axis orthogonal to both coordinate systems.</i></center></figcaption> </figure>\n",
+ "\n",
+ "Consider we want to express the position of point $\\mathbf{P}$ in the Global coordinate system in terms of the local coordinate system knowing only the coordinates at the local coordinate system and the angle of rotation between the two coordinate systems. \n",
+ "\n",
+ "There are different ways of deducing that, we will see three of these methods next. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 1. Using trigonometry\n",
+ "\n",
+ "From figure below, the coordinates of point $\\mathbf{P}$ in the Global coordinate system can be determined finding the sides of the triangles marked in red. \n",
+ "<br>\n",
+ "<figure><img src='./../images/rotation2Db.png' alt='rotation 2D'/> <figcaption><center><i>Figure. The coordinates of a point at the Global coordinate system in terms of the coordinates of this point at the local coordinate system.</i></center></figcaption> </figure>\n",
+ "Where $\\mathbf{R_{Gl}}$ is the rotation matrix that rotates the coordinates from the local to the Global coordinate system:\n",
+ "\n",
+ "$$ \\mathbf{R_{Gl}} = \\begin{bmatrix}\n",
+ "\\cos\\alpha & -\\sin\\alpha \\\\\n",
+ "\\sin\\alpha & \\cos\\alpha \n",
+ "\\end{bmatrix} $$\n",
+ "\n",
+ "So, given any position at the local coordinate system, with the rotation matrix above we are able to determine the position at the Global coordinate system. Let's check that before looking at other methods to obtain this matrix. \n",
+ "\n",
+ "For instance, consider a local coordinate system rotated by $45^o$ in relation to the Global coordinate system, a point in the local coordinate system with position $\\mathbf{P_l}=[1, 1]$ will have the following position at the Global coordinate system:"
+ "PG = np.dot(RGl, Pl) # the function dot() is used for matrix multiplication with arrays\n",
+ "np.around(PG, 4) # round the number due to floating point arithmetic errors"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We have done all the calculations using the array function in NumPy. A [NumPy array is different than a matrix](http://www.scipy.org/NumPy_for_Matlab_Users), if we want to use explicit matrices in NumPy, the calculation above will be:"
+ "PG = RGl*Pl # matrix multiplication in NumPy\n",
+ "np.around(PG, 4) # round the number due to floating point arithmetic errors"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Both array and matrix types work in NumPy, but you should choose only one type and not mix them; the array is preferred because it is [the standard vector/matrix/tensor type of NumPy](http://www.scipy.org/NumPy_for_Matlab_Users)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 2. Using direction cosines\n",
+ "\n",
+ "Another way to determine the rotation matrix is to use the concept of direction cosine. \n",
+ "\n",
+ "> Direction cosines are the cosines of the angles between any two vectors. \n",
+ "\n",
+ "For the present case with two coordinate systems, they are the cosines of the angles between each axis of one coordinate system and each axis of the other coordinate system. The figure below illustrates the directions angles between the two coordinate systems, expressing the local coordinate system in terms of the Global coordinate system. \n",
+ "<br>\n",
+ "<figure><img src='./../images/directioncosine2D.png' alt='direction angles 2D'/> <figcaption><center><i>Figure. Definition of direction angles at the two-dimensional space.</i></center></figcaption> </figure> \n",
+ "<br>\n",
+ "$$ \\mathbf{R_{Gl}} = \\begin{bmatrix}\n",
+ "\\cos\\mathbf{X}x & \\cos\\mathbf{X}y \\\\\n",
+ "\\cos\\mathbf{Y}x & \\cos\\mathbf{Y}y \n",
+ "\\end{bmatrix} = \n",
+ "\\begin{bmatrix}\n",
+ "\\cos(\\alpha) & \\cos(90^o+\\alpha) \\\\\n",
+ "\\cos(90^o-\\alpha) & \\cos(\\alpha)\n",
+ "\\end{bmatrix} = \n",
+ "\\begin{bmatrix}\n",
+ "\\cos\\alpha & -\\sin\\alpha \\\\\n",
+ "\\sin\\alpha & \\cos\\alpha \n",
+ "\\end{bmatrix} $$ \n",
+ "\n",
+ "The same rotation matrix as obtained before.\n",
+ "\n",
+ "Note that the order of the direction cosines is because in our convention, the first row is for the $\\mathbf{X}$ coordinate and the second row for the $\\mathbf{Y}$ coordinate (the outputs). For the inputs, we followed the same order, first column for the $\\mathbf{x}$ coordinate, second column for the $\\mathbf{y}$ coordinate."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 3. Using a basis\n",
+ "\n",
+ "Another way to deduce the rotation matrix is to view the axes of the rotated coordinate system as unit vectors, versors, of a <a href=\"http://en.wikipedia.org/wiki/Basis_(linear_algebra)\">basis</a> as illustrated in the figure below.\n",
+ "\n",
+ "> A basis is a set of linearly independent vectors that can represent every vector in a given vector space, i.e., a basis defines a coordinate system.\n",
+ "\n",
+ "<figure><img src='./../images/basis2D2.png' alt='basis 2D'/> <figcaption><center><i>Figure. Definition of the rotation matrix using a basis at the two-dimensional space.</i></center></figcaption> </figure>\n",
+ "\n",
+ "The coordinates of these two versors at the local coordinate system in terms of the Global coordinate system are:\n",
+ "Note that as unit vectors, each of the versors above should have norm (length) equals to one, which indeed is the case.\n",
+ "\n",
+ "If we express each versor above as different columns of a matrix, we obtain the rotation matrix again: \n",
+ "\n",
+ "$$ \\mathbf{R_{Gl}} = \\begin{bmatrix}\n",
+ "\\cos\\alpha & -\\sin\\alpha \\\\\\\n",
+ "\\sin\\alpha & \\cos\\alpha \n",
+ "\\end{bmatrix} $$\n",
+ "\n",
+ "This means that the rotation matrix can be viewed as the basis of the rotated coordinate system defined by its versors. \n",
+ "\n",
+ "This third way to derive the rotation matrix is in fact the method most commonly used in motion analysis because the coordinates of markers (in the Global/laboratory coordinate system) are what we measure with cameras. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 4. Using the inner (dot or scalar) product between versors\n",
+ "\n",
+ "Yet another way to deduce the rotation matrix is to define it as the dot product between the versors of the bases related to the two coordinate systems:\n",
+ "And the rotation matrix will be equal to the matrix deduced based on the direction cosines."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Local-to-Global and Global-to-local coordinate systems' rotations"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If we want the inverse operation, to express the position of point $\\mathbf{P}$ in the local coordinate system in terms of the Global coordinate system, the figure below illustrates that using trigonometry. \n",
+ "<br>\n",
+ "<figure><img src='./../images/rotation2Dc.png' alt='rotation 2D'/> <figcaption><center><i>Figure. The coordinates of a point at the local coordinate system in terms of the coordinates at the Global coordinate system.</i></center></figcaption> </figure>\n",
+ "Where $\\mathbf{R_{lG}}$ is the rotation matrix that rotates the coordinates from the Global to the local coordinate system (note the inverse order of the subscripts):\n",
+ "\n",
+ "$$ \\mathbf{R_{lG}} = \\begin{bmatrix}\n",
+ "\\cos\\alpha & \\sin\\alpha \\\\\n",
+ "-\\sin\\alpha & \\cos\\alpha \n",
+ "\\end{bmatrix} $$\n",
+ "\n",
+ "If we use the direction cosines to calculate the rotation matrix, because the axes didn't change, the cosines are the same, only the order changes, now $\\mathbf{x, y}$ are the rows (outputs) and $\\mathbf{X, Y}$ are the columns (inputs):\n",
+ "\n",
+ "$$ \\mathbf{R_{lG}} = \\begin{bmatrix}\n",
+ "\\cos\\mathbf{X}x & \\cos\\mathbf{Y}x \\\\\n",
+ "\\cos\\mathbf{X}y & \\cos\\mathbf{Y}y \n",
+ "\\end{bmatrix} = \n",
+ "\\begin{bmatrix}\n",
+ "\\cos(\\alpha) & \\cos(90^o-\\alpha) \\\\\n",
+ "\\cos(90^o+\\alpha) & \\cos(\\alpha)\n",
+ "\\end{bmatrix} = \n",
+ "\\begin{bmatrix}\n",
+ "\\cos\\alpha & \\sin\\alpha \\\\\n",
+ "-\\sin\\alpha & \\cos\\alpha \n",
+ "\\end{bmatrix} $$\n",
+ "\n",
+ "And defining the versors of the axes in the Global coordinate system for a basis in terms of the local coordinate system would also produce this latter rotation matrix.\n",
+ "\n",
+ "The two sets of equations and matrices for the rotations from Global-to-local and local-to-Global coordinate systems are very similar, this is no coincidence. Each of the rotation matrices we deduced, $\\mathbf{R_{Gl}}$ and $\\mathbf{R_{lG}}$, perform the inverse operation in relation to the other. Each matrix is the inverse of the other. \n",
+ "\n",
+ "In other words, the relation between the two rotation matrices means it is equivalent to instead of rotating the local coordinate system by $\\alpha$ in relation to the Global coordinate system, to rotate the Global coordinate system by $-\\alpha$ in relation to the local coordinate system; remember that $\\cos(-\\alpha)=\\cos(\\alpha)$ and $\\sin(-\\alpha)=-\\sin(\\alpha)$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Rotation of a Vector\n",
+ "\n",
+ "We can also use the rotation matrix to rotate a vector by a given angle around an axis of the coordinate system as shown in the figure below. \n",
+ "<br>\n",
+ "<figure><img src='./../images/rotation2Dvector.png' alt='rotation 2D of a vector'/> <figcaption><center><i>Figure. Rotation of a position vector $\\mathbf{P}$ by an angle $\\alpha$ in the two-dimensional space.</i></center></figcaption> </figure>\n",
+ "\n",
+ "We will not prove that we use the same rotation matrix, but think that in this case the vector position rotates by the same angle instead of the coordinate system. The new coordinates of the vector position $\\mathbf{P'}$ rotated by an angle $\\alpha$ is simply the rotation matrix (for the angle $\\alpha$) multiplied by the coordinates of the vector position $\\mathbf{P}$:\n",
+ "**[See here for a review about matrix and its main properties](http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/Matrix.ipynb)**.\n",
+ "\n",
+ "A nice property of the rotation matrix is that its inverse is the transpose of the matrix (because the columns/rows are mutually orthogonal and have norm equal to one). \n",
+ "This property can be shown with the rotation matrices we deduced:\n",
+ "Using the inverse and the transpose mathematical operations, the coordinates at the local coordinate system given the coordinates at the Global coordinate system and the rotation matrix can be obtained by: \n",
+ "Where we referred the inverse of $\\mathbf{R_{Gl}}\\;(\\:\\mathbf{R_{Gl}^{-1}})$ as $\\mathbf{R_{lG}}$ (note the different order of the subscripts). \n",
+ "\n",
+ "Let's show this calculation in NumPy:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Rotation matrix (RGl):\n",
+ " [[ 0.7071 -0.7071]\n",
+ " [ 0.7071 0.7071]]\n",
+ "Position at the local coordinate system (Pl):\n",
+ " [[1]\n",
+ " [1]]\n",
+ "Position at the Global coordinate system (PG=RGl*Pl):\n",
+ " [[0. ]\n",
+ " [1.41]]\n",
+ "Position at the local coordinate system using the inverse of RGl (Pl=RlG*PG):\n",
+ "Pl = np.array([[1, 1]]).T # transpose the array for correct matrix multiplication\n",
+ "print('Position at the local coordinate system (Pl):\\n', Pl)\n",
+ "\n",
+ "PG = np.dot(RGl, Pl) # the function dot() is used for matrix multiplication with arrays\n",
+ "print('Position at the Global coordinate system (PG=RGl*Pl):\\n', np.around(PG,2))\n",
+ "\n",
+ "Pl = np.dot(RGl.T, PG)\n",
+ "print('Position at the local coordinate system using the inverse of RGl (Pl=RlG*PG):\\n', Pl)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**In summary, some of the properties of the rotation matrix are:** \n",
+ "1. The columns of the rotation matrix form a basis of (independent) unit vectors (versors) and the rows are also independent versors since the transpose of the rotation matrix is another rotation matrix. \n",
+ "2. The rotation matrix is orthogonal. There is no linear combination of one of the lines or columns of the matrix that would lead to the other row or column, i.e., the lines and columns of the rotation matrix are independent, orthogonal, to each other (this is property 1 rewritten). Because each row and column have norm equal to one, this matrix is also sometimes said to be orthonormal. \n",
+ "3. The determinant of the rotation matrix is equal to one (or equal to -1 if a left-hand coordinate system was used, but you should rarely use that). For instance, the determinant of the rotation matrix we deduced is $cos\\alpha cos\\alpha - sin\\alpha(-sin\\alpha)=1$.\n",
+ "4. The inverse of the rotation matrix is equals to its transpose.\n",
+ "\n",
+ "**On the different meanings of the rotation matrix:** \n",
+ "- It represents the coordinate transformation between the coordinates of a point expressed in two different coordinate systems. \n",
+ "- It describes the rotation between two coordinate systems. The columns are the direction cosines (versors) of the axes of the rotated coordinate system in relation to the other coordinate system and the rows are also direction cosines (versors) for the inverse rotation. \n",
+ "- It is an operator for the calculation of the rotation of a vector in a coordinate system.\n",
+ "- Rotation matrices provide a means of numerically representing rotations without appealing to angular specification.\n",
+ "\n",
+ "**Which matrix to use, from local to Global or Global to local?** \n",
+ "- A typical use of the transformation is in movement analysis, where there are the fixed Global (laboratory) coordinate system and the local (moving, e.g. anatomical) coordinate system attached to each body segment. Because the movement of the body segment is measured in the Global coordinate system, using cameras for example, and we want to reconstruct the coordinates of the markers at the anatomical coordinate system, we want the transformation leading from the Global coordinate system to the local coordinate system.\n",
+ "- Of course, if you have one matrix, it is simple to get the other; you just have to pay attention to use the right one."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Translation and rotation\n",
+ "\n",
+ "Consider now the case where the local coordinate system is translated and rotated in relation to the Global coordinate system and a point is described in both coordinate systems as illustrated in the figure below (once again, remember that this is equivalent to describing a translation and a rotation between two rigid bodies). \n",
+ "<br>\n",
+ "<figure><img src='./../images/transrot2D.png' alt='translation and rotation 2D'/> <figcaption><center><i>Figure. A point in two-dimensional space represented in two coordinate systems, with one system translated and rotated.</i></center></figcaption> </figure>\n",
+ "\n",
+ "The position of point $\\mathbf{P}$ originally described in the local coordinate system, but now described in the Global coordinate system in vector form is:\n",
+ "This means that we first *disrotate* the local coordinate system and then correct for the translation between the two coordinate systems. Note that we can't invert this order: the point position is expressed in the local coordinate system and we can't add this vector to another vector expressed in the Global coordinate system, first we have to convert the vectors to the same coordinate system.\n",
+ "\n",
+ "If now we want to find the position of a point at the local coordinate system given its position in the Global coordinate system, the rotation matrix and the translation vector, we have to invert the expression above:\n",
+ "The expression above indicates that to perform the inverse operation, to go from the Global to the local coordinate system, we first translate and then rotate the coordinate system."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Transformation matrix\n",
+ "\n",
+ "It is possible to combine the translation and rotation operations in only one matrix, called the transformation matrix (also referred as homogeneous transformation matrix):\n",
+ "However, because $\\mathbf{T_{Gl}}$ is not orthonormal when there is a translation, its inverse is not its transpose. Its inverse in matrix form is given by:\n",
+ "A typical scenario in motion analysis is to calculate the rotation matrix using the position of markers placed on the moving rigid body. With the markers' positions, we create a local basis, which by definition is the rotation matrix for the rigid body with respect to the Global (laboratory) coordinate system. To define a coordinate system using a basís, we also will need to define an origin. \n",
+ "\n",
+ "Let's see how to calculate a basis given the markers' positions. \n",
+ "Consider the markers at m1=[1,1]`, m2=[1,2] and m3=[-1,1] measured in the Global coordinate system as illustrated in the figure below: \n",
+ "<br>\n",
+ "<figure><img src='./../images/transrot2Db.png' alt='translation and rotation 2D'/> <figcaption><center><i>Figure. Three points in the two-dimensional space, two possible vectors given these points, and the corresponding basis.</i></center></figcaption> </figure>\n",
+ "\n",
+ "A possible local coordinate system with origin at the position of m1 is also illustrated in the figure above. Intentionally, the three markers were chosen to form orthogonal vectors. \n",
+ "The translation vector between the two coordinate system is:\n",
+ "\n",
+ "$$\\mathbf{L_{Gl}} = m_1 - [0,0] = [1,1]$$\n",
+ "\n",
+ "The vectors expressing the axes of the local coordinate system are:\n",
+ "Note that these two vectors do not form a basis yet because they are not unit vectors (in fact, only *y* is not a unit vector). Let's normalize these vectors:\n",
+ "Beware that the versors above are not exactly the same as the ones shown in the right plot of the last figure, the versors above if plotted will start at the origin of the coordinate system, not at [1,1] as shown in the figure.\n",
+ "\n",
+ "We could have done this calculation in NumPy (we will need to do that when dealing with real data later):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x = [0. 1.] , y = [-2. 0.] \n",
+ "ex= [0. 1.] , ey= [-1. 0.]\n"
+ ]
+ }
+ ],
+ "source": [
+ "m1 = np.array([1.,1.]) # marker 1\n",
+ "m2 = np.array([1.,2.]) # marker 2\n",
+ "m3 = np.array([-1.,1.]) # marker 3\n",
+ "\n",
+ "x = m2 - m1 # vector x\n",
+ "y = m3 - m1 # vector y\n",
+ "\n",
+ "vx = x/np.linalg.norm(x) # versor x\n",
+ "vy = y/np.linalg.norm(y) # verson y\n",
+ "\n",
+ "print(\"x =\", x, \", y =\", y, \"\\nex=\", vx, \", ey=\", vy)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now, both $\\mathbf{e}_x$ and $\\mathbf{e}_y$ are unit vectors (versors) and they are orthogonal, a basis can be formed with these two versors, and we can represent the rotation matrix using this basis (just place the versors of this basis as columns of the rotation matrix):\n",
+ "\n",
+ "$$ \\mathbf{R_{Gl}} = \\begin{bmatrix}\n",
+ "0 & -1 \\\\\n",
+ "1 & 0 \n",
+ "\\end{bmatrix} $$\n",
+ "\n",
+ "This rotation matrix makes sense because from the figure above we see that the local coordinate system we defined is rotated by 90$^o$ in relation to the Global coordinate system and if we use the general form for the rotation matrix:\n",
+ "\n",
+ "$$ \\mathbf{R} = \\begin{bmatrix}\n",
+ "\\cos\\alpha & -\\sin\\alpha \\\\\n",
+ "\\sin\\alpha & \\cos\\alpha \n",
+ "\\end{bmatrix} = \n",
+ "\\begin{bmatrix}\n",
+ "\\cos90^o & -\\sin90^o \\\\\n",
+ "\\sin90^o & \\cos90^o \n",
+ "\\end{bmatrix} =\n",
+ "\\begin{bmatrix}\n",
+ "0 & -1 \\\\\n",
+ "1 & 0 \n",
+ "\\end{bmatrix} $$\n",
+ "\n",
+ "So, the position of any point in the local coordinate system can be represented in the Global coordinate system by:\n",
+ "For example, the point $\\mathbf{P_l}=[1,1]$ has the following position at the Global coordinate system:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Translation vector:\n",
+ " [[1]\n",
+ " [1]]\n",
+ "Rotation matrix:\n",
+ " [[ 0 -1]\n",
+ " [ 1 0]]\n",
+ "Position at the local coordinate system:\n",
+ " [[1]\n",
+ " [1]]\n",
+ "Position at the Global coordinate system, PG = LGl + RGl*Pl:\n",
+ " [[0]\n",
+ " [2]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "LGl = np.array([[1, 1]]).T\n",
+ "print('Translation vector:\\n', LGl)\n",
+ "\n",
+ "RGl = np.array([[0, -1], [1, 0]])\n",
+ "print('Rotation matrix:\\n', RGl)\n",
+ "\n",
+ "Pl = np.array([[1, 1]]).T\n",
+ "print('Position at the local coordinate system:\\n', Pl)\n",
+ "\n",
+ "PG = LGl + np.dot(RGl, Pl)\n",
+ "print('Position at the Global coordinate system, PG = LGl + RGl*Pl:\\n', PG)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Determination of the unknown angle of rotation\n",
+ "\n",
+ "If we didn't know the angle of rotation between the two coordinate systems, which is the typical situation in motion analysis, we simply would equate one of the terms of the two-dimensional rotation matrix in its algebraic form to its correspondent value in the numerical rotation matrix we calculated.\n",
+ "\n",
+ "For instance, taking the first term of the rotation matrices above: $\\cos\\alpha = 0$ implies that $\\theta$ is 90$^o$ or 270$^o$, but combining with another matrix term, $\\sin\\alpha = 1$, implies that $\\alpha=90^o$. We can solve this problem in one step using the tangent $(\\sin\\alpha/\\cos\\alpha)$ function with two terms of the rotation matrix and calculating the angle with the `arctan2(y, x)` function:"
+ "And this procedure would be repeated for each segment and for each instant of the analyzed movement to find the rotation of each segment."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Joint angle as a sequence of rotations of adjacent segments\n",
+ "\n",
+ "In the notebook about [two-dimensional angular kinematics](http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/AngularKinematics2D.ipynb), we calculated segment and joint angles using simple trigonometric relations. We can also calculate these two-dimensional angles using what we learned here about the rotation matrix.\n",
+ "\n",
+ "The segment angle will be given by the matrix representing the rotation from the laboratory coordinate system (G) to a coordinate system attached to the segment and the joint angle will be given by the matrix representing the rotation from one segment coordinate system (l1) to the other segment coordinate system (l2). So, we have to calculate two basis now, one for each segment and the joint angle will be given by the product between the two rotation matrices. \n",
+ "\n",
+ "To define a two-dimensional basis, we need to calculate vectors perpendicular to each of these lines. Here is a way of doing that. First, let's find three non-collinear points for each basis:"
+ "Same result as obtained in [Angular kinematics in a plane (2D)](http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/AngularKinematics2D.ipynb). "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Kinematic chain in a plain (2D)\n",
+ "\n",
+ "The fact that we simply multiplied the rotation matrices to calculate the rotation matrix of one segment in relation to the other is powerful and can be generalized for any number of segments: given a serial kinematic chain with links 1, 2, ..., n and 0 is the base/laboratory, the rotation matrix between the base and last link is: $\\mathbf{R_{n,n-1}R_{n-1,n-2} \\dots R_{2,1}R_{1,0}}$, where each matrix in this product (calculated from right to left) is the rotation of one link with respect to the next one. \n",
+ "\n",
+ "For instance, consider a kinematic chain with two links, the link 1 is rotated by $\\alpha_1$ with respect to the base (0) and the link 2 is rotated by $\\alpha_2$ with respect to the link 1. \n",
+ "Using Sympy, the rotation matrices for link 2 w.r.t. link 1 $(R_{12})$ and for link 1 w.r.t. base 0 $(R_{01})$ are: "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from IPython.display import display, Math\n",
+ "from sympy import sin, cos, Matrix, simplify, latex, symbols\n",
+ "The typical use of all these concepts is in the three-dimensional motion analysis where we will have to deal with angles in different planes, which needs a special manipulation as we will see next."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problems\n",
+ "\n",
+ "1. A local coordinate system is rotated 30$^o$ clockwise in relation to the Global reference system. \n",
+ " A. Determine the matrices for rotating one coordinate system to another (two-dimensional). \n",
+ " B. What are the coordinates of the point [1, 1] (local coordinate system) at the global coordinate system? \n",
+ " C. And if this point is at the Global coordinate system and we want the coordinates at the local coordinate system? \n",
+ " D. Consider that the local coordinate system, besides the rotation is also translated by [2, 2]. What are the matrices for rotation, translation, and transformation from one coordinate system to another (two-dimensional)? \n",
+ " E. Repeat B and C considering this translation.\n",
+ " \n",
+ "2. Consider a local coordinate system U rotated 45$^o$ clockwise in relation to the Global reference system and another local coordinate system V rotated 45$^o$ clockwise in relation to the local reference system U. \n",
+ " A. Determine the rotation matrices of all possible transformations between the coordinate systems. \n",
+ " B. For the point [1, 1] in the coordinate system U, what are its coordinates in coordinate system V and in the Global coordinate system? \n",
+ " \n",
+ "3. Using the rotation matrix, deduce the new coordinates of a square figure with coordinates [0, 0], [1, 0], [1, 1], and [0, 1] when rotated by 0$^o$, 45$^o$, 90$^o$, 135$^o$, and 180$^o$ (always clockwise).\n",
+ " \n",
+ "4. Solve the problem 2 of [Angular kinematics in a plane (2D)](http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/AngularKinematics2D.ipynb) but now using the concept of two-dimensional transformations. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## References\n",
+ "\n",
+ "- Robertson G, Caldwell G, Hamill J, Kamen G (2013) [Research Methods in Biomechanics](http://books.google.com.br/books?id=gRn8AAAAQBAJ). 2nd Edition. Human Kinetics. \n",
+ "- Ruina A, Rudra P (2013) [Introduction to Statics and Dynamics](http://ruina.tam.cornell.edu/Book/index.html). Oxford University Press. \n",
+ "- Winter DA (2009) [Biomechanics and motor control of human movement](http://books.google.com.br/books?id=_bFHL08IWfwC). 4 ed. Hoboken, EUA: Wiley. \n",
+ "- Zatsiorsky VM (1997) [Kinematics of Human Motion](http://books.google.com.br/books/about/Kinematics_of_Human_Motion.html?id=Pql_xXdbrMcC&redir_esc=y). Champaign, Human Kinetics."
+ "Author: J. R. Johansson (robert@riken.jp), http://dml.riken.jp/~rob/\n"
+ "\n"
+ "The latest version of this [IPython notebook](http://ipython.org/ipython-doc/dev/interactive/htmlnotebook.html) lecture is available at [http://github.com/jrjohansson/qutip-lectures](http://github.com/jrjohansson/qutip-lectures).\n"
+ "\n"
+ "The other notebooks in this lecture series are indexed at [http://jrjohansson.github.com](http://jrjohansson.github.com)."
+ "Consider a single atom coupled to a single cavity mode, as illustrated in the figure below. If there atom excitation rate $\\Gamma$ exceeds the relaxation rate, a population inversion can occur in the atom, and if coupled to the cavity the atom can then act as a photon pump on the cavity."
+ "Here we evolve the system with the Lindblad master equation solver, and we request that the expectation values of the operators $a^\\dagger a$ and $\\sigma_+\\sigma_-$ are returned by the solver by passing the list `[a.dag()*a, sm.dag()*sm]` as the fifth argument to the solver."
+ "- using the online [IPython notebook slide viewer](https://slideviewer.herokuapp.com/) (some slides of the notebook might not be properly rendered).\n"
+ "\n"
+ "This and other related IPython notebooks can be found at the course github repository:\n"
+ "* [Random variable](http://en.wikipedia.org/wiki/Random_variable): a variable whose value is subject to variations due to __chance__. A random variable can take on a set of possible different values, each with an associated probability, in contrast to other mathematical variables.\n"
+ "\n"
+ "* [Probability distribution](http://en.wikipedia.org/wiki/Probability_distribution): mathematical function describing the possible values of a random variable and their associated probabilities.\n"
+ "\n"
+ "* [Probability density function (pdf)](http://en.wikipedia.org/wiki/Probability_density_function) of a __continuous random variable__ is a function that describes the relative likelihood for this random variable to take on a given value. \n"
+ " * The probability of the random variable falling within a particular range of values is given by the integral of this variable’s density over that range.\n"
+ " * The probability density function is nonnegative everywhere, and its integral over the entire space is equal to one.\n"
+ "The probability distribution of a random variable is often characterised by a small number of parameters, which also have a practical interpretation.\n"
+ "\n"
+ "* [Mean](http://en.wikipedia.org/wiki/Mean) (a.k.a expected value) refers to one measure of the central tendency either of a probability distribution or of the random variable characterized by that distribution.\n"
+ " * population mean: $\\mu = \\operatorname{E}[X]$.\n"
+ " * estimation of sample mean: $\\bar{x}$.\n"
+ "* [Standard deviation](http://en.wikipedia.org/wiki/Standard_deviation) measures the amount of variation or dispersion from the mean.\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n"
+ " warnings.warn(\"The 'normed' kwarg is deprecated, and has been \""
+ "In applied mathematics, [test functions](http://en.wikipedia.org/wiki/Test_functions_for_optimization), also known as artificial landscapes, are useful to evaluate characteristics of optimization algorithms, such as:\n"
+ "\n"
+ "* Velocity of convergence.\n"
+ "* Precision.\n"
+ "* Robustness.\n"
+ "* General performance.\n"
+ "\n"
+ "DEAP has a number of test problems already implemented. See http://deap.readthedocs.org/en/latest/api/benchmarks.html"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)"
+ ));
+ testImageResult(entry, 1);
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### Mutation of an evolution strategy individual according to its strategy attribute. \n"
+ "First the strategy is mutated according to an extended log normal rule, \n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)"
+ ));
+ testImageResult(entry, 1);
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "The plot (most probably) shows a \"dark blue\" ellipse as all individuals are overlapping. "
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "Let's see how the evolutionary process took place in animated form."
+ ));
+
+ qDebug() << "command entry 39";
+ testCommandEntry(entry, 39, QString::fromUtf8(
+ "from matplotlib import animation\n"
+ "from IPython.display import HTML"
+ ));
+
+ qDebug() << "command entry 40";
+ testCommandEntry(entry, 40, QString::fromUtf8(
+ "def animate(i):\n"
+ " 'Updates all plots to match frame _i_ of the animation.'\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)"
+ ));
+ testHtmlResult(entry, 1, QString::fromUtf8(
+ "<IPython.core.display.HTML object>"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "How the population progressed as the evolution proceeded?"
+ "* In an evolution strategy, new candidate solutions are sampled according to a multivariate normal distribution in the $\\mathbb{R}^n$. \n"
+ "* Recombination amounts to selecting a new mean value for the distribution. \n"
+ "* Mutation amounts to adding a random vector, a perturbation with zero mean. \n"
+ "* Pairwise dependencies between the variables in the distribution are represented by a covariance matrix. \n"
+ "\n"
+ "### The covariance matrix adaptation (CMA) is a method to update the covariance matrix of this distribution. \n"
+ "\n"
+ "> This is particularly useful, if the objective function $f()$ is ill-conditioned."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### CMA-ES features\n"
+ "\n"
+ "* Adaptation of the covariance matrix amounts to learning a second order model of the underlying objective function.\n"
+ "* This is similar to the approximation of the inverse Hessian matrix in the Quasi-Newton method in classical optimization. \n"
+ "* In contrast to most classical methods, fewer assumptions on the nature of the underlying objective function are made. \n"
+ "* *Only the ranking between candidate solutions is exploited* for learning the sample distribution and neither derivatives nor even the function values themselves are required by the method."
+ "/home/mmmm1998/.local/lib/python3.6/site-packages/deap/creator.py:141: RuntimeWarning: A class named 'Individual' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.\n"
+ " RuntimeWarning)"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "We will place our start point by hand at $(5,5)$."
+ "### OK, but wouldn't it be nice to have an animated plot of how CMA-ES progressed? \n"
+ "\n"
+ "* We need to do some coding to make this animation work.\n"
+ "* We are going to create a class named `PlotableStrategy` that inherits from `deap.cma.Strategy`. This class logs the features we need to make the plots as evolution takes place. That is, for every iteration we store:\n"
+ " * Current centroid and covariance ellipsoid.\n"
+ " * Updated centroid and covariance.\n"
+ " * Sampled individuals.\n"
+ " * Evolution path.\n"
+ " \n"
+ "_Note_: I think that DEAP's implementation of CMA-ES has the drawback of storing information that should be stored as part of \"individuals\". I leave this for an afternoon hack."
+ ));
+
+ qDebug() << "command entry 50";
+ testCommandEntry(entry, 50, QString::fromUtf8(
+ "from math import sqrt, log, exp\n"
+ "class PlotableStrategy(cma.Strategy):\n"
+ " \"\"\"This is a modification of deap.cma.Strategy class.\n"
+ " We store the execution data in order to plot it.\n"
+ " **Note:** This class should not be used for other uses than\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ "/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:960: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride', 'linewidth'\n"
+ " s)"
+ ));
+ testHtmlResult(entry, 1, QString::fromUtf8(
+ "<IPython.core.display.HTML object>"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "* Current centroid and covariance: **red**.\n"
+ "* Updated centroid and covariance: **green**. \n"
+ "* Sampled individuals: **shades of gray representing their corresponding weight**.\n"
+ "* Evolution path: **blue line starting in (0,0)**. "
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Homework\n"
+ "\n"
+ "1. Make an animated plot with the covariance update process. You can rely on the notebook of the previous demonstration class.\n"
+ "2. Compare ES, CMA-ES and a genetic algortihm.\n"
+ "2. How do you think that evolutionary strategies and CMA-ES should be modified in order to cope with combinatorial problems?\n"
+ "3. How can evolution strategies be improved?\n"
+ " This work is licensed under a [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License](http://creativecommons.org/licenses/by-nc-sa/4.0/).\n"
+ "In order to study population genetics, we first need a model of a population. And even before that, we need to define what we mean by *population*. Populations can be defined on many levels and with many diffferent criteria. For our purposes, we will simply say that a population is a set of individuals sharing a common environment. And because this is population *genetics* we can think of individuals as entities comprising of specific genes or chromosomes. \n"
+ "\n"
+ "So where do we get a population from? As you may have discussed in previous workshops, there are very large datasets containing sequencing information from different populations. So we could download one of these datasets and perform some analysis on it. But I find this can be dry and tedious. So why download data when we can simply create our own?\n"
+ "\n"
+ "In this workshop we're going to be creating and studying our own \"artificial\" populations to illustrate some important population genetics concepts and methodologies. Not only will this help you learn population genetics, but you will get a lot more programming practice than if we were to simply parse data files and go from there. \n"
+ "\n"
+ "More specifically, we're going to build our own RNA world.\n"
+ "\n"
+ "As you may know, RNA is widely thought to be the first self replicating life-form to arise x billion years ago. One of the strongest arguments for this theory is that RNA is able to carry information in its nucleotides like DNA, and like protein, it is able to adopt higher order structures to catalyze reactions, such as self replication. So it is likely, and there is growing evidence that this is the case, that the first form of replicating life was RNA. And because of this dual property of RNA as an information vessel as well as a structural/functional element we can use RNA molecules to build very nice population models. \n"
+ "\n"
+ "So in this notebook, I'll be walking you through building genetic populations, simulating their evolution, and using statistics and other mathematical tools for understanding key properties of populations.\n"
+ "\n"
+ "### Building an RNA population\n"
+ "\n"
+ "As we saw earlier, RNA has the nice property of posessing a strong mapping between information carrying (sequence) and function (structure). This is analogous to what is known in evolutionary terms as a genotype and a phenotype. With these properties, we have everything we need to model a population, and simulate its evolution.\n"
+ "\n"
+ "#### RNA sequence-structure\n"
+ "\n"
+ "We can think of the genotype as a sequence $s$ consisting of letters/nucleotides from the alphabet $\\{U,A,C,G\\}$. The corresponding phenotype $\\omega$ is the secondary structure of $s$ which can be thought of as a pairing between nucleotides in the primary sequence that give rise to a 2D architecture. Because it has been shown that the function of many biomolecules, including RNA, is driven by structure this gives us a good proxy for phenotype. \n"
+ "\n"
+ "Below is an example of what an RNA secondary structure, or pairing, looks like."
+ "As you can see, unparied positions are forming loop-like structures, and paired positions are forming stem-like structures. It is this spatial arrangement of nucleotides that drives RNA's function. Therefore, another sequence that adopts a similar shape, is likely to behave in a similar manner. Another thing to notice is that, although in reality this is often not the case, in general we only allow pairs between $\\{C,G\\}$ and $\\{A, U\\}$ nucleotides, most modern approaches allow for non-canonical pairings and you will find some examples of this in the above structure.\n"
+ "\n"
+ "*How do we go from a sequence to a structure?*\n"
+ "\n"
+ "So a secondary structure is just a list of pairings between positions. How do we get the optimal pairing?\n"
+ "\n"
+ "The algorithm we're going to be using in our simulations is known as the Nussinov Algorithm. The Nussinov algorithm is one of the first and simplest attempts at predicting RNA structure. Because bonds tend to stabilize RNA, the algorithm tries to maximize the number of pairs in the structure and return that as its solution. Current approaches achieve more accurate solutions by using energy models based one experimental values to then obtain a structure that minimizes free energy. But since we're not really concerned with the accuracy of our predictions, Nussinov is a good entry point. Furthermore, the main algorithmic concepts are the same between Nussinov and state of the art RNA structure prediction algorithms. I implemented the algorithm in a separate file called `fold.py` that we can import and use its functions. I'm not going to go into detail here on how the algorithm works because it is beyond the scope of this workshop but there is a bonus exercise at the end if you're curious.\n"
+ "\n"
+ "You can predict a secondary structure by calling `nussinov()` with a sequence string and it will return a tuple in the form `(structure, pairs)`."
+ "You will see a funny dot-bracket string in the output. This is a representation of the structure of an RNA. Quite simply, a matching parir of parentheses (open and close) correspond to the nucleotides at those positions being paried. Whereas, a dot means that that position is unpaired in the structure. Feel free to play around with the input sequence to get a better understanding of the notation.\n"
+ "\n"
+ "So that's enough about RNA structure prediction. Let's move on to building our populations.\n"
+ "\n"
+ "### Fitness of a sequence: Target Structure\n"
+ "\n"
+ "Now that we have a good way of getting a phenotype (secondary structure), we need a way to evaluate the fitness of that phenotype. If we think in real life terms, fitness is the ability of a genotype to replicate into the next generation. If you have a gene carrying a mutation that causes some kind of disease, your fitness is decreased and you have a lower chance of contributing offspring to the next generation. On a molecular level the same concept applies. A molecule needs to accomplish a certain function, i.e. bind to some other molecule or send some kind of signal. And as we've seen before, the most important factor that determines how well it can carry out this function is its structure. So we can imagine that a certain structure, we can call this a 'target' structure, is required in order to accomplish a certain function. So a sequence that folds correctly to a target structure is seen as having a greater fitness than one that does not. Since we've encoded structures as simple dot-bracket strings, we can easily compare structures and thus evaluate the fitness between a given structure and the target, or 'correct' structure. \n"
+ "\n"
+ "There are many ways to compare structures $w_{1}$ and $w_{2}$, but we're going to use one of the simplest ways, which is base-pair distance. This is just the number of pairs in $w_{1}$ that are not in $w_{2}$. Again, this is beyond the scope of this workshop so I'll just give you the code for it and if you would like to know more you can ask me."
+ "#give the list of pairs to bp_distance and see what the distance is.\n"
+ "print(bp_distance(w1[-1], w2[-1]))"
+ ));
+ testTextResult(entry, 0, QString::fromUtf8(
+ "('((....))', [(0, 7), (1, 6)])\n"
+ "('.(....).', [(1, 6)])\n"
+ "1"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Defining a cell: a little bit of Object Oriented Programming (OOP)\n"
+ "\n"
+ "Since we're going to be playing aroudn with sequences and structures and fitness values a lot, it's best to package it all nicely into an object. As you'll have seen with Vlad, objects are just a nice way of grouping data into an easily accessible form. \n"
+ "\n"
+ "We're trying to simulate evolution on a very simple kind of organism, or cell. It contains two copies of a RNA gene, each with a corresponding structure. "
+ "Now we've defined a 'Cell'. Since a population is a collection of individuals our populations will naturally consist of **lists** of 'Cell' objects, each with their own sequences. Here we initialize all the Cells with random sequences and add them to the 'population' list."
+ ));
+
+ qDebug() << "command entry 5";
+ testCommandEntry(entry, 5, QString::fromUtf8(
+ "### 5\n"
+ "import random\n"
+ "\n"
+ "def populate(target, pop_size=100):\n"
+ " \n"
+ " population = []\n"
+ "\n"
+ " for i in range(pop_size):\n"
+ " #get a random sequence to start with\n"
+ " sequence = \"\".join([random.choice(\"AUCG\") for _ in range(len(target))])\n"
+ " #use nussinov to get the secondary structure for the sequence\n"
+ " structure = nussinov(sequence)\n"
+ " #add a new Cell object to the population list\n"
+ "Now that we can store populatoins of cells, we need a way to evaluate the fitness of a given Cell. Recall that a Cell is simply an object that contains two RNA sequences (think of it as two copies of a gene on each chromosome). \n"
+ "\n"
+ "So we simply need to loop through each Cell in a population and compute base pair distance to the target structure. However, simply using base-pair distance is not a very good measure of fitness. There are two reasons for this: \n"
+ "\n"
+ "1. We want fitness to represent a *probability* that a cell will reproduce, and base pair distance is an integer.\n"
+ "2. We want this probability to be a *relative* measure. That is, we want to be the fitness to be proportional to how good a cell is with respect to all others in the population. This touches on an important principle in evolution where we only need to be 'better' than the competition and not good in some absolute measure. For example, if you and I are being chased by a bear. In order to survive, I only need to be faster than you, and not necessarily some absolute level of fitness.\n"
+ "\n"
+ "In order to get a probability (number between 0 and 1) we use the following equation to define the fitness of a structure $\\omega$ on a target structure $T$:\n"
+ "As you can see, it's a very simple function that evaluates to 1 (highest fitness) if the base pair distance is 0, and decreases as the structures get further and further away from the target. I didn't include the $N$ in the plotting as it will be a bit more annoying to compute, but it is simply a scaling factor so the shape and main idea won't be different.\n"
+ "\n"
+ "Now we can use this function to get a fitness value for each Cell in our population."
+ "Evolution would go nowhere without random mutations. While mutations are technically just random errors in the copying of genetic material, they are essential in the process of evolution. This is because they introduce novel diversity to populatons, which with a low frequency can be beneficial. And when a beneficial mutation arises (i.e. a mutation that increases fitness, or replication probability) it quickly takes over the population and the populatioin as a whole has a higher fitness.\n"
+ "\n"
+ "Implementing mutations in our model will be quite straightforward. Since mutations happen at the genotype/sequence level, we simply have to iterate through our strings of nucleotides (sequences) and randomly introduce changes."
+ "The final process in this evolution model is selection. Once you have populations with a diverse range of fitnesses, we need to select the fittest individuals and let them replicate and contribute offspring to the next generation. In real populations this is just the process of reproduction. If you're fit enough you will be likely to reproduce more than another individual who is not as well suited to the environment.\n"
+ "\n"
+ "In order to represent this process in our model, we will use the fitness values that we assigned to each Cell earlier and use that to select replicating Cells. This is equivalent to sampling from a population with the sampling being weighted by the fitness of each Cell. Thankfully, `numpy.random.choice` comes to the rescue here. Once we have sampled enough Cells to build our next generation, we introduce mutations and compute the fitness values of the new generation."
+ "Here we simply store some statistics (in a dictionary) on the population at each generation such as the average base pair distance and the average fitness of the populations. No coding to do here, it's not a very interesting function but feel free to give it a look."
+ ));
+
+ qDebug() << "command entry 12";
+ testCommandEntry(entry, 12, QString::fromUtf8(
+ "def record_stats(pop, population_stats):\n"
+ " \"\"\"\n"
+ " Takes a population list and a dictionary and updates it with stats on the population.\n"
+ " \"\"\"\n"
+ " generation_bp_distance_1 = [rna.bp_distance_1 for rna in pop]\n"
+ " generation_bp_distance_2 = [rna.bp_distance_2 for rna in pop]\n"
+ "Let's see if it actually worked by plotting the average base pair distance as a function of generations for both genes in each cell. We should expect a gradual decrease as the populations get closer to the target structure."
+ " Plot base pair distance for each chromosome over generations.\n"
+ " \"\"\"\n"
+ " for m in ['mean_bp_distance_1', 'mean_bp_distance_2']:\n"
+ " plt.plot(pops_stats[m], label=m)\n"
+ " plt.legend()\n"
+ " plt.xlabel(\"Generations\")\n"
+ " plt.ylabel(\"Mean Base Pair Distance\")\n"
+ " \n"
+ "evo_plot(pops_stats)"
+ ));
+ testImageResult(entry, 0);
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "You should see a nice drop in base pair distance! Another way of visualizing this is by plotting a histogram of the base pair distance of all Cells in the initial population versus the final population."
+ " Plots histograms of base pair distance in initial and final populations.\n"
+ " \"\"\"\n"
+ " #plot bp_distance_1 for rnas in first population\n"
+ " g = sns.distplot([rna.bp_distance_1 for rna in pops[0]], label='initial population')\n"
+ " #plot bp_distance_1 for rnas in first population\n"
+ " g = sns.distplot([rna.bp_distance_1 for rna in pops[-1]], label='final population')\n"
+ " g.set(xlabel='Mean Base Pair Distance')\n"
+ " g.legend()\n"
+ "bp_distance_distributions(pops)"
+ ));
+ testImageResult(entry, 0);
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Studying our evolved sequences with some Population Genetics tools\n"
+ "\n"
+ "Now that we've generated some sequences, we can analyze them!\n"
+ "\n"
+ "So after several rounds of selection, what do we get? We have a bunch of different sequences. We would like a way to characterize this diversity. One important tool for doing this is by making what is known as phylogenetic trees. \n"
+ "\n"
+ "Phylogenetic trees tell us about which groups of similar sequences are present and how they are likely related in evolutionary time. \n"
+ "\n"
+ "There are several ways of building phylogenetic trees using BioPython. Here we will go over one type and I'll leave another one as an exercise.\n"
+ "\n"
+ "### UPGMA (Unweighted Pair Group Method with Arithmetic Means)\n"
+ "\n"
+ "This is basically a clustering method based on the distance (or number of differences) between every pair of sequences. It assumes that sequences that are more similar are more likely to be related than the other way around. \n"
+ "\n"
+ "For $N$ sequences, the algorithm builds an $NxN$ matrix that stores the distance between each sequence to every other sequence. The algorithm goes through this matrix and finds the pair of sequences that is most similar and merges it into a 'cluster' or in tree terms, connects them to a common node. This process is repeated until all the sequences have been assigned to a group. Refer to the wikipedia article on [UPGMA](https://en.wikipedia.org/wiki/UPGMA) for a more detailed explanation. \n"
+ ""
+ ));
+
+ qDebug() << "command entry 18";
+ testCommandEntry(entry, 18, QString::fromUtf8(
+ "from Bio import SeqIO\n"
+ "from Bio.Seq import Seq\n"
+ "from Bio.SeqRecord import SeqRecord\n"
+ "from Bio import AlignIO"
+ ));
+
+ qDebug() << "command entry 19";
+ testCommandEntry(entry, 19, QString::fromUtf8(
+ "sequences = []\n"
+ "#let's take the first 10 sequences of our population to keep things simple\n"
+ "for seq in pops[-1][:10]:\n"
+ " #store each sequence in the sequences list as a SeqRecord object\n"
+ "The UPGMA algorithm requires a `MultipleSeqAlignment` object to build the distance matrix. So now that we have the `seq.fasta` file, we can give it to an online multiple sequence alignment tool. We can do this through BioPython but it requires some installation and setup so we will skip that for now. Go to the [MUSCLE Web Server](http://www.ebi.ac.uk/Tools/msa/muscle/) and give it the `seq.fasta` file. It will take a few seconds and it will give you an alignment and click *Download Alignment File*, copy paste the whole thing to a new file called `aln.clustal`. This is the alignment we will use to build our tree."
+ ));
+
+ qDebug() << "command entry 21";
+ testCommandEntry(entry, 21, QString::fromUtf8(
+ "#open the alignmnent file\n"
+ "with open(\"aln.clustal\", \"r\") as aln:\n"
+ " #use AlignIO to read the alignment file in 'clustal' format\n"
+ "The populations we generated evolved asexually. This means that individuals do not mate or exchange genetic information. So to make our simulation a bit more interesting let's let the Cells mate. This is going to require a few small changes in the `selection()` function. Previously, when we selected sequences to go into the next generation we just let them provide one offspring which was a copy of itself and introduced mutations. Now instead of choosing one Cell at a time, we will randomly choose two 'parents' that will mate. When they mate, each parent will contribute one of its chromosomes to the child. We'll repeat this process until we have filled the next generation."
+ "When we are presented with data from a population we don't know much about. It is often useful to try to learn whether there are any evolutionary or behavioural influences that are shaping population dynamics. This could be in the form of selective pressure, mating preference, genetic drift, mutations, gene flow, etc. So in order to detect if something like this is happening we need to develop a test. This is where Hardy Weinberg comes in. \n"
+ "\n"
+ "The Hardy Weinberg equilibrium states that \"allele and genotype frequencies remain constant in the absence of other evolutionary influences. (such as the ones we mentioned above)\" - Wikipedia.\n"
+ "\n"
+ "So if we can measure allele/genotype frequencies (which we can do because we have sequences), we can see whether the HW principle holds true. If it does not, then we can do more digging to see what could be happening to shift populations away from equilibrium.\n"
+ "\n"
+ "In order to do this we need to define an 'allele'. An allele (for our purproses) will be a locus (position in a sequence) that can take one of two states, a *reference* state or an *alternate* state. For example, we can look at locus number **5** (position 5 in our RNA sequences) and call reference **C**, and alternate **G**. If we are in HW we can predict the frequency of each allele in our population.\n"
+ "\n"
+ "To simplify our notation we will call the alternate allele *A* and the reference allele *a*. We can write the probability of each allele as $p_{A} + p_{a} = 1$. Since we are dealing with diploid populations, each individual will have two copies of each locus so it can be $p_{AA}, p{Aa}, p{aA}, p{aa}$. By simple probability laws we can get an expression for the probability of each genotype based on the probabilities of the single loci $p_{a}$ and $p_{A}$.\n"
+ "\n"
+ "$$p_{aa}\\simeq p_{a}^2$$\n"
+ "\n"
+ "$$p_{AA}\\simeq p_{A}^2$$\n"
+ "\n"
+ "$$p_{Aa,~aA} \\simeq 2 p_{a} p_{A}.$$\n"
+ "\n"
+ "Since it is hard to know what the true probability of observing either $p_{a}$ and $p_{A}$ we can estimate this probability from our data as follows:\n"
+ "Now we have expected values, and observed values. We need a test to determine whether we have a significant departure from the hypothesis of Hardy Weinberg equilibrium. The statistical test that is commonly used is known as the $\\chi^{2}$ test. If you take a look at the equation you'll see that the statistic simply takes the squared difference between our observed value and the expected value (divided by expected) and sums this for each possible genotype. The reason we take the squared difference is because we want to deal only with positive values, hence the name $\\chi^{2}$.\n"
+ "The first thing we need to do is get alleles from our sequence data. This boils down to going through each sequence at the position of interest and counting the number of $AA$, $Aa$, $aa$ we get.\n"
+ "\n"
+ "\n"
+ "\\** the sections on Hardy Weinberg and F-statistics are adapted from Simon Gravel's HGEN 661 Notes"
+ ));
+
+ qDebug() << "command entry 29";
+ testCommandEntry(entry, 29, QString::fromUtf8(
+ "def allele_finder(pop, locus, ref, alt):\n"
+ " genotypes = []\n"
+ " for p in pop:\n"
+ " #get the nucleotide at the locus from the first chromosome \n"
+ " locus_1 = p.sequence_1[locus].upper()\n"
+ " #same for the second\n"
+ " locus_2 = p.sequence_2[locus].upper()\n"
+ " \n"
+ " #check that it is either ref or alt, we don't care about other alleles for now.\n"
+ " if locus_1 in (ref, alt) and locus_2 in (ref, alt):\n"
+ " #if the alelle is ref, store a value of 1 in allele_1, and 0 otherwise\n"
+ " allele_1 = int(locus_1 == ref)\n"
+ " #same for the second allele\n"
+ " allele_2 = int(locus_2 == ref)\n"
+ " \n"
+ " #add allele to our list of alleles as a tuple. \n"
+ "Now that we have alleles represented in the right form, we can see if our population is at Hardy Weinberg equilibrium using the $\\chi_{2}$ test and the equations above."
+ "Can we say that our population is at equilibrium? Can you find parameters for `evolution_with_mating()` that will give us populations outside of the HW equilibrium?"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## A brief interlude on the p-value\n"
+ "\n"
+ "Let's take a minute to understand what the p-value means. The p-value is a probability. Specifically, it is the probability of observing a value equal to or more extreme than that our statistic given the test distribution. So in our case, it is the probability of observing a $X^2$ greater than or equal to the one the test gives us under a $\\chi^2$ distribution. When this value is very small, it suggests that it is unlikely that we are sampling from our assumed 'null' distribution and that some other alternate distribution is the true distribution. So a low p-value here would be evidence against the neutral Hardy Weinberg model and would suggest that our population is experiencing some influences such as mating preference, selection, mutation etc.\n"
+ "\n"
+ "A lot of research bases its conclusions solely on p-value and it is important to be very wary of this bad practice. It has become a bad convention that people say a p-value lower than some arbitrary threshold means one's findings are significant. However, very often the p-value does not give us the whole story and we need to know about things like sample size, size of impact, reproducibility, power of the test, etc. (check this out [American Statistical Association statement on p-values](http://www.nature.com/news/statisticians-issue-warning-over-misuse-of-p-values-1.19503), [p-hacking](http://fivethirtyeight.com/features/science-isnt-broken/#part1), and [this](http://allendowney.blogspot.ca/2016/06/there-is-still-only-one-test.html))\n"
+ "\n"
+ "Let's just visualize this very quickly using the $\\chi^{2}_{1}$ distribution. You will see that the p-value corresponds to the shaded red area under the curve. That area is the probability of observing a value as extreme or more than the one we found. When that is a very small area, we can be more confident that our assumption of HW is false."
+ "The last topic we'll cover is F-statistics. \n"
+ "\n"
+ "Once we find that our population strays from the HW condition we can begin to ask why that is the case. Often this deviation from the expected allele frequencies under HW is due to mating preference. Hardy Weinberg assumes that all individuals in a population have an equal probability of mating with any other individual (random mating). However, when certain individuals prefer to mate with specific others (in real populations this can be due to culture, race, geographic barriers, etc.), you get what is known as population structure. Population structure means that we begin to see *sub-populations* within our total population where individuals prefer to mate within their sub-population. This biased mating will result in a higher number of homozygotes than we would expect under Hardy-Weinberg equilibrium. Simply because mating preferences will tend to drive populations toward similar genotypes. So if this is the case, and no other factors are biasing allele dynamics, within sub-populations we should have Hardy-Weinberg like conditions. \n"
+ "\n"
+ "For example, if Raptors fans prefer to mate with other Raptors fans, then when we consider only Raptors fans, we should observe random mating. Simply because if the mating preference criterion is 'being a Raptor's fan' then any Raptor's fan will be equally likely to mate with any other Raptor's fan so we have Hardy Weinberg again.\n"
+ "\n"
+ "Let's express this in quantities we can measure.\n"
+ "\n"
+ "From before we calculated the observed and expected number of heterozygotes in a population. Let's call these $\\hat H$ and $H_{IT}$ respectively. $\\hat H$ is just the count of heterozygotes, and $H_{IT}$ is the same as the expected number of heterozygotes we calculated earlier.\n"
+ "\n"
+ "We define a quantity $e_{IT}$ as a measure of the 'excess heterozygosity' in the population when we consider all individuals $I$ in the total population $T$. $e_{IT} > 1$ when we have more heterozygotes than we expect under HW. And $0 < e_{IT} < 1$ if we have less heterozygotes than we would expect under HW.\n"
+ "\n"
+ "\n"
+ "$$e_{IT}=\\frac{\\mbox{observed proportion of hets}}{\\mbox{expected proportion of hets}}=\\frac{ H_{obs}}{H_{IT}}$$\n"
+ "\n"
+ "We use $e_{IT}$ to define the statistic $F_{IT}$\n"
+ "\n"
+ "$$F_{IT}=1-e_{IT}$$\n"
+ "\n"
+ "So $F_{IT} > 0$ when we have a lack of heterozygotes and $F_{IT} < 0$ when we have an excess of heterozygotes. $F_{IT} = 0$ under random mating.\n"
+ "\n"
+ "When we have a subpropulation $S$ we can calculate the equivalent quantity but instead of considering heterozygosity in the whole population we only take a sub-population into account.\n"
+ "\n"
+ "$$e_{IS} = \\frac{H_{obs}}{H_{IS}}$$\n"
+ "\n"
+ "And lastly, we have $F_{ST}$. This one is not as intuitive to derive so I'm not including the derivation here. But basically it measure the excess heterozygosity in the total population due to the presence of two subpopulations with allele frequencies $p_{1}$ and $p_{2}$.\n"
+ "Let's get some data for our F-tests. First we need to evolve two populations indepenently of each other, to simulate isolated mating. Then to simulate the total population we combine the two sub-populations. We then use our `allele_finder()` function to get all the alleles, and the `hardy_weinberg_chi_2_test()` function to get our expected and observed counts. Finally we plug those into the `f_statistics()` function."
+ "Try playing with different evolution parameters and see the effect on the different F-statistics. This workshop is a work in progress so there may be some biases in our simulation scheme that can make for come confusing F-statistics. If you come up with anything interesting I would love to know about it."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Exercises / Extra Material\n"
+ "\n"
+ "### Programming Exercises\n"
+ "\n"
+ "i. *Heatmap of mutation rates vs. population sizes.* (short) \n"
+ "\n"
+ "Make a heatmap that plots the base pair distance of the average base pair distance of the population at generation `-1` for mutation rates $\\mu = \\{0, 0.001, 0.01, 0.1, 0.5\\}$ and population sizes $N=\\{10, 100, 1000, 10000\\}$. The resulting heatmap will be `5x4` dimensions. You may choose how many generations to evolve your populations, just plot the last one in the heatmap."
+ "#xticks/yticks takes a list of numbers that specify the position of the ticks and a list with the tick labels\n"
+ "plt.xticks([i + .5 for i in range(len(population_sizes))], population_sizes)\n"
+ "plt.ylabel(\"Mutation Rate\")\n"
+ "plt.yticks([i + .5 for i in range(len(mutation_rates))], mutation_rates)"
+ ));
+ testTextResult(entry, 0, QString::fromLatin1(
+ "([<matplotlib.axis.YTick at 0x7fa02ba18ac8>,\n"
+ " <matplotlib.axis.YTick at 0x7fa02ba183c8>,\n"
+ " <matplotlib.axis.YTick at 0x7fa02ba59a58>,\n"
+ " <matplotlib.axis.YTick at 0x7fa02acbe978>,\n"
+ " <matplotlib.axis.YTick at 0x7fa02acc5048>],\n"
+ " <a list of 5 Text yticklabel objects>)"
+ ));
+ testImageResult(entry, 1);
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "ii. *Introduce mating preferences within a population.* (medium length) \n"
+ "\n"
+ "Modify the `selection_with_mating()` function to allow for mating preferences within a population. In our example above we were just running two independent simulations to study barriers to gene flow. But now you will implement mating preferences within a single simulation. Your function will assign each Cell a new attribute called `self.preference` which will take a string value denoting the mating type the current cell prefers to mate with. For example we can have a population with three mating types: $\\{A, B, C\\}$. Your function will randomly assign preferences to each cell in the initial population. We will define a preference between types $A$ and $B$ as the probability that two cells of those given types will mate if selected. \n"
+ "\n"
+ "$$\n"
+ "preferences(A,B,C) = \n"
+ "\\begin{bmatrix}\n"
+ " 0.7 & 0.1 & 0.2 \\\\\n"
+ " 0.1 & 0.9 & 0 \\\\\n"
+ " 0.2 & 0 & 0.8 \\\\\n"
+ "\\end{bmatrix}\n"
+ "$$\n"
+ "\n"
+ "Once you selected two potential parents for mating (as we did earlier) you will use the matrix to evaluate whether or not the two parents will mate and contribute an offspring to the next generation. "
+ "### Population Genetics / Bioinformatics Exercises\n"
+ "\n"
+ "*Exercise 1. Make a tree using maximum parsimony.*\n"
+ "\n"
+ "We saw how to make trees using a distance score. Another popular method is known as the maximum parsimony approach. I won't go into too much detail on this since we are short on time, but I will give a quick intro and we'll look at how ot make a tree using maximum parsimony.\n"
+ "\n"
+ "This approach is based on the principle of parsimony, which states that the simplest explanation for our data is the most likely to be true. So given an alignment, we assume that the best tree is the one that minimizes the number of changes, or mutations. This is often a reasonable assumption to make since mutation rates in real populations are generally low, and things like back-mutations (e.g. A --> C --> A) are unlikely. Computing the tree that that maximizes parsimony directly is a difficult task, but evaluating the parsimony score of a tree given the tree is easy. So this approach basically generates many random trees for the data and scores them based on parsimony keeping the most parsimonious tree. Take a look at [the biopython manual to work through this example](http://biopython.org/wiki/Phylo), and [this one](http://biopython.org/DIST/docs/api/Bio.Phylo.TreeConstruction.ParsimonyTreeConstructor-class.html).\n"
+ "\n"
+ "Since we already have an alignment (`aln.clustal`) we will just re-use it and make a maximum parsimony tree instead. "
+ "We just saw two methods of growing phylogenetic trees given an alignment. However, as we saw with the maximum parsimony approach, there can be many different trees for a single data set. How do we know our tree is a good representation of the data? By 'good' here we will instead use the word 'robust'. Is the tree we use too sensitive to the particularities of the data we gave it? If we make a small change in the sequence will we get a very different tree? Normally these problems would be addressed by re-sampling and seeing if we obtain similar results. But we can't really re-sample evolution. It happened once and we can't make it happen again. So we use something called *bootstrapping* which is a technique often used in statistics where instead of generating new data, you re-sample from your present data.\n"
+ "\n"
+ "So we have a multiple sequence alignment with $M$ sequences (rows) each with sequences of length $N$ nucleotides (columns). For each row, we can randomly sample $N$ nucleotides with replacement to make a new 'bootstrapped' sequence also of length $N$. Think of it as a kind of shuffling of the data. This gives us a whole new alignment that we can again use to make a new tree.\n"
+ "\n"
+ "This process is repeated many times to obtain many trees. The differences in topology (shape/structure) of the trees we obtained are assessed. If after this shuffling/perturbations we still get similar enough looking trees we can say that our final tree is robust to small changes in the data. ([some more reading on this](http://projecteuclid.org/download/pdf_1/euclid.ss/1063994979))\n"
+ "\n"
+ "Let's run a small example of this using the bootstrapping functions in `BioPython`."
+ "Similarly to the $\\chi^{2}$ test we saw for testing deviations from HW equilibrium, we can use a T-test to compare differences in means between two independent samples. We can use this to revisit a the first programming question in the exercsies section. Does mutation rate and population size have an effect on the fitness of populations? We can translate this question to, is there a difference in the mean base pair distance between populations under different mutation and population size regimes?\n"
+ "\n"
+ "Scipy has a very useful function that implements the T-test called `scipy.stats.ttest_ind`. Run two independent simulations (with different mutation rates) and compute the difference in mean bp distance between the two at their final generation. Store the populations in two different variables. Give a list of `bp_distance_1` values for each memeber of the population to `ttest_ind()`. \n"
+ "\n"
+ "Make sure to read teh `ttest_ind()` documentation, particularly about the argumetn `equal_var`. What should we set it to?"
+ "1. *Nussinov Algorithm (Only try this if you are feeling brave and are done with the other exercises or are interested in getting a taste of Computer Science. It is beyond the scope of this workshop.)*\n"
+ "\n"
+ "There are several approaches for solving this problem, we will look at the simplest one here which is known as the Nussinov Algorithm. This algorithm is a popular example of a class of algorithms know as dynamic programming algorithms. The main idea behind these algorithms is that we can break down the problem into many subproblems which are easier to compute than the full problem. Once we have obtained the solution for the subproblems, we can retrieve the solution to the full problem by doing something called a backtrace (more on the backtrace later). \n"
+ "\n"
+ "Here, the problem is obtaining the optimal pairing on a string of nucleotides. In order to know how good our structure is, we assign a score to it. One possible scoring scheme could be adding 1 to the score per paired set of nucleotides, and 0 otherwise. So in other words, we want a pairing that will give us the highest possible score. We can write this quantity as $OPT(i, j)$ where $i$ and $j$ are the indices of the sequence between which we obtain the pairing score. Our algorithm is therefore going to compute a folding score for all substrings bound by $i$ and $j$ and store the value in what is known as a dynamic programming table. Our dynamic programming table will be a $N$ x $N$ array where $N$ is the length of our sequence. So now that we have a way of measuring how good a structure is, we need a way to evaluate scores given a subsequence. To do this, we set some rules on the structure of an RNA sequence:\n"
+ "\n"
+ "\n"
+ "If $i$ and $j$ form a pair:\n"
+ "1. The pair $i$ and $j$ must form a valid watson-crick pair.\n"
+ "2. $i < j-4$. This ensures that bonding is not happening between positions that are too close to each other, which would produce steric clashes.\n"
+ "3. If pair $(i,j)$ and $(k, l)$ are in the structure, then $i < k < j < l$. This ensures that there is no crossing over of pairs which would result in pseudoknots.\n"
+ "4. No base appears in more than one pair.\n"
+ "\n"
+ "Using these rules we can begin to build our algorithm. The first part of our algorithm needs to take as input indices $i$ and $j$ and return the value $OPT(i,j)$ which is the optimal score of a structure between $i$ and $j$. We start by thinking about values of $i$ and $j$ for which we can immediately know the solution, this is known as a 'base case'. This is a case where the solution is known and no further recursion is required. Once the algorithm reaches the base case, it can return a solution and propagate it upward to the first recursive call. So once we have reached $i$ and $j$ that are too close to form a structure (rule number 2), we know that the score is 0. \n"
+ "\n"
+ "Otherwise, we must weigh the possibility of forming a pair or not forming a pair. If $i$ and $j$ are unpaired, then $OPT(i,j)$ is just $OPT(i, j-1)$ since the score will not increase for unpaired indices. \n"
+ "\n"
+ "The other case is that $i$ is paired to some index $t$ on the interval $[i,j]$. We then add 1 to the score and consider the structure formed before and after the pairing between $i$ and $t$. We can write these two cases as $OPT(i, t-1)$ and $OPT(t+1, j)$. But how do we know which $t$ to pair $i$ with? Well we simply try all possible values of $t$ within the allowed range and choose the best one. \n"
+ "\n"
+ "All of this can be summed up as follows:\n"
+ "\n"
+ "$$ OPT(i,j) = max\\begin{cases}\n"
+ " OPT(i, j-1) \\quad \\text{If $i$ and $j$ are not paired with each other.}\\\\\n"
+ " max(1 + OPT(i, t-1) + OPT(t+1, j)) \\quad \\text{Where we try all values of $t$ < j - 4}\n"
+ " \\end{cases}$$\n"
+ "\n"
+ "\n"
+ "We can now use this recursion to fill our dynamic programming table. Once we have filled the table with scores, we can retrieve the optimal folding by a process called backtracking. We won't go into detail on how this works, but the main idea is that we can start by looking at the entry containing the score for the full sequence $OPT[0][N]$. We can then look at adjacent entries and deduce which case (pairing or not pairing) resulted in the current value. We can continue like this for the full table until we have retrieved the full structure."
+ "where $N$ and $J$ are the number of discrete time and space points in our grid respectively.\n"
+ "$\\Delta t$ and $\\Delta x$ are the time step and space step respectively and defined as follows:\n"
+ "\n"
+ "$$\\Delta t = T / N,$$\n"
+ "\n"
+ "$$\\Delta x = L / J,$$\n"
+ "\n"
+ "where $T$ is the point in time up to which we will integrate $u$ numerically.\n"
+ "\n"
+ "Our ultimate goal is to construct a numerical method that allows us to approximate the unknonwn analytic solution $u(x,t)$\n"
+ "reasonably well in these discrete grid points.\n"
+ "\n"
+ "That is we want construct a method that computes values $U(j \\Delta x, n \\Delta t)$ (note: capital $U$) so that\n"
+ "\n"
+ "$$U(j \\Delta x, n \\Delta t) \\approx u(j \\Delta x, n \\Delta t)$$\n"
+ "\n"
+ "As a shorthand we will write $U_j^n = U(j \\Delta x, n \\Delta t)$ and $(j,n)$ to refer to grid point $(j \\Delta x, n \\Delta t)$."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## The Crank-Nicolson Stencil"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "Based on the two-dimensional grid we construct we then approximate the operators of our reaction-diffusion system.\n"
+ "\n"
+ "For instance, to approximate the time derivative on the left-hand side in grid point $(j,n)$ we use the values of $U$ in two specific grid points:\n"
+ "\n"
+ "$$\\frac{\\partial u}{\\partial t}\\Bigg|_{x = j \\Delta x, t = n \\Delta t} \\approx \\frac{U_j^{n+1} - U_j^n}{\\Delta t}.$$\n"
+ "\n"
+ "We can think of this scheme as a stencil that we superimpose on our $(x,t)$-grid and this particular stencil is\n"
+ "commonly referred to as [forward difference](http://en.wikipedia.org/wiki/Finite_difference#Forward.2C_backward.2C_and_central_differences).\n"
+ "\n"
+ "The spatial part of the [Crank-Nicolson stencil](http://journals.cambridge.org/abstract_S0305004100023197)\n"
+ "(or see [Table 3.2.2 of Trefethen](http://people.maths.ox.ac.uk/trefethen/3all.pdf))\n"
+ "for the heat equation ($u_t = u_{xx}$) approximates the \n"
+ "[Laplace operator](http://en.wikipedia.org/wiki/Laplace_operator) of our equation and takes the following form\n"
+ "To approximate $f(u(j \\Delta x, n \\Delta t))$ we write simply $f(U_j^n)$.\n"
+ "\n"
+ "These approximations define the stencil for our numerical method as pictured on [Wikipedia](http://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method).\n"
+ "Note that since our numerical integration starts with a well-defined initial condition at $n=0$, $\\mathbf{U}^0$, the\n"
+ "vector $\\mathbf{U}^{n+1}$ on the left-hand side is the only unknown in this system of linear equations.\n"
+ "\n"
+ "Thus, to integrate numerically our reaction-diffusion system from time point $n$ to $n+1$ we need to solve numerically for vector $\\mathbf{U}^{n+1}$.\n"
+ "\n"
+ "Let us call the matrix on the left-hand side $A$, the one on the right-hand side $B$,\n"
+ "and the vector on the right-hand side $\\mathbf{f}^n$.\n"
+ "Using this notation we can write the above system as\n"
+ "\n"
+ "$$A \\mathbf{U}^{n+1} = B \\mathbf{U}^n + f^n.$$\n"
+ "\n"
+ "In this linear equation, matrices $A$ and $B$ are defined by our problem: we need to specify these matrices once for our\n"
+ "problem and incorporate our boundary conditions in them.\n"
+ "Vector $\\mathbf{f}^n$ is a function of $\\mathbf{U}^n$ and so needs to be reevaluated in every time point $n$.\n"
+ "We also need to carry out one matrix-vector multiplication every time point, $B \\mathbf{U}^n$, and\n"
+ "Our one-dimensional domain has unit length and we define `J = 100` equally spaced\n"
+ "grid points in this domain.\n"
+ "This divides our domain into `J-1` subintervals, each of length `dx`."
+ ));
+
+ qDebug() << "command entry 3";
+ testCommandEntry(entry, 3, QString::fromUtf8(
+ "L = 1.\n"
+ "J = 100\n"
+ "dx = float(L)/float(J-1)\n"
+ "x_grid = numpy.array([j*dx for j in range(J)])"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "Equally, we define `N = 1000` equally spaced grid points on our time domain of length `T = 200` thus dividing our time domain into `N-1` intervals of length `dt`."
+ ));
+
+ qDebug() << "command entry 4";
+ testCommandEntry(entry, 4, QString::fromUtf8(
+ "T = 200\n"
+ "N = 1000\n"
+ "dt = float(T)/float(N-1)\n"
+ "t_grid = numpy.array([n*dt for n in range(N)])"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### Specify System Parameters and the Reaction Term"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "We choose our parameter values based on the work by\n"
+ "To advance our system by one time step, we need to do one matrix-vector multiplication followed by one vector-vector addition on the right hand side.\n"
+ "\n"
+ "To facilitate this, we rewrite our reaction term so that it accepts concentration vectors $\\mathbf{U}^n$ and $\\mathbf{V}^n$ as arguments\n"
+ "and returns vector $\\mathbf{f}^n$.\n"
+ "\n"
+ "As a reminder, this is our non-vectorial definition of $f$\n"
+ "\n"
+ " f = lambda u, v: v*(k0 + float(u*u)/float(1. + u*u)) - u"
+ ));
+
+ qDebug() << "command entry 12";
+ testCommandEntry(entry, 12, QString::fromUtf8(
+ "f_vec = lambda U, V: numpy.multiply(dt, numpy.subtract(numpy.multiply(V, \n"
+ "# Jove helps teach models of computation using Jupyter \n"
+ "\n"
+ "Included are modules on:\n"
+ "\n"
+ "* Sets, strings and languages\n"
+ "* Language operations\n"
+ "* Construction of and operations on DFA and NFA\n"
+ "* Regular expression parsing and automata inter-conversion\n"
+ "* Derivate-based parsing\n"
+ "* Pushdown automata\n"
+ "* The construction of parsers using context-free productions, including\n"
+ " a full lexer/parser for Jove's own markdown syntax\n"
+ "* Studies of parsing: ambiguity, associativity, precedence\n"
+ "* Turing machines (including one for the Collatz problem)\n"
+ "\n"
+ "For a complete Jove top-level reference, kindly refer to https://github.com/ganeshutah/Jove from where you can download and obtain Jove. You can also visit this Github link now and poke around (the NBViewer will display the contents).\n"
+ "\n"
+ "Once you are in the top-level Gallery link we provide, feel free to explore the hierarchy of modules found there.\n"
+ "\n"
+ "These notebooks should give you an idea of the contents.\n"
+ "\n"
+ "* [DFA Illustrations (has a Youtube)](http://nbviewer.jupyter.org/github/ganeshutah/Jove/blob/master/notebooks/tutorial/DFAUnit2.ipynb)\n"
+ "sys.path[0:0] = ['/home/mmmm1998/Документы/Репозитории/Jove','/home/mmmm1998/Документы/Репозитории/Jove/3rdparty'] # Put these at the head of the search path\n"
+ "from jove.DotBashers import *\n"
+ "from jove.Def_DFA import *\n"
+ "from jove.Def_NFA import *\n"
+ "from jove.Def_RE2NFA import *\n"
+ "from jove.Def_NFA2RE import *\n"
+ "from jove.Def_md2mc import *"
+ ));
+ testTextResult(entry, 0, QString::fromUtf8(
+ "You may use any of these help commands:\n"
+ "help(ResetStNum)\n"
+ "help(NxtStateStr)\n"
+ "\n"
+ "You may use any of these help commands:\n"
+ "help(mkp_dfa)\n"
+ "help(mk_dfa)\n"
+ "help(totalize_dfa)\n"
+ "help(addtosigma_delta)\n"
+ "help(step_dfa)\n"
+ "help(run_dfa)\n"
+ "help(accepts_dfa)\n"
+ "help(comp_dfa)\n"
+ "help(union_dfa)\n"
+ "help(intersect_dfa)\n"
+ "help(pruneUnreach)\n"
+ "help(iso_dfa)\n"
+ "help(langeq_dfa)\n"
+ "help(same_status)\n"
+ "help(h_langeq_dfa)\n"
+ "help(fixptDist)\n"
+ "help(min_dfa)\n"
+ "help(pairFR)\n"
+ "help(state_combos)\n"
+ "help(sepFinNonFin)\n"
+ "help(bash_eql_classes)\n"
+ "help(listminus)\n"
+ "help(bash_1)\n"
+ "help(mk_rep_eqc)\n"
+ "help(F_of)\n"
+ "help(rep_of_s)\n"
+ "help(q0_of)\n"
+ "help(Delta_of)\n"
+ "help(mk_state_eqc_name)\n"
+ "\n"
+ "You may use any of these help commands:\n"
+ "help(mk_nfa)\n"
+ "help(totalize_nfa)\n"
+ "help(step_nfa)\n"
+ "help(run_nfa)\n"
+ "help(ec_step_nfa)\n"
+ "help(Eclosure)\n"
+ "help(Echelp)\n"
+ "help(accepts_nfa)\n"
+ "help(nfa2dfa)\n"
+ "help(n2d)\n"
+ "help(inSets)\n"
+ "help(rev_dfa)\n"
+ "help(min_dfa_brz)\n"
+ "\n"
+ "You may use any of these help commands:\n"
+ "help(re2nfa)\n"
+ "\n"
+ "You may use any of these help commands:\n"
+ "help(RE2Str)\n"
+ "help(mk_gnfa)\n"
+ "help(mk_gnfa_from_D)\n"
+ "help(dfa2nfa)\n"
+ "help(del_gnfa_states)\n"
+ "help(gnfa_w_REStr)\n"
+ "help(del_one_gnfa_state)\n"
+ "help(Edges_Exist_Via)\n"
+ "help(choose_state_to_del)\n"
+ "help(form_alt_RE)\n"
+ "help(form_concat_RE)\n"
+ "help(form_kleene_RE)\n"
+ "\n"
+ "You may use any of these help commands:\n"
+ "help(md2mc)\n"
+ ".. and if you want to dig more, then ..\n"
+ "help(default_line_attr)\n"
+ "help(length_ok_input_items)\n"
+ "help(union_line_attr_list_fld)\n"
+ "help(extend_rsltdict)\n"
+ "help(form_delta)\n"
+ "help(get_machine_components)"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ " # Jove allows you to set problems in markdown and have students solve"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "1) LOdd1Three0 : Set of strings over {0,1} with an odd # of 1s OR exactly three 0s. \n"
+ "\n"
+ "* Hint on how to arrive at the language:\n"
+ "\n"
+ " - develop NFAs for the two cases and perform their union. Obtain DFA\n"
+ "\n"
+ " - develop REs for the two cases and perform the union. \n"
+ "\n"
+ " - Testing the creations:\n"
+ "\n"
+ " . Come up with language for even # of 1s and separately for \"other than three 0s\". \n"
+ " \n"
+ " . Do two intersections. \n"
+ " \n"
+ " . Is the language empty?\n"
+ "\n"
+ "\n"
+ "2) Language of strings over {0,1} with exactly two occurrences of 0101 in it.\n"
+ "\n"
+ " * Come up with it directly (take overlaps into account, i.e. 010101 has two occurrences in it\n"
+ "\n"
+ " * Come up in another way\n"
+ "\n"
+ "Notes:\n"
+ "\n"
+ "* Most of the problem students will have in this course is interpreting English (technical English)\n"
+ "\n"
+ "* So again, read the writeup at the beginning of Module6 (should be ready soon today) and work on using the tool.\n"
+ "\n"
+ "\n"
+ "\n"
+ ""
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "__Solutions__\n"
+ "\n"
+ "1) LOdd1Three0 : Set of strings over {0,1} with an odd # of 1s OR exactly three 0s. \n"
+ "\n"
+ "* Hint on how to arrive at the language:\n"
+ "\n"
+ " - develop NFAs for the two cases and perform their union. Obtain DFA\n"
+ "\n"
+ " - develop REs for the two cases and perform the union. \n"
+ "\n"
+ " - Testing the creations:\n"
+ "\n"
+ " . Come up with language for even # of 1s and separately for \"other than three 0s\". \n"
+ " \n"
+ " . Do two intersections. \n"
+ " \n"
+ " . Is the language empty?\n"
+ "\n"
+ "\n"
+ "2) Language of strings over {0,1} with exactly two occurrences of 0101 in it.\n"
+ "\n"
+ " * Come up with it directly (take overlaps into account, i.e. 010101 has two occurrences in it\n"
+ "\n"
+ " * Come up in another way\n"
+ "\n"
+ "Notes:\n"
+ "\n"
+ "* Most of the problem students will have in this course is interpreting English (technical English)\n"
+ "\n"
+ "* So again, read the writeup at the beginning of Module6 (should be ready soon today) and work on using the tool.\n"
+ "<p>We live in a complex environment and must constantly integrate sensory information to interact with the world around us. Inputs from different modalities might not always be congruent with each other, but dissociating the true nature of the stimulus may be a matter of life or death for an organism.</p>\n"
+ "<p>You hear and see evidence of a rattlesnake in tall grass near you. You get an auditory and a visual cue of the snake's location $s$. Both cues are associated with a likelihood function indicating the probability of that cue for all possible locations of the snake. The likelihood function associated with the visual cue, $p(c_\\mathrm{V}|s)$, has high uncertainty, because of the tall grass. The auditory cue is easier to localize, so its associated likelihood function, $p(c_\\mathrm{A}|s)$, is sharper. In accordance with Bayes' Rule, and assuming a flat prior over the snake's location, an optimal estimate of the location of the snake can be computed by multiplying the two likelihoods. This joint likelihood will be between the two cues but closer to the less uncertain cue, and will have less uncertainty than both unimodal likelihood functions.</p>"
+ "Behavioral experiments have demonstrated that humans perform near-optimal Bayesian inference on ambiguous sensory information (van Beers *et al.*, 1999; Ernst & Banks, 2002; Kording & Wolpert, 2004; Stocker & Simoncelli, 2006). This has been demonstrated in cue combination experiments in which subjects report a near-optimal estimate of the stimulus given two noisy measurements of that stimulus. However, the neural basis for how humans might perform these computations is unclear. \n"
+ "\n"
+ "Ma *et. al.* (2006) propose that variance in cortical activity, rather than impairing sensory systems, is an adaptive mechanism to encode uncertainty in sensory measurements. They provide theory showing how the brain might use probabilistic population codes to perform near-optimal cue combination. We will re-derive the theory in here, and demonstrate it by simulating and decoding neural populations.\n"
+ "\n"
+ "## Cues can be represented by neural populations\n"
+ "\n"
+ "To return to our deadly rattlesnake, let's now assume that $c_\\mathrm{V}$ and $c_\\mathrm{A}$ are represented by populations of neurons $\\mathbf{r}_\\mathrm{V}$ and $\\mathbf{r}_\\mathrm{A}$, respectively. For our math and simulations, we assume that $\\mathbf{r}_\\mathrm{V}$ and $\\mathbf{r}_\\mathrm{A}$ are each composed of $N$ neurons that:\n"
+ "\n"
+ "* have independent Poisson variability\n"
+ "* have regularly spaced Gaussian tuning curves that are identical in mean and variance for neurons with the same index in both populations\n"
+ "\n"
+ "The populations may have different gains, $g_\\mathrm{V}$ and $g_\\mathrm{A}$.\n"
+ "\n"
+ "These are the tuning curves for the neurons in $\\mathbf{r}_\\mathrm{V}$ (purple) and $\\mathbf{r}_\\mathrm{A}$ (red). Each curve represents the mean firing rate of a single neuron given a location $s$. Each neuron thus has a preferred location, which is where its tuning curve peaks."
+ "The tuning curves are dense enough that we can also assume that $\\sum_{i=0}^N f_i(s) = k$ (*i.e.*, the sum of the tuning curves in a population is constant.)"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "\n"
+ "\n"
+ ""
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "\n"
+ "\n"
+ ""
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "First, we will show how the brain can decode a likelihood over stimulus from neural activity. Then we will ask how the brain can compute joint likelihoods.\n"
+ "### How can the brain decode $p(\\mathbf{r_\\mathrm{V}}|s)$?\n"
+ "## How can the brain compute the joint likelihood $p(\\mathbf{r}_\\mathrm{V}|s)\\ p(\\mathbf{r}_\\mathrm{A}|s)$?\n"
+ "The fact that we see neurons from $\\mathbf{r}_\\mathrm{V}$ and $\\mathbf{r}_\\mathrm{A}$ being added on a neuron-by-neuron basis in the exponent above suggests that we could construct a third population vector, $\\mathbf{r}_\\mathrm{V}+\\mathbf{r}_\\mathrm{A}$, and decode that.\n"
+ "\n"
+ "### First, we must prove that the sum of two Poisson-distributed random variables $X+Y$ is again Poisson-distributed.\n"
+ "Since equations $(9)$ and $(27)$ are proportional, we have shown that optimal cue combination can be executed by decoding linear sums of populations."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "$$x = 2$$"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "\n"
+ "\n"
+ ""
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "\n"
+ "\n"
+ ""
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Simulation\n"
+ "Here are the spike counts (during 1 s) from the two populations on one trial. Depicted in blue is a third population vector that is the sum of $\\mathbf{r}_\\mathrm{V}$ and $\\mathbf{r}_\\mathrm{A}$."
+ "Here are the decoded likelihoods for each population alone $(6)$, the joint likelihood $(9)$, and the likelihood for the summed population $(27)$. Note that the joint likelihood (gold) is less uncertain than either unimodal likelihood. Also note that it is identical to the likelihood for the summed population (blue)."
+ "Here, we break the assumption that the two populations have the same tuning curve width. Note that the joint likelihood (gold) is no longer identical to the likelihood for the summed population (blue)."
+ "Now you can play interactively with the parameters of the simulation using these sliders, and watch the decoded likelihoods shift around. Every time you change a parameter, new sets of spikes are generated and used to infer $s$.\n"
+ "\n"
+ "For the simulation to be interactive, you'll have to download this notebook."
+ "i = ipywidgets.interactive(spikes_and_inference,\n"
+ " true_stimulus = (-40, 40, .1),\n"
+ " number_of_neurons = (2, 200, 1),\n"
+ " r_V_gain = (0, 100, 1),\n"
+ " r_A_gain = (0, 100, 1),\n"
+ " r_V_tuning_curve_sigma = (0.1, 50, .1),\n"
+ " r_A_tuning_curve_sigma = (0.1, 50, .1),\n"
+ " tuning_curve_baseline = (0, 20, .1));\n"
+ "display(ipywidgets.VBox(i.children[2:-1]))"
+ ));
+ testImageResult(entry, 0);
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "\n"
+ "\n"
+ ""
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "\n"
+ "\n"
+ ""
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Conclusion\n"
+ "\n"
+ "It has been shown behaviorally that humans perform near-optimal Bayesian inference on ambiguous sensory information. As suggested by Ma *et. al.* (2006) and shown here, it is possible that the brain does this operation by simply performing linear combinations of populations of Poisson neurons receiving various sensory input. Cortical neurons may be particularly well suited for this task because they have Poisson-like firing rates, displaying reliable variability from trial to trial (Tolhurst, Movshon & Dean, 1982; Softky & Koch, 1993).\n"
+ "\n"
+ "High levels of noise in these populations might at first be difficult to reconcile considering highly precise behavioral data. However, variability in neural populations might be direcly representative of uncertainty in environmental stimuli. Variability in cortical populations would then be critical for precise neural coding.\n"
+ "\n"
+ "## References\n"
+ "\n"
+ "* Ernst MO, Banks MS. (2002). Humans integrate visual and haptic information in a statistically optimal fashion. *Nature.*\n"
+ "* Ma WJ, Beck JM, Latham PE, Pouget A. (2006). Bayesian inference with probabilistic population codes. *Nature Neuroscience.*\n"
+ "* Softky WR, Koch C. (1993). The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. *Journal of Neuroscience.*\n"
+ "* Stocker AA, Simoncelli EP. (2006). Noise characteristics and prior expectations in human visual speed perception. *Nature Neuroscience.*\n"
+ "* Tolhurst, DJ, Movshon JA, Dean AF. (1983). The statistical reliability of signals in single neurons in cat and monkey visual cortex. *Vision Research.*\n"
+ "* van Beers RJ, Sittig AC, Gon JJ. (1999). Integration of proprioceptive and visual position-information: An experimentally supported model. *Journal of Neurophysiology.*"
+ "# Rigid-body transformations in a plane (2D)\n"
+ "\n"
+ "> Marcos Duarte \n"
+ "> Laboratory of Biomechanics and Motor Control ([http://demotu.org/](http://demotu.org/)) \n"
+ "> Federal University of ABC, Brazil"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "The kinematics of a rigid body is completely described by its pose, i.e., its position and orientation in space (and the corresponding changes are translation and rotation). The translation and rotation of a rigid body are also known as rigid-body transformations (or simply, rigid transformations).\n"
+ "\n"
+ "Remember that in physics, a [rigid body](https://en.wikipedia.org/wiki/Rigid_body) is a model (an idealization) for a body in which deformation is neglected, i.e., the distance between every pair of points in the body is considered constant. Consequently, the position and orientation of a rigid body can be completely described by a corresponding coordinate system attached to it. For instance, two (or more) coordinate systems can be used to represent the same rigid body at two (or more) instants or two (or more) rigid bodies in space.\n"
+ "\n"
+ "Rigid-body transformations are used in motion analysis (e.g., of the human body) to describe the position and orientation of each segment (using a local (anatomical) coordinate system defined for each segment) in relation to a global coordinate system fixed at the laboratory. Furthermore, one can define an additional coordinate system called technical coordinate system also fixed at the rigid body but not based on anatomical landmarks. In this case, the position of the technical markers is first described in the laboratory coordinate system, and then the technical coordinate system is calculated to recreate the anatomical landmarks position in order to finally calculate the original anatomical coordinate system (and obtain its unknown position and orientation through time).\n"
+ "\n"
+ "In what follows, we will study rigid-body transformations by looking at the transformations between two coordinate systems. For simplicity, let's first analyze planar (two-dimensional) rigid-body transformations and later we will extend these concepts to three dimensions (where the study of rotations are more complicated)."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Affine transformations\n"
+ "\n"
+ "Translation and rotation are two examples of [affine transformations](https://en.wikipedia.org/wiki/Affine_transformation). Affine transformations preserve straight lines, but not necessarily the distance between points. Other examples of affine transformations are scaling, shear, and reflection. The figure below illustrates different affine transformations in a plane. Note that a 3x3 matrix is shown on top of each transformation; these matrices are known as the transformation matrices and are the mathematical representation of the physical transformations. Next, we will study how to use this approach to describe the translation and rotation of a rigid-body. \n"
+ "<br>\n"
+ "<figure><img src='https://upload.wikimedia.org/wikipedia/commons/thumb/2/2c/2D_affine_transformation_matrix.svg/360px-2D_affine_transformation_matrix.svg.png' alt='Affine transformations'/> <figcaption><center><i>Figure. Examples of affine transformations in a plane applied to a square (with the letter <b>F</b> in it) and the corresponding transformation matrices (<a href=\"https://en.wikipedia.org/wiki/Affine_transformation\">image from Wikipedia</a>).</i></center></figcaption> </figure>"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Translation\n"
+ "\n"
+ "In a two-dimensional space, two coordinates and one angle are sufficient to describe the pose of the rigid body, totalizing three degrees of freedom for a rigid body. Let's see first the transformation for translation, then for rotation, and combine them at last.\n"
+ "\n"
+ "A pure two-dimensional translation of a coordinate system in relation to other coordinate system and the representation of a point in these two coordinate systems are illustrated in the figure below (remember that this is equivalent to describing a translation between two rigid bodies). \n"
+ "<br>\n"
+ "<figure><img src='./../images/translation2D.png' alt='translation 2D'/> <figcaption><center><i>Figure. A point in two-dimensional space represented in two coordinate systems (Global and local), with one system translated.</i></center></figcaption> </figure>\n"
+ "\n"
+ "The position of point $\\mathbf{P}$ originally described in the local coordinate system but now described in the Global coordinate system in vector form is:\n"
+ "Because position and translation can be treated as vectors, the inverse operation, to describe the position at the local coordinate system in terms of the Global coordinate system, is simply:\n"
+ "From classical mechanics, this transformation is an example of [Galilean transformation](http://en.wikipedia.org/wiki/Galilean_transformation). \n"
+ "\n"
+ "For example, if the local coordinate system is translated by $\\mathbf{L_G}=[2, 3]$ in relation to the Global coordinate system, a point with coordinates $\\mathbf{P_l}=[4, 5]$ at the local coordinate system will have the position $\\mathbf{P_G}=[6, 8]$ at the Global coordinate system:"
+ "Pl = np.array([[4, 5], [6, 7], [8, 9]]) # 2D array with 3 rows and two columns\n"
+ "PG = LG + Pl\n"
+ "PG"
+ ));
+ testTextResult(entry, 0, QString::fromLatin1(
+ "array([[ 6, 8],\n"
+ " [ 8, 10],\n"
+ " [10, 12]])"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Rotation\n"
+ "\n"
+ "A pure two-dimensional rotation of a coordinate system in relation to other coordinate system and the representation of a point in these two coordinate systems are illustrated in the figure below (remember that this is equivalent to describing a rotation between two rigid bodies). The rotation is around an axis orthogonal to this page, not shown in the figure (for a three-dimensional coordinate system the rotation would be around the $\\mathbf{Z}$ axis). \n"
+ "<br>\n"
+ "<figure><img src='./../images/rotation2D.png' alt='rotation 2D'/> <figcaption><center><i>Figure. A point in the two-dimensional space represented in two coordinate systems (Global and local), with one system rotated in relation to the other around an axis orthogonal to both coordinate systems.</i></center></figcaption> </figure>\n"
+ "\n"
+ "Consider we want to express the position of point $\\mathbf{P}$ in the Global coordinate system in terms of the local coordinate system knowing only the coordinates at the local coordinate system and the angle of rotation between the two coordinate systems. \n"
+ "\n"
+ "There are different ways of deducing that, we will see three of these methods next. "
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### 1. Using trigonometry\n"
+ "\n"
+ "From figure below, the coordinates of point $\\mathbf{P}$ in the Global coordinate system can be determined finding the sides of the triangles marked in red. \n"
+ "<br>\n"
+ "<figure><img src='./../images/rotation2Db.png' alt='rotation 2D'/> <figcaption><center><i>Figure. The coordinates of a point at the Global coordinate system in terms of the coordinates of this point at the local coordinate system.</i></center></figcaption> </figure>\n"
+ "Where $\\mathbf{R_{Gl}}$ is the rotation matrix that rotates the coordinates from the local to the Global coordinate system:\n"
+ "\n"
+ "$$ \\mathbf{R_{Gl}} = \\begin{bmatrix}\n"
+ "\\cos\\alpha & -\\sin\\alpha \\\\\n"
+ "\\sin\\alpha & \\cos\\alpha \n"
+ "\\end{bmatrix} $$\n"
+ "\n"
+ "So, given any position at the local coordinate system, with the rotation matrix above we are able to determine the position at the Global coordinate system. Let's check that before looking at other methods to obtain this matrix. \n"
+ "\n"
+ "For instance, consider a local coordinate system rotated by $45^o$ in relation to the Global coordinate system, a point in the local coordinate system with position $\\mathbf{P_l}=[1, 1]$ will have the following position at the Global coordinate system:"
+ "PG = np.dot(RGl, Pl) # the function dot() is used for matrix multiplication with arrays\n"
+ "np.around(PG, 4) # round the number due to floating point arithmetic errors"
+ ));
+ testTextResult(entry, 0, QString::fromLatin1(
+ "array([[ 0. , -0.7071, 0.7071],\n"
+ " [ 1.4142, 0.7071, 0.7071]])"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "We have done all the calculations using the array function in NumPy. A [NumPy array is different than a matrix](http://www.scipy.org/NumPy_for_Matlab_Users), if we want to use explicit matrices in NumPy, the calculation above will be:"
+ "PG = RGl*Pl # matrix multiplication in NumPy\n"
+ "np.around(PG, 4) # round the number due to floating point arithmetic errors"
+ ));
+ testTextResult(entry, 0, QString::fromLatin1(
+ "array([[ 0. , -0.7071, 0.7071],\n"
+ " [ 1.4142, 0.7071, 0.7071]])"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "Both array and matrix types work in NumPy, but you should choose only one type and not mix them; the array is preferred because it is [the standard vector/matrix/tensor type of NumPy](http://www.scipy.org/NumPy_for_Matlab_Users)."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### 2. Using direction cosines\n"
+ "\n"
+ "Another way to determine the rotation matrix is to use the concept of direction cosine. \n"
+ "\n"
+ "> Direction cosines are the cosines of the angles between any two vectors. \n"
+ "\n"
+ "For the present case with two coordinate systems, they are the cosines of the angles between each axis of one coordinate system and each axis of the other coordinate system. The figure below illustrates the directions angles between the two coordinate systems, expressing the local coordinate system in terms of the Global coordinate system. \n"
+ "<br>\n"
+ "<figure><img src='./../images/directioncosine2D.png' alt='direction angles 2D'/> <figcaption><center><i>Figure. Definition of direction angles at the two-dimensional space.</i></center></figcaption> </figure> \n"
+ "<br>\n"
+ "$$ \\mathbf{R_{Gl}} = \\begin{bmatrix}\n"
+ "\\cos\\mathbf{X}x & \\cos\\mathbf{X}y \\\\\n"
+ "\\cos\\mathbf{Y}x & \\cos\\mathbf{Y}y \n"
+ "\\end{bmatrix} = \n"
+ "\\begin{bmatrix}\n"
+ "\\cos(\\alpha) & \\cos(90^o+\\alpha) \\\\\n"
+ "\\cos(90^o-\\alpha) & \\cos(\\alpha)\n"
+ "\\end{bmatrix} = \n"
+ "\\begin{bmatrix}\n"
+ "\\cos\\alpha & -\\sin\\alpha \\\\\n"
+ "\\sin\\alpha & \\cos\\alpha \n"
+ "\\end{bmatrix} $$ \n"
+ "\n"
+ "The same rotation matrix as obtained before.\n"
+ "\n"
+ "Note that the order of the direction cosines is because in our convention, the first row is for the $\\mathbf{X}$ coordinate and the second row for the $\\mathbf{Y}$ coordinate (the outputs). For the inputs, we followed the same order, first column for the $\\mathbf{x}$ coordinate, second column for the $\\mathbf{y}$ coordinate."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### 3. Using a basis\n"
+ "\n"
+ "Another way to deduce the rotation matrix is to view the axes of the rotated coordinate system as unit vectors, versors, of a <a href=\"http://en.wikipedia.org/wiki/Basis_(linear_algebra)\">basis</a> as illustrated in the figure below.\n"
+ "\n"
+ "> A basis is a set of linearly independent vectors that can represent every vector in a given vector space, i.e., a basis defines a coordinate system.\n"
+ "\n"
+ "<figure><img src='./../images/basis2D2.png' alt='basis 2D'/> <figcaption><center><i>Figure. Definition of the rotation matrix using a basis at the two-dimensional space.</i></center></figcaption> </figure>\n"
+ "\n"
+ "The coordinates of these two versors at the local coordinate system in terms of the Global coordinate system are:\n"
+ "Note that as unit vectors, each of the versors above should have norm (length) equals to one, which indeed is the case.\n"
+ "\n"
+ "If we express each versor above as different columns of a matrix, we obtain the rotation matrix again: \n"
+ "\n"
+ "$$ \\mathbf{R_{Gl}} = \\begin{bmatrix}\n"
+ "\\cos\\alpha & -\\sin\\alpha \\\\\\\n"
+ "\\sin\\alpha & \\cos\\alpha \n"
+ "\\end{bmatrix} $$\n"
+ "\n"
+ "This means that the rotation matrix can be viewed as the basis of the rotated coordinate system defined by its versors. \n"
+ "\n"
+ "This third way to derive the rotation matrix is in fact the method most commonly used in motion analysis because the coordinates of markers (in the Global/laboratory coordinate system) are what we measure with cameras. "
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### 4. Using the inner (dot or scalar) product between versors\n"
+ "\n"
+ "Yet another way to deduce the rotation matrix is to define it as the dot product between the versors of the bases related to the two coordinate systems:\n"
+ "And the rotation matrix will be equal to the matrix deduced based on the direction cosines."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### Local-to-Global and Global-to-local coordinate systems' rotations"
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "If we want the inverse operation, to express the position of point $\\mathbf{P}$ in the local coordinate system in terms of the Global coordinate system, the figure below illustrates that using trigonometry. \n"
+ "<br>\n"
+ "<figure><img src='./../images/rotation2Dc.png' alt='rotation 2D'/> <figcaption><center><i>Figure. The coordinates of a point at the local coordinate system in terms of the coordinates at the Global coordinate system.</i></center></figcaption> </figure>\n"
+ "Where $\\mathbf{R_{lG}}$ is the rotation matrix that rotates the coordinates from the Global to the local coordinate system (note the inverse order of the subscripts):\n"
+ "\n"
+ "$$ \\mathbf{R_{lG}} = \\begin{bmatrix}\n"
+ "\\cos\\alpha & \\sin\\alpha \\\\\n"
+ "-\\sin\\alpha & \\cos\\alpha \n"
+ "\\end{bmatrix} $$\n"
+ "\n"
+ "If we use the direction cosines to calculate the rotation matrix, because the axes didn't change, the cosines are the same, only the order changes, now $\\mathbf{x, y}$ are the rows (outputs) and $\\mathbf{X, Y}$ are the columns (inputs):\n"
+ "\n"
+ "$$ \\mathbf{R_{lG}} = \\begin{bmatrix}\n"
+ "\\cos\\mathbf{X}x & \\cos\\mathbf{Y}x \\\\\n"
+ "\\cos\\mathbf{X}y & \\cos\\mathbf{Y}y \n"
+ "\\end{bmatrix} = \n"
+ "\\begin{bmatrix}\n"
+ "\\cos(\\alpha) & \\cos(90^o-\\alpha) \\\\\n"
+ "\\cos(90^o+\\alpha) & \\cos(\\alpha)\n"
+ "\\end{bmatrix} = \n"
+ "\\begin{bmatrix}\n"
+ "\\cos\\alpha & \\sin\\alpha \\\\\n"
+ "-\\sin\\alpha & \\cos\\alpha \n"
+ "\\end{bmatrix} $$\n"
+ "\n"
+ "And defining the versors of the axes in the Global coordinate system for a basis in terms of the local coordinate system would also produce this latter rotation matrix.\n"
+ "\n"
+ "The two sets of equations and matrices for the rotations from Global-to-local and local-to-Global coordinate systems are very similar, this is no coincidence. Each of the rotation matrices we deduced, $\\mathbf{R_{Gl}}$ and $\\mathbf{R_{lG}}$, perform the inverse operation in relation to the other. Each matrix is the inverse of the other. \n"
+ "\n"
+ "In other words, the relation between the two rotation matrices means it is equivalent to instead of rotating the local coordinate system by $\\alpha$ in relation to the Global coordinate system, to rotate the Global coordinate system by $-\\alpha$ in relation to the local coordinate system; remember that $\\cos(-\\alpha)=\\cos(\\alpha)$ and $\\sin(-\\alpha)=-\\sin(\\alpha)$."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### Rotation of a Vector\n"
+ "\n"
+ "We can also use the rotation matrix to rotate a vector by a given angle around an axis of the coordinate system as shown in the figure below. \n"
+ "<br>\n"
+ "<figure><img src='./../images/rotation2Dvector.png' alt='rotation 2D of a vector'/> <figcaption><center><i>Figure. Rotation of a position vector $\\mathbf{P}$ by an angle $\\alpha$ in the two-dimensional space.</i></center></figcaption> </figure>\n"
+ "\n"
+ "We will not prove that we use the same rotation matrix, but think that in this case the vector position rotates by the same angle instead of the coordinate system. The new coordinates of the vector position $\\mathbf{P'}$ rotated by an angle $\\alpha$ is simply the rotation matrix (for the angle $\\alpha$) multiplied by the coordinates of the vector position $\\mathbf{P}$:\n"
+ "**[See here for a review about matrix and its main properties](http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/Matrix.ipynb)**.\n"
+ "\n"
+ "A nice property of the rotation matrix is that its inverse is the transpose of the matrix (because the columns/rows are mutually orthogonal and have norm equal to one). \n"
+ "This property can be shown with the rotation matrices we deduced:\n"
+ "Using the inverse and the transpose mathematical operations, the coordinates at the local coordinate system given the coordinates at the Global coordinate system and the rotation matrix can be obtained by: \n"
+ "Where we referred the inverse of $\\mathbf{R_{Gl}}\\;(\\:\\mathbf{R_{Gl}^{-1}})$ as $\\mathbf{R_{lG}}$ (note the different order of the subscripts). \n"
+ "Pl = np.array([[1, 1]]).T # transpose the array for correct matrix multiplication\n"
+ "print('Position at the local coordinate system (Pl):\\n', Pl)\n"
+ "\n"
+ "PG = np.dot(RGl, Pl) # the function dot() is used for matrix multiplication with arrays\n"
+ "print('Position at the Global coordinate system (PG=RGl*Pl):\\n', np.around(PG,2))\n"
+ "\n"
+ "Pl = np.dot(RGl.T, PG)\n"
+ "print('Position at the local coordinate system using the inverse of RGl (Pl=RlG*PG):\\n', Pl)"
+ ));
+ testTextResult(entry, 0, QString::fromUtf8(
+ "Rotation matrix (RGl):\n"
+ " [[ 0.7071 -0.7071]\n"
+ " [ 0.7071 0.7071]]\n"
+ "Position at the local coordinate system (Pl):\n"
+ " [[1]\n"
+ " [1]]\n"
+ "Position at the Global coordinate system (PG=RGl*Pl):\n"
+ " [[0. ]\n"
+ " [1.41]]\n"
+ "Position at the local coordinate system using the inverse of RGl (Pl=RlG*PG):\n"
+ " [[1.]\n"
+ " [1.]]"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "**In summary, some of the properties of the rotation matrix are:** \n"
+ "1. The columns of the rotation matrix form a basis of (independent) unit vectors (versors) and the rows are also independent versors since the transpose of the rotation matrix is another rotation matrix. \n"
+ "2. The rotation matrix is orthogonal. There is no linear combination of one of the lines or columns of the matrix that would lead to the other row or column, i.e., the lines and columns of the rotation matrix are independent, orthogonal, to each other (this is property 1 rewritten). Because each row and column have norm equal to one, this matrix is also sometimes said to be orthonormal. \n"
+ "3. The determinant of the rotation matrix is equal to one (or equal to -1 if a left-hand coordinate system was used, but you should rarely use that). For instance, the determinant of the rotation matrix we deduced is $cos\\alpha cos\\alpha - sin\\alpha(-sin\\alpha)=1$.\n"
+ "4. The inverse of the rotation matrix is equals to its transpose.\n"
+ "\n"
+ "**On the different meanings of the rotation matrix:** \n"
+ "- It represents the coordinate transformation between the coordinates of a point expressed in two different coordinate systems. \n"
+ "- It describes the rotation between two coordinate systems. The columns are the direction cosines (versors) of the axes of the rotated coordinate system in relation to the other coordinate system and the rows are also direction cosines (versors) for the inverse rotation. \n"
+ "- It is an operator for the calculation of the rotation of a vector in a coordinate system.\n"
+ "- Rotation matrices provide a means of numerically representing rotations without appealing to angular specification.\n"
+ "\n"
+ "**Which matrix to use, from local to Global or Global to local?** \n"
+ "- A typical use of the transformation is in movement analysis, where there are the fixed Global (laboratory) coordinate system and the local (moving, e.g. anatomical) coordinate system attached to each body segment. Because the movement of the body segment is measured in the Global coordinate system, using cameras for example, and we want to reconstruct the coordinates of the markers at the anatomical coordinate system, we want the transformation leading from the Global coordinate system to the local coordinate system.\n"
+ "- Of course, if you have one matrix, it is simple to get the other; you just have to pay attention to use the right one."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Translation and rotation\n"
+ "\n"
+ "Consider now the case where the local coordinate system is translated and rotated in relation to the Global coordinate system and a point is described in both coordinate systems as illustrated in the figure below (once again, remember that this is equivalent to describing a translation and a rotation between two rigid bodies). \n"
+ "<br>\n"
+ "<figure><img src='./../images/transrot2D.png' alt='translation and rotation 2D'/> <figcaption><center><i>Figure. A point in two-dimensional space represented in two coordinate systems, with one system translated and rotated.</i></center></figcaption> </figure>\n"
+ "\n"
+ "The position of point $\\mathbf{P}$ originally described in the local coordinate system, but now described in the Global coordinate system in vector form is:\n"
+ "This means that we first *disrotate* the local coordinate system and then correct for the translation between the two coordinate systems. Note that we can't invert this order: the point position is expressed in the local coordinate system and we can't add this vector to another vector expressed in the Global coordinate system, first we have to convert the vectors to the same coordinate system.\n"
+ "\n"
+ "If now we want to find the position of a point at the local coordinate system given its position in the Global coordinate system, the rotation matrix and the translation vector, we have to invert the expression above:\n"
+ "The expression above indicates that to perform the inverse operation, to go from the Global to the local coordinate system, we first translate and then rotate the coordinate system."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### Transformation matrix\n"
+ "\n"
+ "It is possible to combine the translation and rotation operations in only one matrix, called the transformation matrix (also referred as homogeneous transformation matrix):\n"
+ "However, because $\\mathbf{T_{Gl}}$ is not orthonormal when there is a translation, its inverse is not its transpose. Its inverse in matrix form is given by:\n"
+ "A typical scenario in motion analysis is to calculate the rotation matrix using the position of markers placed on the moving rigid body. With the markers' positions, we create a local basis, which by definition is the rotation matrix for the rigid body with respect to the Global (laboratory) coordinate system. To define a coordinate system using a basís, we also will need to define an origin. \n"
+ "\n"
+ "Let's see how to calculate a basis given the markers' positions. \n"
+ "Consider the markers at m1=[1,1]`, m2=[1,2] and m3=[-1,1] measured in the Global coordinate system as illustrated in the figure below: \n"
+ "<br>\n"
+ "<figure><img src='./../images/transrot2Db.png' alt='translation and rotation 2D'/> <figcaption><center><i>Figure. Three points in the two-dimensional space, two possible vectors given these points, and the corresponding basis.</i></center></figcaption> </figure>\n"
+ "\n"
+ "A possible local coordinate system with origin at the position of m1 is also illustrated in the figure above. Intentionally, the three markers were chosen to form orthogonal vectors. \n"
+ "The translation vector between the two coordinate system is:\n"
+ "\n"
+ "$$\\mathbf{L_{Gl}} = m_1 - [0,0] = [1,1]$$\n"
+ "\n"
+ "The vectors expressing the axes of the local coordinate system are:\n"
+ "Note that these two vectors do not form a basis yet because they are not unit vectors (in fact, only *y* is not a unit vector). Let's normalize these vectors:\n"
+ "Beware that the versors above are not exactly the same as the ones shown in the right plot of the last figure, the versors above if plotted will start at the origin of the coordinate system, not at [1,1] as shown in the figure.\n"
+ "\n"
+ "We could have done this calculation in NumPy (we will need to do that when dealing with real data later):"
+ "print(\"x =\", x, \", y =\", y, \"\\nex=\", vx, \", ey=\", vy)"
+ ));
+ testTextResult(entry, 0, QString::fromUtf8(
+ "x = [0. 1.] , y = [-2. 0.] \n"
+ "ex= [0. 1.] , ey= [-1. 0.]"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "Now, both $\\mathbf{e}_x$ and $\\mathbf{e}_y$ are unit vectors (versors) and they are orthogonal, a basis can be formed with these two versors, and we can represent the rotation matrix using this basis (just place the versors of this basis as columns of the rotation matrix):\n"
+ "\n"
+ "$$ \\mathbf{R_{Gl}} = \\begin{bmatrix}\n"
+ "0 & -1 \\\\\n"
+ "1 & 0 \n"
+ "\\end{bmatrix} $$\n"
+ "\n"
+ "This rotation matrix makes sense because from the figure above we see that the local coordinate system we defined is rotated by 90$^o$ in relation to the Global coordinate system and if we use the general form for the rotation matrix:\n"
+ "\n"
+ "$$ \\mathbf{R} = \\begin{bmatrix}\n"
+ "\\cos\\alpha & -\\sin\\alpha \\\\\n"
+ "\\sin\\alpha & \\cos\\alpha \n"
+ "\\end{bmatrix} = \n"
+ "\\begin{bmatrix}\n"
+ "\\cos90^o & -\\sin90^o \\\\\n"
+ "\\sin90^o & \\cos90^o \n"
+ "\\end{bmatrix} =\n"
+ "\\begin{bmatrix}\n"
+ "0 & -1 \\\\\n"
+ "1 & 0 \n"
+ "\\end{bmatrix} $$\n"
+ "\n"
+ "So, the position of any point in the local coordinate system can be represented in the Global coordinate system by:\n"
+ "print('Position at the local coordinate system:\\n', Pl)\n"
+ "\n"
+ "PG = LGl + np.dot(RGl, Pl)\n"
+ "print('Position at the Global coordinate system, PG = LGl + RGl*Pl:\\n', PG)"
+ ));
+ testTextResult(entry, 0, QString::fromUtf8(
+ "Translation vector:\n"
+ " [[1]\n"
+ " [1]]\n"
+ "Rotation matrix:\n"
+ " [[ 0 -1]\n"
+ " [ 1 0]]\n"
+ "Position at the local coordinate system:\n"
+ " [[1]\n"
+ " [1]]\n"
+ "Position at the Global coordinate system, PG = LGl + RGl*Pl:\n"
+ " [[0]\n"
+ " [2]]"
+ ));
+ entry = entry->next();
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### Determination of the unknown angle of rotation\n"
+ "\n"
+ "If we didn't know the angle of rotation between the two coordinate systems, which is the typical situation in motion analysis, we simply would equate one of the terms of the two-dimensional rotation matrix in its algebraic form to its correspondent value in the numerical rotation matrix we calculated.\n"
+ "\n"
+ "For instance, taking the first term of the rotation matrices above: $\\cos\\alpha = 0$ implies that $\\theta$ is 90$^o$ or 270$^o$, but combining with another matrix term, $\\sin\\alpha = 1$, implies that $\\alpha=90^o$. We can solve this problem in one step using the tangent $(\\sin\\alpha/\\cos\\alpha)$ function with two terms of the rotation matrix and calculating the angle with the `arctan2(y, x)` function:"
+ "And this procedure would be repeated for each segment and for each instant of the analyzed movement to find the rotation of each segment."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "#### Joint angle as a sequence of rotations of adjacent segments\n"
+ "\n"
+ "In the notebook about [two-dimensional angular kinematics](http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/AngularKinematics2D.ipynb), we calculated segment and joint angles using simple trigonometric relations. We can also calculate these two-dimensional angles using what we learned here about the rotation matrix.\n"
+ "\n"
+ "The segment angle will be given by the matrix representing the rotation from the laboratory coordinate system (G) to a coordinate system attached to the segment and the joint angle will be given by the matrix representing the rotation from one segment coordinate system (l1) to the other segment coordinate system (l2). So, we have to calculate two basis now, one for each segment and the joint angle will be given by the product between the two rotation matrices. \n"
+ "\n"
+ "To define a two-dimensional basis, we need to calculate vectors perpendicular to each of these lines. Here is a way of doing that. First, let's find three non-collinear points for each basis:"
+ "Same result as obtained in [Angular kinematics in a plane (2D)](http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/AngularKinematics2D.ipynb). "
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "### Kinematic chain in a plain (2D)\n"
+ "\n"
+ "The fact that we simply multiplied the rotation matrices to calculate the rotation matrix of one segment in relation to the other is powerful and can be generalized for any number of segments: given a serial kinematic chain with links 1, 2, ..., n and 0 is the base/laboratory, the rotation matrix between the base and last link is: $\\mathbf{R_{n,n-1}R_{n-1,n-2} \\dots R_{2,1}R_{1,0}}$, where each matrix in this product (calculated from right to left) is the rotation of one link with respect to the next one. \n"
+ "\n"
+ "For instance, consider a kinematic chain with two links, the link 1 is rotated by $\\alpha_1$ with respect to the base (0) and the link 2 is rotated by $\\alpha_2$ with respect to the link 1. \n"
+ "Using Sympy, the rotation matrices for link 2 w.r.t. link 1 $(R_{12})$ and for link 1 w.r.t. base 0 $(R_{01})$ are: "
+ ));
+
+ qDebug() << "command entry 16";
+ testCommandEntry(entry, 16, QString::fromUtf8(
+ "from IPython.display import display, Math\n"
+ "from sympy import sin, cos, Matrix, simplify, latex, symbols\n"
+ "The typical use of all these concepts is in the three-dimensional motion analysis where we will have to deal with angles in different planes, which needs a special manipulation as we will see next."
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## Problems\n"
+ "\n"
+ "1. A local coordinate system is rotated 30$^o$ clockwise in relation to the Global reference system. \n"
+ " A. Determine the matrices for rotating one coordinate system to another (two-dimensional). \n"
+ " B. What are the coordinates of the point [1, 1] (local coordinate system) at the global coordinate system? \n"
+ " C. And if this point is at the Global coordinate system and we want the coordinates at the local coordinate system? \n"
+ " D. Consider that the local coordinate system, besides the rotation is also translated by [2, 2]. What are the matrices for rotation, translation, and transformation from one coordinate system to another (two-dimensional)? \n"
+ " E. Repeat B and C considering this translation.\n"
+ " \n"
+ "2. Consider a local coordinate system U rotated 45$^o$ clockwise in relation to the Global reference system and another local coordinate system V rotated 45$^o$ clockwise in relation to the local reference system U. \n"
+ " A. Determine the rotation matrices of all possible transformations between the coordinate systems. \n"
+ " B. For the point [1, 1] in the coordinate system U, what are its coordinates in coordinate system V and in the Global coordinate system? \n"
+ " \n"
+ "3. Using the rotation matrix, deduce the new coordinates of a square figure with coordinates [0, 0], [1, 0], [1, 1], and [0, 1] when rotated by 0$^o$, 45$^o$, 90$^o$, 135$^o$, and 180$^o$ (always clockwise).\n"
+ " \n"
+ "4. Solve the problem 2 of [Angular kinematics in a plane (2D)](http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/AngularKinematics2D.ipynb) but now using the concept of two-dimensional transformations. "
+ ));
+
+ testMarkdown(entry, QString::fromUtf8(
+ "## References\n"
+ "\n"
+ "- Robertson G, Caldwell G, Hamill J, Kamen G (2013) [Research Methods in Biomechanics](http://books.google.com.br/books?id=gRn8AAAAQBAJ). 2nd Edition. Human Kinetics. \n"
+ "- Ruina A, Rudra P (2013) [Introduction to Statics and Dynamics](http://ruina.tam.cornell.edu/Book/index.html). Oxford University Press. \n"
+ "- Winter DA (2009) [Biomechanics and motor control of human movement](http://books.google.com.br/books?id=_bFHL08IWfwC). 4 ed. Hoboken, EUA: Wiley. \n"
+ "- Zatsiorsky VM (1997) [Kinematics of Human Motion](http://books.google.com.br/books/about/Kinematics_of_Human_Motion.html?id=Pql_xXdbrMcC&redir_esc=y). Champaign, Human Kinetics."
KMessageBox::information(worksheetView(), i18n("In read-only mode Cantor couldn't guarantee, that the export will be valid for %1", m_backendName), i18n("Cantor"));
KMessageBox::information(worksheetView(), i18n("%1 backend was not found. Editing and executing entries is not possible", m_backendName), i18n("Cantor"));
m_readOnly = true;
}
else
m_readOnly = false;
if(!m_readOnly && !b->isEnabled())
{
QApplication::restoreOverrideCursor();
KMessageBox::information(worksheetView(), i18n("There are some problems with the %1 backend,\n"\
"please check your configuration or install the needed packages.\n"
"You will only be able to view this worksheet.", m_backendName), i18n("Cantor"));
m_readOnly = true;
}
if (m_readOnly)
{
// TODO: Handle this here?
for (QAction* action : m_richTextActionList)
action->setEnabled(false);
}
m_isLoadingFromFile = true;
//cleanup the worksheet and all it contains
delete m_session;
m_session=nullptr;
m_loginDone = false;
//file can only be loaded in a worksheet that was not eidted/modified yet (s.a. CantorShell::load())
//in this case on the default "first entry" is available -> delete it.
+ KMessageBox::error(worksheetView(), i18n("The file is old Jupyter notebook (found version %1.%2), which isn't supported by Cantor",nbformatMajor, nbformatMinor ), i18n("Cantor"));
+ if (QT_VERSION_CHECK(nbformatMajor, nbformatMinor, 0) < QT_VERSION_CHECK(4,0,0))
+ {
+ QApplication::restoreOverrideCursor();
+ KMessageBox::error(worksheetView(), i18n("Cantor doesn't support import Jupyter notebooks with version lower that 4.0 (detected %1.%2)",nbformatMajor, nbformatMinor), i18n("Cantor"));
+ KMessageBox::information(worksheetView(), i18n("%1 backend was not found. Editing and executing entries is not possible", m_backendName), i18n("Cantor"));
+ m_readOnly = true;
+ }
+ else
+ m_readOnly = false;
+
+ if(!m_readOnly && !backend->isEnabled())
+ {
+ QApplication::restoreOverrideCursor();
+ KMessageBox::information(worksheetView(), i18n("There are some problems with the %1 backend,\n"\
+ "please check your configuration or install the needed packages.\n"
+ "You will only be able to view this worksheet.", m_backendName), i18n("Cantor"));
+ m_readOnly = true;
+ }
+
+ if (m_readOnly)
+ {
+ for (QAction* action : m_richTextActionList)
+ action->setEnabled(false);
+ }
+
+
+ m_isLoadingFromFile = true;
+
+ if (m_session)
+ delete m_session;
+ m_session = nullptr;
+ m_loginDone = false;
+
+ if (m_firstEntry) {
+ delete m_firstEntry;
+ m_firstEntry = nullptr;
+ }
+
+ resetEntryCursor();
+
+ if (!m_readOnly)
+ m_session=backend->createSession();
+
+ qDebug() << "loading jupyter entries";
+
+ WorksheetEntry* entry = nullptr;
+ for (QJsonArray::const_iterator iter = cells.begin(); iter != cells.end(); iter++) {
+ if (!JupyterUtils::isJupyterCell(*iter))
+ {
+ QApplication::restoreOverrideCursor();
+ QString explanation;
+ if (iter->isObject())
+ explanation = i18n("an object with keys: %1", iter->toObject().keys().join(QLatin1String(", ")));
+ else
+ explanation = i18n("non object JSON value");
+
+ m_isLoadingFromFile = false;
+ showInvalidNotebookSchemeError(i18n("found incorrect data (%1) that is not Jupyter cell", explanation));