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.
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
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.
sys.path to import the Python modules defined by
import sys auto_directory = !echo $AUTO_DIR if auto_directory == ['']: home = !echo $HOME auto_directory = home+'/auto/07p/' sys.path.append(auto_directory+'/python') else: sys.path.append(auto_directory + '/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.
import AUTOCommands as ac import AUTOclui as acl import interactiveBindings as ib import runAUTO as ra
Start AUTO and catch the returned
runner = ra.runAUTO()
AUTOCommands defines multiple methods that return solution objects. Method
load is one of them.
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
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:
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.oin 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.exe
By 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.
autotoolchain seems broken by our use of the IPython Notebook
596 raise AUTOExceptions.AUTORuntimeError("Error running AUTO") 597 598 def test(): AUTORuntimeError: Error running AUTO
Despite the fact that something seems broken, we notice that
auto still ran and created three output files.
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.
import pandas as pd from matplotlib import pylab as pl
content = None with open('b.lpa', 'r') as f: content = f.readlines()
Line 15 (zero-indexed) onwards are the branches.
content_csv = [[el for el in content.split(' ') if len(el) > 0 and el != '\n']] content_csv = 'branch' column_names = content_csv
These are our column names found in
for line in content: dummy = line.split(' ') dummy = [el for el in dummy if len(el) > 0 and el != '\n'] if dummy == '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
df = pd.DataFrame(content_csv, columns=column_names)
Plotting branches is easy now.
scatter(df[df['branch'] == 1].tot,df[df['branch'] == 1].aL) scatter(df[df['branch'] == 2].tot,df[df['branch'] == 2].aL)
The manual says that files
fort.9 are equivalent to
d.* (in that order).
fort.* naming scheme may be older as the latest version of AUTO appears to produce files in the latter naming scheme.
auto/07p/python/ is probably a good start for what we want to do here.
parseB.parseBMixin provides a method called
readFilename so let's hope this allows us to read and parse our
pb_obj = parseB.parseBMixin()
Hmmm ... nope.
pb_obj = parseB.parseB()
b_lpa = open('b.lpa', 'r') ab_obj.read(b_lpa)
pb_obj seems to store an array of the branches found in
Branches are saved as instances of
Points.Pointset which inherits from class
Points.py is found in
Let's see what members the class
pb_obj.branches 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!
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.
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 (
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.
import parseD as pd
pd_obj = pd.parseD() pd_obj.read(open('d.lpa', 'r'))
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!).
for i in range(60,100): print 'Branch number',pd_obj[i]['Branch number'],'eigenvalues',pd_obj[i]['Eigenvalues']
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.
import parseS as ps
ps_obj = ps.parseS()
What's still missing are stability data (eigenvalues along the solution branches) and bifurcation points.
d.lpa generated by
AUTO seems to contain that information.
dlpa = None with open('d.lpa', 'r') as f: dlpa = f.readlines()