I have used both the fantastic XPP/XPPAUT software package and the matcont MATLAB toolbox for numerical continuation and bifurcation analysis of ordinary differential equations.
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.
The software, AUTO, that XPPAUT relies on (the auto part of XPP used for continuation and bifurcation analysis) is a bit trickier to use
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.
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 technical aspects of defining your ODE and setting AUTO's numerous parameters and variables.
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.
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) experience but left me confused regarding some aspects:
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 the way I wanted.
I started this notebook to just play around with the Python interface that ships with the latest AUTO version.
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 stable steady state into AUTO (zeal abound).
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.
This notebook lives here and is part of a repoistory all things auto.
The structure of this notebook:
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,
then I'll parse and plot some of AUTO's output manually (using pandas), and
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.
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.
Add $AUTO_DIR/python to sys.path to import the Python modules defined by AUTO 07p.
In [ ]:
import sys
auto_directory = !echo $AUTO_DIR
if auto_directory == ['']:
home = !echo $HOME
auto_directory = home[0]+'/auto/07p/'
sys.path.append(auto_directory+'/python')
else:
sys.path.append(auto_directory[0] + '/python')
Not all of these are necessary -- in fact we'll use just one of these modules for now -- however it's good to see what Python modules have been written by the authors AUTO.
In [ ]:
import AUTOCommands as ac
import AUTOclui as acl
import interactiveBindings as ib
import runAUTO as ra
Start AUTO and catch the returned runner object.
In [ ]:
runner = ra.runAUTO()
AUTOCommands defines multiple methods that return solution objects. Method load is one of them.
In [ ]:
lpa = ac.load('lpa', runner=runner);
The above is shorthand for
ac.load(e='lpa', c='lpa', runner=runner)Going through the library code, it is hard to understand what the following does exactly.
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.
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:
$ ./lpa.exe
In [ ]:
lpa.run()
Going through the output of the above lpa.run() command you'll notice at least two things:
os.system() run the following commands (delete lpa.o in the current directory if you don't notice these command invocations at the top) gfortran -fopenmp -O -c lpa.f90 -o lpa.o
gfortran -fopenmp -O lpa.o -o lpa.exe $HOME/auto/07p/lib/*.o
./lpa.exeBy that I mean, you could probably just write a Python script that does
os.system('gfortran -fopenmp -O -c lpa.f90 -o lpa.o')
os.system('gfortran -fopenmp -O lpa.o -o lpa.exe $HOME/auto/07p/lib/*.o')
os.system('./lpa.exe')
This should generate the same output files in your directory.
auto toolchain seems broken by our use of the IPython Notebook
596 raise AUTOExceptions.AUTORuntimeError("Error running AUTO")
597
598 def test():
AUTORuntimeError: Error running AUTODespite the fact that something seems broken, we notice that auto still ran and created three output files.
In [ ]:
!ls
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.
Let's import a few tools that will come in handy reading b.lpa and plotting the branches described therein.
In [ ]:
import pandas as pd
from matplotlib import pylab as pl
In [ ]:
content = None
with open('b.lpa', 'r') as f:
content = f.readlines()
Line 15 (zero-indexed) onwards are the branches.
In [ ]:
content[15:20]
In [ ]:
content_csv = [[el for el in content[15].split(' ') if len(el) > 0 and el != '\n']]
content_csv[0][0] = 'branch'
column_names = content_csv[0]
These are our column names found in content[15].
In [ ]:
column_names
In [ ]:
for line in content:
dummy = line.split(' ')
dummy = [el for el in dummy if len(el) > 0 and el != '\n']
if dummy[0] == '0':
continue
for el_i, el in enumerate(dummy):
if el_i < 4:
dummy[el_i] = int(el)
else:
dummy[el_i] = float(el)
if len(dummy) > 1:
content_csv.append(dummy)
Load the resulting list of lists into a pandas DataFrame.
In [ ]:
df = pd.DataFrame(content_csv, columns=column_names)
In [ ]:
df
Plotting branches is easy now.
In [ ]:
scatter(df[df['branch'] == 1].tot,df[df['branch'] == 1].aL)
scatter(df[df['branch'] == 2].tot,df[df['branch'] == 2].aL)
b.lpaThe manual says that files fort.7, fort.8 and fort.9 are equivalent to b.*, s.*, and d.* (in that order).
The fort.* naming scheme may be older as the latest version of AUTO appears to produce files in the latter naming scheme.
parseB.py in auto/07p/python/ is probably a good start for what we want to do here.
In [ ]:
import parseB
The class parseB.parseBMixin provides a method called readFilename so let's hope this allows us to read and parse our b.lpa file.
In [ ]:
pb_obj = parseB.parseBMixin()
In [ ]:
pb_obj.readFilename('b.lpa')
Hmmm ... nope.
In [ ]:
pb_obj = parseB.parseB()
In [ ]:
b_lpa = open('b.lpa', 'r')
ab_obj.read(b_lpa)
In [ ]:
pb_obj.read(b_lpa)
The object pb_obj seems to store an array of the branches found in b.lpa.
In [ ]:
pb_obj.branches[:2]
Branches are saved as instances of Points.Pointset which inherits from class Points (Points.py is found in auto/07p/python).
Let's see what members the class Pointset defines.
In [ ]:
pb_obj.branches[0].name
In [ ]:
pb_obj.branches[0].keys()
In [ ]:
pb_obj.branches[0]['tot']
In [ ]:
pb_obj.branches[0].todict()
The object pb_obj.branches[0] holds the coordinates for the first branch in b.lpa in an accessible format -- perfect for plotting.
All we need now are stability properties and bifurcation points on this branch!
In [ ]:
pb_obj.branches[0].labels
This branch is not very eventful: All we get are one endpoint at index 0 and another endpoint at index 46.
Let's check out the next branch.
In [ ]:
pb_obj.branches[1].labels
In [ ]:
pb_obj.branches[1]['tot'][37]
This branch point BP (in this case, a transcritical bifurcation point) is one of four bifurcation points described in our publication -- 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.
Hence, in this publcation we reported T = 23.0 as bifurcation point whereas here we observe tot = 2.3 -- these bifurcation points are equivalent.
We now have the coordinates of all branches in b.lpa and special points (bifurcation points) along these branches.
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.
A Python module parseD.py exists in auto/07p/python so let's see what this offers.
In [ ]:
import parseD as pd
In [ ]:
pd_obj = pd.parseD()
pd_obj.read(open('d.lpa', 'r'))
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!).
In [ ]:
pd_obj[0]
In [ ]:
for i in range(60,100):
print 'Branch number',pd_obj[i]['Branch number'],'eigenvalues',pd_obj[i]['Eigenvalues']
In [ ]:
pd_obj[71]
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.
In [ ]:
import parseS as ps
In [ ]:
ps_obj = ps.parseS()
In [ ]:
ps_obj.read(open('s.lpa', 'r'))
In [ ]:
ps_obj[0]
In [ ]:
ps_obj[0].__str__()
In [ ]:
ps_obj[0].data
The above bifurcation plot is close to the one shown in Figure 3 of a previous article where the same analysis was done on the same system.
What's still missing are stability data (eigenvalues along the solution branches) and bifurcation points.
The file d.lpa generated by AUTO seems to contain that information.
In [ ]:
dlpa = None
with open('d.lpa', 'r') as f:
dlpa = f.readlines()
In [ ]:
dlpa[:20]