{
 "metadata": {
  "name": "Mori_LPA_1D"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# AUTO 07p in the IPython Notebook"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Introduction"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I have used both the fantastic [XPP/XPPAUT](http://www.math.pitt.edu/~bard/xpp/xpp.html) software package and the [matcont](http://www.matcont.ugent.be/) MATLAB toolbox\n",
      "for numerical continuation and bifurcation analysis of ordinary differential equations.\n",
      "\n",
      "As a user, both of these packages have you do mostly GUI operations: clicking buttons and typing keyboard shortcuts -- although there may be ways of automating usage of both packages.\n",
      "\n",
      "The software, [AUTO](http://indy.cs.concordia.ca/auto/), that XPPAUT relies on (the `auto` part of XPP used for continuation and bifurcation analysis) is a bit trickier to use\n",
      "than either one of the aforementioned packages -- although, in my very limited experience, using AUTO directly seems to be a more stable eperience than calling it through XPPAUT.\n",
      "\n",
      "While both XPP and matcont provide their own, user-friendly notation for defining ODEs and relatively straight-forward interfaces, AUTO requires the user to dig into some\n",
      "technical aspects of defining your ODE and setting AUTO's numerous parameters and variables.\n",
      "\n",
      "Numerous examples provided with the latest version of AUTO give newcomers to AUTO a reasonbly fast and pain-free entry point and you'll be up and running relatively quickly.\n",
      "\n",
      "Using AUTO with its current AUTO CLUI (a Python-based command-line tool) and the various methods defined in AUTO's Python interface has been a pleasurable (whatever I mean by that)\n",
      "experience but left me confused regarding some aspects:\n",
      "\n",
      "I generally found it hard to grasp what data structure the AUTO output had been saved in and how to plot the resultant continuation branches and bifurcation points in precisely\n",
      "the way I wanted.\n",
      "\n",
      "I started this notebook to just play around with the Python interface that ships with the latest AUTO version.\n",
      "\n",
      "For now, my goal is to parse and plot (using matplotlib) AUTO output. At some later stage I would love to find a way of integrating an ODE numerically and then passing the discovered \n",
      "stable steady state into AUTO (zeal abound).\n",
      "\n",
      "In the unlikely event that the person reading this is not me and you have any words of advice or criticism, I would be delighted to [hear from you](http://georg.io).\n",
      "\n",
      "This notebook [lives here](https://github.com/waltherg/auto/blob/master/examples/LPA_Mori_1D/Mori_LPA_1D.ipynb) and is part of a [repoistory all things auto](https://github.com/waltherg/auto/)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## This Notebook\n",
      "\n",
      "The structure of this notebook:\n",
      "\n",
      "* First I'll load a few Python modules (most of which I won't use for now) and then invoke AUTO for [a sample system](https://github.com/waltherg/auto/tree/master/examples/LPA_Mori_1D),\n",
      "\n",
      "* then I'll parse and plot some of AUTO's output manually (using `pandas`), and\n",
      "\n",
      "* at last I'll use some of the Python modules that were written by the AUTO developers that can be used to parse AUTO's output.\n",
      "\n",
      "Again, eveything I'm doing here can be done with the AUTO CLUI with just a few keystrokes -- here I'm just trying to start digging into parsing AUTO output and hopefully, eventually how AUTO works."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Load AUTO 07p Python Moduls"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Add `$AUTO_DIR/python` to `sys.path` to import the Python modules defined by `AUTO 07p`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sys\n",
      "auto_directory = !echo $AUTO_DIR\n",
      "if auto_directory == ['']:\n",
      "    home = !echo $HOME\n",
      "    auto_directory = home[0]+'/auto/07p/'\n",
      "    sys.path.append(auto_directory+'/python')\n",
      "else:\n",
      "    sys.path.append(auto_directory[0] + '/python')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Not all of these are necessary -- in fact we'll use just one of these modules for now --\n",
      "however it's good to see what Python modules have been written by the authors AUTO."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import AUTOCommands as ac\n",
      "import AUTOclui as acl\n",
      "import interactiveBindings as ib\n",
      "import runAUTO as ra"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Start AUTO, Load System, and Run"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Start AUTO and catch the returned `runner` object."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "runner = ra.runAUTO()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`AUTOCommands` defines multiple methods that return *solution* objects. Method `load` is one of them."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lpa = ac.load('lpa', runner=runner);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The above is shorthand for\n",
      "\n",
      "    ac.load(e='lpa', c='lpa', runner=runner)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Going through the library code, it is hard to understand what the following does exactly.\n",
      "\n",
      "For now, this compiles our system defined in `lpa.f90` into an object `lpa.o` and then links `lpa.o` and AUTO's FORTRAN library objects in `auto/07p/lib/` into an executable `lpa.exe`.\n",
      "\n",
      "After executing the next command, you'll notice this `lpa.exe` in your directory and you will be able to rerun AUTO just by doing this in your shell:\n",
      "\n",
      "    $ ./lpa.exe"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "lpa.run()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Going through the output of the above `lpa.run()` command you'll notice at least two things:**\n",
      "\n",
      "1. It would probably suffice to `os.system()` run the following commands (delete `lpa.o` in the current directory if you don't notice these command invocations at the top)\n",
      "\n",
      "```\n",
      "    gfortran -fopenmp -O -c lpa.f90 -o lpa.o\n",
      "    gfortran -fopenmp -O lpa.o -o lpa.exe $HOME/auto/07p/lib/*.o\n",
      "    ./lpa.exe\n",
      "```\n",
      "\n",
      "   By that I mean, you could probably just write a Python script that does\n",
      "\n",
      "    os.system('gfortran -fopenmp -O -c lpa.f90 -o lpa.o')\n",
      "    os.system('gfortran -fopenmp -O lpa.o -o lpa.exe $HOME/auto/07p/lib/*.o')\n",
      "    os.system('./lpa.exe')\n",
      "\n",
      "   This should generate the same output files in your directory.\n",
      "\n",
      "2. Some part of the `auto` toolchain seems broken by our use of the IPython Notebook\n",
      "\n",
      "```\n",
      "\n",
      "       596             raise AUTOExceptions.AUTORuntimeError(\"Error running AUTO\")\n",
      "\n",
      "       597 \n",
      "\n",
      "       598 def test():\n",
      "\n",
      "   AUTORuntimeError: Error running AUTO\n",
      "```\n",
      "\n",
      "**Despite the fact that something seems broken, we notice that `auto` still ran and created three output files.**"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Read Output Files and Plot Solution Branches\n",
      "\n",
      "### Parse and Plot By Hand"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!ls"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `b.lpa` file (see comment on naming schemes below) seems to hold all branches and some `auto`-specific info at the top of the file.\n",
      "\n",
      "Let's import a few tools that will come in handy reading `b.lpa` and plotting the branches described therein."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pandas as pd\n",
      "from matplotlib import pylab as pl"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Parse `b.lpa`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "content = None\n",
      "with open('b.lpa', 'r') as f:\n",
      "    content = f.readlines()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Line 15 (zero-indexed) onwards are the branches."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "content[15:20]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "content_csv = [[el for el in content[15].split(' ') if len(el) > 0 and el != '\\n']]\n",
      "content_csv[0][0] = 'branch'\n",
      "column_names = content_csv[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These are our column names found in `content[15]`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "column_names"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for line in content:\n",
      "    dummy = line.split(' ')\n",
      "    dummy = [el for el in dummy if len(el) > 0 and el != '\\n']\n",
      "    if dummy[0] == '0':\n",
      "        continue\n",
      "        \n",
      "    for el_i, el in enumerate(dummy):\n",
      "        if el_i < 4:\n",
      "            dummy[el_i] = int(el)\n",
      "        else:\n",
      "            dummy[el_i] = float(el)\n",
      "            \n",
      "    if len(dummy) > 1:\n",
      "        content_csv.append(dummy)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Load the resulting list of lists into a `pandas` `DataFrame`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df = pd.DataFrame(content_csv, columns=column_names)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Plot branches in `b.lpa`"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plotting branches is easy now."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "scatter(df[df['branch'] == 1].tot,df[df['branch'] == 1].aL)\n",
      "scatter(df[df['branch'] == 2].tot,df[df['branch'] == 2].aL)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Parse and Plot with Built-In AUTO CLUI methods"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Parse `b.lpa`\n",
      "\n",
      "The manual says that files `fort.7`, `fort.8` and `fort.9` are equivalent to `b.*`, `s.*`, and `d.*` (in that order).\n",
      "The `fort.*` naming scheme may be older as the latest version of AUTO appears to produce files in the latter naming scheme.\n",
      "\n",
      "`parseB.py` in `auto/07p/python/` is probably a good start for what we want to do here."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import parseB"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The class `parseB.parseBMixin` provides a method called `readFilename` so let's hope this allows us to read and parse our `b.lpa` file."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj = parseB.parseBMixin()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.readFilename('b.lpa')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Hmmm ... nope."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj = parseB.parseB()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "b_lpa = open('b.lpa', 'r')\n",
      "ab_obj.read(b_lpa)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.read(b_lpa)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The object `pb_obj` seems to store an array of the branches found in `b.lpa`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.branches[:2]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Branches are saved as instances of `Points.Pointset` which inherits from class `Points` (`Points.py` is found in `auto/07p/python`).\n",
      "\n",
      "Let's see what members the class `Pointset` defines."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.branches[0].name"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.branches[0].keys()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.branches[0]['tot']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.branches[0].todict()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The object `pb_obj.branches[0]` holds the coordinates for the first branch in `b.lpa` in an accessible format -- perfect for plotting.\n",
      "\n",
      "All we need now are stability properties and bifurcation points on this branch!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.branches[0].labels"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This branch is not very eventful: All we get are one endpoint at index 0 and another endpoint at index 46.\n",
      "\n",
      "Let's check out the next branch."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.branches[1].labels"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pb_obj.branches[1]['tot'][37]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This branch point `BP` (in this case, a transcritical bifurcation point) is one of four bifurcation points described in our [publication](http://link.springer.com/article/10.1007%2Fs11538-012-9766-5/fulltext.html) -- note that the `tot` parameter used in this notebook is defined an order of magnitude smaller than the same parameter (`tot` or `T` as we called it) in this publication.\n",
      "\n",
      "Hence, in this publcation we reported `T = 23.0` as bifurcation point whereas here we observe `tot = 2.3` -- these bifurcation points are equivalent."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We now have the coordinates of all branches in `b.lpa` and special points (bifurcation points) along these branches.\n",
      "\n",
      "What we still want is the stability of these branches. Inspecting our current directory, we notice a file `d.lpa` (this file may be called `fort.9` in the older naming scheme) which contains the eigenvalues at some points along these branches.\n",
      "\n",
      "A Python module `parseD.py` exists in `auto/07p/python` so let's see what this offers."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import parseD as pd"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pd_obj = pd.parseD()\n",
      "pd_obj.read(open('d.lpa', 'r'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The object `pd_obj` appears to be a list of dictionaries and while the overall parsing of `d.lpa` does not appear to be implemented to perfection (yet) we are given access to the eigenvalues and branch number (let's hope that branch numbers are consistent throughout!)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pd_obj[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(60,100):\n",
      "    print 'Branch number',pd_obj[i]['Branch number'],'eigenvalues',pd_obj[i]['Eigenvalues']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pd_obj[71]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It is not entirely clear to me how the diagnostic output in `d.lpa` and the branches in `b.lpa` can be combined into one bifurcation diagram."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import parseS as ps"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ps_obj = ps.parseS()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ps_obj.read(open('s.lpa', 'r'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ps_obj[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ps_obj[0].__str__()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ps_obj[0].data"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Stability Data and Bifurcation Points"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The above bifurcation plot is close to the one shown in [Figure 3](http://link.springer.com/static-content/images/354/art%253A10.1007%252Fs11538-012-9766-5/MediaObjects/11538_2012_9766_Fig3_HTML.gif) of a [previous article](http://link.springer.com/article/10.1007%2Fs11538-012-9766-5/fulltext.html) where the same analysis was done on the same system.\n",
      "\n",
      "What's still missing are stability data (eigenvalues along the solution branches) and bifurcation points.\n",
      "\n",
      "The file `d.lpa` generated by `AUTO` seems to contain that information."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dlpa = None\n",
      "with open('d.lpa', 'r') as f:\n",
      "    dlpa = f.readlines()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dlpa[:20]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}