06_BestPractice_Maintenance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/JamesFergusson/Introduction-to-Research-Computing/blob/master/06_BestPractice_Maintenance.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Good practice for coding\n",
    "\n",
    "Now we know how to control the development of our code we need to think about how to write the code itself.\n",
    "\n",
    "The main goal of this section it to help you to write code that is readable and maintainable.  There is a well known addage which is: \"The older code gets the more it starts to stink\".  When we start projects it is always tempting to just get on with it.  Initally we know what we are doing and while we are starting out the code is fairly short and usually only does one or two things.  The problems only emerge a year or two down the line.  By now you will have modified your routines several times and added multiple new features.  The code will be long, messy and will likely contain redundant routines you are keeping \"just in case\".  Most likely there will be multiple sections that only work properly if you get the input just right.  Now editing the code has become a nightmare.  The code is too large to keep in your head so you can't quite remember exactly how every part works and which parts depend on each other.  Hence when you edit it you tend to introduce bugs that are then hard to track down.  Also you may wish to share the code with a collaborator so they can use for some new project but you have to spend a significant portion of every day trying to explain how it works and what bits they can edit safely.  The following ideas are all about stopping this happening.  \n",
    "\n",
    "There are three main things we will focus on \n",
    "- Documentation\n",
    "- Modularity\n",
    "- Prototyping\n",
    "\n",
    "The first two ensure readablity of your code so that others (and yourself) can follow what it does.  One of the great strengths of python is that its syntax closely follows normal English so it is easily readable, your code should be the same.\n",
    "\n",
    "## Documentation\n",
    "\n",
    "The first, and most important, topic in this section is documentation.  When you write code you should always so that it is clear exactly what it is that you are doing and why.  The goal for good documentation should be this:\n",
    "\n",
    "<b>\"If I die tomorrow, no-one should care\"</b>\n",
    "\n",
    "In essence your code should be so easy to follow anyone could read and understand it without any help from you.  Let's look at functions first.\n",
    "\n",
    "Python has something called the 'Docstring' which is a comment immediatly after the `def` line and is what is returned with `functionname?`.  This should be a short description of the function and how to use it.  Here is a popular format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def harmonic_log(n,x):\n",
    "    \"\"\"\n",
    "    Harmonic Logarithm\n",
    "    \n",
    "    Calculate and return the order 1 harmonic logarithm\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    n: int\n",
    "        the degree, must be positive\n",
    "    x: float\n",
    "        must be positive\n",
    "        \n",
    "    Returns\n",
    "    ----------\n",
    "    float\n",
    "        the harmonic logarithm of x with degree n\n",
    "    \"\"\"\n",
    "    harm_num = sum([1.0/i for i in range(1,n)])\n",
    "    return x**n*(math.log(x) - harm_num)\n",
    "\n",
    "harmonic_log?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should always do this when creating your own functions.  Use complete sentences and describe the function, what the input should be, and what it returns. Docstrings should also exist for modules. For packages they can be included in the `__init__.py` file.  The module/package docstring should describe its purpose and should list all classes, exceptions and functions that are exported from it. Here is the docstring in `__init__.py` for numpy as an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "NumPy\n",
    "=====\n",
    "\n",
    "Provides\n",
    "  1. An array object of arbitrary homogeneous items\n",
    "  2. Fast mathematical operations over arrays\n",
    "  3. Linear Algebra, Fourier Transforms, Random Number Generation\n",
    "\n",
    "How to use the documentation\n",
    "----------------------------\n",
    "Documentation is available in two forms: docstrings provided\n",
    "with the code, and a loose standing reference guide, available from\n",
    "`the NumPy homepage <http://www.scipy.org>`_.\n",
    "\n",
    "We recommend exploring the docstrings using\n",
    "`IPython <http://ipython.scipy.org>`_, an advanced Python shell with\n",
    "TAB-completion and introspection capabilities.  See below for further\n",
    "instructions.\n",
    "\n",
    "The docstring examples assume that `numpy` has been imported as `np`::\n",
    "\n",
    "  >>> import numpy as np\n",
    "\n",
    "Code snippets are indicated by three greater-than signs::\n",
    "\n",
    "  >>> x = 42\n",
    "  >>> x = x + 1\n",
    "\n",
    "Use the built-in ``help`` function to view a function's docstring::\n",
    "\n",
    "  >>> help(np.sort)\n",
    "  ... # doctest: +SKIP\n",
    "\n",
    "For some objects, ``np.info(obj)`` may provide additional help.  This is\n",
    "particularly true if you see the line \"Help on ufunc object:\" at the top\n",
    "of the help() page.  Ufuncs are implemented in C, not Python, for speed.\n",
    "The native Python help() does not know how to view their help, but our\n",
    "np.info() function does.\n",
    "\n",
    "To search for documents containing a keyword, do::\n",
    "\n",
    "  >>> np.lookfor('keyword')\n",
    "  ... # doctest: +SKIP\n",
    "\n",
    "General-purpose documents like a glossary and help on the basic concepts\n",
    "of numpy are available under the ``doc`` sub-module::\n",
    "\n",
    "  >>> from numpy import doc\n",
    "  >>> help(doc)\n",
    "  ... # doctest: +SKIP\n",
    "\n",
    "Available subpackages\n",
    "---------------------\n",
    "doc\n",
    "    Topical documentation on broadcasting, indexing, etc.\n",
    "lib\n",
    "    Basic functions used by several sub-packages.\n",
    "random\n",
    "    Core Random Tools\n",
    "linalg\n",
    "    Core Linear Algebra Tools\n",
    "fft\n",
    "    Core FFT routines\n",
    "polynomial\n",
    "    Polynomial tools\n",
    "testing\n",
    "    NumPy testing tools\n",
    "f2py\n",
    "    Fortran to Python Interface Generator.\n",
    "distutils\n",
    "    Enhancements to distutils with support for\n",
    "    Fortran compilers support and more.\n",
    "\n",
    "Utilities\n",
    "---------\n",
    "test\n",
    "    Run numpy unittests\n",
    "show_config\n",
    "    Show numpy build configuration\n",
    "dual\n",
    "    Overwrite certain functions with high-performance Scipy tools\n",
    "matlib\n",
    "    Make everything matrices.\n",
    "__version__\n",
    "    NumPy version string\n",
    "\n",
    "Viewing documentation using IPython\n",
    "-----------------------------------\n",
    "Start IPython with the NumPy profile (``ipython -p numpy``), which will\n",
    "import `numpy` under the alias `np`.  Then, use the ``cpaste`` command to\n",
    "paste examples into the shell.  To see which functions are available in\n",
    "`numpy`, type ``np.<TAB>`` (where ``<TAB>`` refers to the TAB key), or use\n",
    "``np.*cos*?<ENTER>`` (where ``<ENTER>`` refers to the ENTER key) to narrow\n",
    "down the list.  To view the docstring for a function, use\n",
    "``np.cos?<ENTER>`` (to view the docstring) and ``np.cos??<ENTER>`` (to view\n",
    "the source code).\n",
    "\n",
    "Copies vs. in-place operation\n",
    "-----------------------------\n",
    "Most of the functions in `numpy` return a copy of the array argument\n",
    "(e.g., `np.sort`).  In-place versions of these functions are often\n",
    "available as array methods, i.e. ``x = np.array([1,2,3]); x.sort()``.\n",
    "Exceptions to this rule are documented.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Naming \n",
    "The next thing is to spend a bit of time deciding on naming conventions for your code to make sure that the variable/function/class/etc.. names are sensible and explanatory.  For example compare the following two routines:   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sqb(x,e):\n",
    "    l=0.0\n",
    "    h=max(1.0,x)\n",
    "    a=0.5*(h+l)\n",
    "    while abs(a*a-x)>=e:\n",
    "        if a*a<x:\n",
    "            l=a\n",
    "        else:\n",
    "            h=a\n",
    "        a=0.5*(h+l)\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def square_root_bisection(x,tolerance):\n",
    "    \n",
    "    low_guess=0.0\n",
    "    high_guess=max(1.0,x)\n",
    "    \n",
    "    guess=0.5*(high_guess+low_guess)\n",
    "    \n",
    "    while abs(guess**2 - x) >= tolerance:\n",
    "        if guess**2 < x:\n",
    "            low_guess = guess\n",
    "        else:\n",
    "            high_guess = guess\n",
    "        guess = 0.5*(high_guess+low_guess)\n",
    "        \n",
    "    return guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sqb(25,0.01))\n",
    "print(square_root_bisection(25,0.01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are some basic conventions for naming things in python.  `i,j,k` are dummy integers, `x,y,z` are dummy floats, variables and functions are named in lower case with underscores (i.e. `obj_count` or `matrix_inv`), classes use camelcase (i.e. PointParticle).  You should avoid using `l`, `I` or `O` as they look like `1` and `0` in some fonts. You can't start a variable with a number and if some variables are related then you could give them the same prefix (e.g `loop_start, loop_end, loop_count`).  You should also use white space to make things clearer.  Generally put space around binary operators and space around lower priority operations e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = x*2 - 1\n",
    "hypot2 = x*x + y*y\n",
    "c = (a+b) * (a-b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next you should also include comments to describe what each code section does, again use complete sentences.  Comments begin with a `#` and in general should be on the preceding line. For the `harmonic_log` function we could have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def harmonic_log(n,x):\n",
    "    \"\"\"\n",
    "    Harmonic Logarithm\n",
    "    \n",
    "    Calculate and return the order 1 harmonic logarithm\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    n: int\n",
    "        the degree, must be positive\n",
    "    x: float\n",
    "        must be positive\n",
    "        \n",
    "    Returns\n",
    "    ----------\n",
    "    float\n",
    "        the harmonic logarithm of x with degree n\n",
    "    \"\"\"\n",
    "    \n",
    "    # Calculate the harmonic number H_n\n",
    "    harm_num = sum([1.0/i for i in range(1,n)])\n",
    "    \n",
    "    #return the harmonic logarithm \n",
    "    return x**n*(math.log(x) - harm_num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a little bit overkill. Generally docstrings should explain how to use funtions and modules, comments should explain what sections of code do and everything else should be named to that it is clear what it is. Don't be afraid to be verbose, your goal it to make it as simple as possible for someone unfamiliar with the code to understand it at first reading and with pythons auto-complete you can't complain that they take too long to type."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modularity\n",
    "\n",
    "The next thing you should do to make your code readable is to break it up into discreate blocks each of which excecutes one 'idea' by moving detailed calculation into functions and modules.  We met this idea when we first met modules earlier.  This makes the code more readable as you can see its structure without getting bogged down in detail.  It also generally makes it more stable and easier to edit as you can work on discrete blocks in isolation.  Finally it allows code reuse.  Once you have written a module that say, reads and processes your particular data set, you can use it everywhere and if you then find a bug in it you only have to fix it in one place.\n",
    "\n",
    "It is tricky to write a guide as to the best way to modulise your code but here are three main rules to guide you:\n",
    "\n",
    "1. Each routine should to be not more than one screen long.\n",
    "2. Any code that you will use in more than one place should be moved to a seperate routine.\n",
    "3. Routines should only have one main task.\n",
    "\n",
    "The first makes your code more readable and makes it harder to get 'lost'. It also makes editing safer as each part is understandable on its own.  The second not only saves time but means that if you want to change something you only need to do it once.  The third means that you only have to worry about what else you might be breaking when you edit something.  All of these rules can be broken if it makes your code simpler or clearer. There is also a risk that if code is very heavily modulised then it can be hard to find where calculations actually happen and can even slow down your code.\n",
    "\n",
    "\n",
    "## Prototyping\n",
    "\n",
    "One other important coding practice is prototyping.  Prototyping is simply deciding what you are going to do, then checking that your ideas work, before launching into writing the detailed code for your problem.  The temptation is always to \"just get to work\" and start coding.  The problem with this approach is that once you start coding in a particular direction when you meet issues it can be hard to go back and fix our original bad decisions rather than just finding a workaround for the particular issue.  The analogy is if you start building a car but discover halfday through that where you want to drive is very muddy. So you go back and add bigger grippier tyres but then you find the car doesn't have the tourqe to turn them. To solve this you install a low range gearbox but then the engine can't handle it so you add a second engine to boost power...   At the end of this you will end up with a weird mish-mash car when really you should have gone back to the start and made a tractor instead. \n",
    "\n",
    "The best place to start prototyping is with a pen and paper.  As almost all programming we are going to do is related to data processing the first thing to do is to think about what data objects you are going to need starting from input data and the final results you want to output. Next you have to decide how to get between the two in the most efficent manner and what data objects you will need on the way.  With this we can slowly build up a plan for the code.  Here is a list of things you should consider when you prototype:\n",
    "\n",
    "1. What is the size of the problem? Do I need to worry about memory/parallelisation/architecture constraints for this code?\n",
    "2. What does my input data look like and what data do I want to output?\n",
    "3. What logical operations do I need to perform on the data?\n",
    "4. Can I map the operations onto standard optimised routines from numpy/scipy/etc...?\n",
    "5. What temporary data objects will I need to create?\n",
    "6. What is a sensible naming convention for my data objects/variables? (more important than it sounds)\n",
    "7. How should I modularise this, what might I want to reuse from here in other future codes?\n",
    "8. What should the inputs and outputs of each module be?\n",
    "\n",
    "Try to be reasonably detailed, the goal is to identify roadblocks with your method and if they were obvious from a rough plan you would have already thought of them.\n",
    "\n",
    "The next step is to make small models to test parts of the code to make sure your ideas work.  Our tests should check the method works, that is is accurate and also what the memory and CPU requirments wiil be for the full problem size.  Not everything needs this kind of check but you should try to do it with any section you are unsure about.  \n",
    "\n",
    "So far this has been very theoretical, let's try a worked example to see the ideas in action.\n",
    "\n",
    "### Example\n",
    "Suppose you want to calculate a multidimensional sum over some combination of functions, $X^i_{l}$, but with some particular constraints:\n",
    "\n",
    "$E_{ijk} = \\sum_{l_1}\\sum_{l_2}\\sum_{l_3} \\left( \\int_{-1}^{1} d\\mu P_{l_1}(\\mu) P_{l_2}(\\mu) P_{l_3}(\\mu)  \\right) X^i (l_1) X^j (l_2) X^k (l_3) $\n",
    "\n",
    "Here in the triple sum the integral over Legendre polynomials enforces $l_1+l_2+l_3$ is even (parity) and that no $l_i$ is larger than the sum of the other two (conservation of angular momentum).  The integral representation lets us re-factor the problem as:\n",
    "\n",
    "$E_{ijk} = \\int_{-1}^{1} d\\mu Y^i(\\mu)Y^j(\\mu)Y^k(\\mu)$\n",
    "\n",
    "$Y^i(\\mu) = \\sum_l P_{l}(\\mu)  X^i (l)$\n",
    "\n",
    "So our inputs are $X^i_{l}$, and we want to output $E_{ijk}$. Let's take $X^i_l$ to be Chebyshev polynomials of the first kind, $T_i(l)$, as an example which we will generate by recursion. We could just as easily read them from a file or whatever, they're not important.  \n",
    "\n",
    "![](Plots/Proto1.png)\n",
    "\n",
    "The next steps are to fill in the operations, a sum and integration\n",
    "\n",
    "![](Plots/Proto2.png)\n",
    "\n",
    "How will we do the operations? The sum is easy, this is a matrix multiplication so we can use numpy `np.dot` which will be efficent, we will also need to generate the Legendre polynomials in a matrix, recursion is best for this as iut's both fast and accurate.  For the integral we want to stick with the array form so the simple trapazoidal rule sounds pretty good as a first step\n",
    "\n",
    "![](Plots/Proto3.png)\n",
    "\n",
    "The only thing is the integral over $\\mu$ is a bit tricky, we need it to be zero for both $l_1+l_2+l_3$ odd and is one $l_i$ is larger than the sum of the other two.  This has to happen <b>after</b> we sum over $l$ so there is a real chance of numerical issues here.  Let's do a quick check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Test\n",
    "import numpy as np\n",
    "import scipy.integrate as sp\n",
    "\n",
    "mu = 50\n",
    "lmax = 5\n",
    "P = np.zeros((lmax,mu), dtype='float')\n",
    "\n",
    "X = np.linspace(-1,1,mu)\n",
    "\n",
    "P[0,:]=1.0\n",
    "P[1,:]=X\n",
    "for n in range(2,lmax):\n",
    "    P[n,:] = ((2*n-1)*X*P[n-1,:] - (n-1)*P[n-2,:])/n\n",
    "\n",
    "PPP1 = P[2,:]*P[2,:]*P[2,:] # 4/35\n",
    "PPP2 = P[0,:]*P[1,:]*P[3,:] # 0\n",
    "PPP3 = P[0,:]*P[1,:]*P[2,:] # 0\n",
    "PPP4 = P[2,:]*P[4,:]*P[4,:] # 40/693\n",
    "print(sp.trapz(PPP1,X),' vs ',4/35)\n",
    "print(sp.trapz(PPP2,X),' vs ',0)\n",
    "print(sp.trapz(PPP3,X),' vs ',0)\n",
    "print(sp.trapz(PPP4,X),' vs ',40/693)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmmm even with $l$'s as low as 4 using 50 points is pretty bad.  We could increase the number of points more but as it's a polynomial integral we should probably use Gauss-Legendre quadriture as that is exact with (degree-1)/2 points which here would only be 6 points (much better than 50).  Let's check it works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "lmax = 5\n",
    "npts = 3*lmax//2+1\n",
    "\n",
    "X, W = np.polynomial.legendre.leggauss(npts)\n",
    "P = np.zeros((lmax,npts), dtype='float')\n",
    "\n",
    "P[0,:]=1.0\n",
    "P[1,:]=X\n",
    "for n in range(2,lmax):\n",
    "    P[n,:] = ((2*n-1)*X*P[n-1,:] - (n-1)*P[n-2,:])/n\n",
    "\n",
    "PPP1 = W[:]*P[2,:]*P[2,:]*P[2,:] # 4/35\n",
    "PPP2 = W[:]*P[0,:]*P[1,:]*P[3,:] # 0\n",
    "PPP3 = W[:]*P[0,:]*P[1,:]*P[2,:] # 0\n",
    "PPP4 = W[:]*P[2,:]*P[4,:]*P[4,:] # 40/693\n",
    "print(np.sum(PPP1),' vs ',4/35)\n",
    "print(np.sum(PPP2),' vs ',0)\n",
    "print(np.sum(PPP3),' vs ',0)\n",
    "print(np.sum(PPP4),' vs ',40/693)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks much better! (and faster!).  Now our plan looks like this:\n",
    "\n",
    "![](Plots/Proto4.png)\n",
    "\n",
    "Now we can write our code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from itertools import combinations_with_replacement\n",
    "lmax = 50\n",
    "imax = 5\n",
    "npts = 3*lmax//2+1\n",
    "\n",
    "X, W = np.polynomial.legendre.leggauss(npts)\n",
    "P = np.zeros((lmax,npts), dtype='float')\n",
    "C = np.zeros((imax,lmax), dtype='float')\n",
    "\n",
    "\n",
    "P[0,:]=1.0\n",
    "P[1,:]=X\n",
    "for n in range(2,lmax):\n",
    "    P[n,:] = ((2*n-1)*X*P[n-1,:] - (n-1)*P[n-2,:])/n\n",
    "\n",
    "\n",
    "L = np.linspace(-1,1,lmax)\n",
    "\n",
    "C[0,:]=1.0\n",
    "C[1,:]=L\n",
    "for n in range(2,imax):\n",
    "    C[n,:] = 2*L*C[n-1,:] - C[n-2,:]\n",
    "\n",
    "T = np.dot(C,P)\n",
    "\n",
    "comb = combinations_with_replacement(range(imax),3)\n",
    "\n",
    "for i,j,k in comb:\n",
    "    z = np.sum(W[:]*T[i,:]*T[j,:]*T[k,:])\n",
    "    print(\"({},{},{}) produces {:5g}\".format(i,j,k,z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's pretty short so there is not much need for functions, however we may want to move our Legendre and Chebyshev polynomial generators to a module so we can use them elsewhere. Some documentation would also be great.\n",
    "\n",
    "Here is the final module file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "This Module contains functions to create arrays of common polynomials\n",
    "\n",
    "Contains:\n",
    "----------------------------------------\n",
    "    legendre_array\n",
    "        Creates an array of the legendre polynomials\n",
    "    chebyshev_array\n",
    "        Creates an array of the chebyshev polynomials\n",
    "----------------------------------------\n",
    "        \n",
    "        \n",
    "All functions take two arguments:\n",
    "----------------------------------------\n",
    "    lmax\n",
    "        the maximum dregree to include\n",
    "    X\n",
    "        The array of X points to calculate the polynomials on\n",
    "----------------------------------------\n",
    "\n",
    "Calculation is done by recursion\n",
    "\n",
    "Written by James Fergusson: J.Fergusson@DAMTP.cam.ac.uk \n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "def legendre_array(X,lmax):\n",
    "    \"\"\"\n",
    "    Calculates an array of Legendre polynomials\n",
    "    with size: len(X) by lmax\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lmax: int\n",
    "        the maximum degree, must be positive\n",
    "    X: array(float)\n",
    "        points to calculate the polynomials on\n",
    "        \n",
    "    Returns\n",
    "    ----------\n",
    "    array(float)\n",
    "        A table of the Legendre polynomials\n",
    "    \"\"\"\n",
    "    \n",
    "    # Initilise the array\n",
    "    P = np.zeros((lmax,len(X)), dtype='float')\n",
    "    \n",
    "    # Set P_0 and P_l\n",
    "    P[0,:]=1.0\n",
    "    P[1,:]=X\n",
    "    \n",
    "    # Calculate the rest via iteration\n",
    "    for n in range(2,lmax):\n",
    "        P[n,:] = ((2*n-1)*X*P[n-1,:] - (n-1)*P[n-2,:])/n\n",
    "    \n",
    "    # return the array of Legendre polynomials\n",
    "    return P\n",
    "\n",
    "def chebyshev_array(X,lmax):\n",
    "    \"\"\"\n",
    "    Calculates an array of Chebyshev polynomials\n",
    "    with size: len(X) by lmax\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lmax: int\n",
    "        the maximum degree, must be positive\n",
    "    X: array(float)\n",
    "        points to calculate the polynomials on\n",
    "        \n",
    "    Returns\n",
    "    ----------\n",
    "    array(float)\n",
    "        A table of the Chebyshev polynomials\n",
    "    \"\"\"\n",
    "    \n",
    "    # Initilise the array\n",
    "    C = np.zeros((lmax,len(X)), dtype='float')\n",
    "    \n",
    "    # Set C_0 and C_l\n",
    "    C[0,:]=1.0\n",
    "    C[1,:]=X\n",
    "    \n",
    "    # Calculate the rest via iteration\n",
    "    for n in range(2,lmax):\n",
    "        C[n,:] = 2*X*C[n-1,:] - C[n-2,:]\n",
    "     \n",
    "    # return the array of Chebyshev polynomials   \n",
    "    return C\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here is the main code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Code in triangle_int.py  \n",
    "\"\"\"\n",
    " This code computes the value of the integral:                        \n",
    "     $E_{i,j,k} = \\sum_{l_1}\\sum_{l_2}\\sum_{l_3}                      \n",
    "         \\left( \\int_{-1}^{1} d\\mu P_{l_1}(\\mu) P_{l_2}(\\mu)          \n",
    "         P_{l_3}(\\mu)\\right) X^i(l_1) X^j(l_2) X^k(l_3) $          \n",
    "                                                                      \n",
    " The code takes some set of $X^i(l)$ for i element of [0,imax]       \n",
    " and calculates $E_{i,j,k}$ for every possible triple of [0,imax]     \n",
    "                                                                      \n",
    " The code works by reordering the integration and sum to form:        \n",
    "     $Y^i(\\mu) = \\sum_l P_{l}(\\mu)  X^i(l)$                          \n",
    "     $E_{i,j,k} = \\int_{-1}^{1} d\\mu Y^i(\\mu)Y^j(\\mu)Y^k(\\mu)$       \n",
    "                                                                      \n",
    " The first is computed as a matix product                             \n",
    " The second is computed using Gauss Legendre quadriture               \n",
    "     as it is exact for polynomial arguments                          \n",
    "                                                                      \n",
    " Written by James Fergusson: J.Fergusson@DAMTP.cam.ac.uk              \n",
    "\"\"\"\n",
    "\n",
    "# Load modules:\n",
    "import numpy as np\n",
    "from itertools import combinations_with_replacement\n",
    "from Tools import polynomials as poly\n",
    "\n",
    "# Set maximums for l and i.\n",
    "lmax = 50\n",
    "imax = 5\n",
    "\n",
    "# Compute number of points requiredfor the Gauss-Legendre \n",
    "# quadtiture to be exact for a given lmax.\n",
    "npts = 3*lmax//2+1\n",
    "\n",
    "# Calculate Gauss-Legendre quadriture points, X, and weights, W.\n",
    "X, W = np.polynomial.legendre.leggauss(npts)\n",
    "\n",
    "# Create arrays for Legenre polynomials\n",
    "P = poly.legendre_array(X,lmax)\n",
    "\n",
    "    \n",
    "# Create Chebyshev Polynomials \n",
    "# This could be changed to accept a general function or read from disk\n",
    "# in future versions\n",
    "L = np.linspace(-1,1,lmax)\n",
    "C = poly.chebyshev_array(L,imax)\n",
    "\n",
    "# Compute the sum over L using a matrix product\n",
    "Y = np.dot(C,P)\n",
    "\n",
    "# Construct an iterator over possible triples of i\n",
    "comb = combinations_with_replacement(range(imax),3)\n",
    "\n",
    "# Calculate the integral (just a sum here) for each triple and print \n",
    "# the result \n",
    "for i,j,k in comb:\n",
    "    z = np.sum(W[:]*Y[i,:]*Y[j,:]*Y[k,:])\n",
    "    print(\"({},{},{}) produces {:5g}\".format(i,j,k,z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Now it's your turn. I want you to create a programme that can take the noisy data in the files `Data/period?.txt` extract the periodic features and stack them up on a plot.  You will want to use routines in the package scipy signal: https://docs.scipy.org/doc/scipy/reference/signal.html. You should plan you strategy on paper then test and implement it with proper ducumentation.  Once you have finished you will swap places with someone else and try to understand each others code without any explaination. The result should take the left and turn it into the right in the figues below:\n",
    "\n",
    "![](Plots/Period1.png)\n",
    "![](Plots/Period2.png)\n",
    "![](Plots/Period3.png)\n",
    "![](Plots/Period4.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally I will note that when you are starting with python it can be helpful to use \"style checkers\" which check if your code confirms to the standard style guide for python.  The most used one is `pylint` and you can learn about it here:\n",
    "\n",
    "https://pylint.readthedocs.io/en/latest/index.html\n",
    "\n",
    "It will analyse python files and suggest stylistic changes which can help make the code clearer."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}