\documentclass[12pt]{article}

%% To process:
%% htlatex coast.tex "/home/rob/conversion/tex2sws/tex4ht-sage.cfg" " -cunihtf -utf8"
%% ~/conversion/tex2sws/tex2sws.py

\usepackage{sagetex}
\usepackage{graphicx}
\setlength{\parindent}{0pt}

\title{Spicing Up Your Research in Combinatorics with Sage}
\date{February 12, 2011}

\begin{document}
\maketitle
\begin{center}
\Large
Rob Beezer\\
University of Puget Sound\\
\ \newline
14th Coast Combinatorics Conference\\
University of Victoria\\
\end{center}
%
\begin{sageverbatim}

\end{sageverbatim}
%
\section{Disclaimers}
%
$\bullet$ I may be the wrong person to give this talk!\\
$\bullet$ Open source fanatic\\
$\bullet$ Member, Sage developer community\\
$\bullet$ Sage development: linear algebra, graph theory, group theory\\
$\bullet$ But I will \emph{try} to be objective
%
\section{Introduction}
% List OK, section ends right away
\begin{itemize}
\item  Sage: Open source software for mathematics
\item  Mission: ``Creating a viable free open source alternative to Magma, Maple, Mathematica and Matlab.''
\item  Consolidates over 100 packages, such as:
    \begin{itemize}
    \item NetworkX --- Graph Theory
    \item GAP --- Groups, Algorithms, Programming
    \item ATLAS --- Automatically Tuned Linear Algebra System
    \item IML --- Integer Matrix Library
    \end{itemize}
\item Command Line, Batch, Notebook Interface
\item FREE!
\end{itemize}

\section{Graph Theory}

Create graphs in a natural way:
%
\begin{sageverbatim}
harary = Graph([(0,1), (1,2), (2,3), (3,0), (1,3)])
\end{sageverbatim}
%
\begin{sageverbatim}
harary
\end{sageverbatim}
%
\begin{sageverbatim}
harary.plot()
\end{sageverbatim}
%
\begin{sageverbatim}
harary.num_verts(), harary.num_edges()
\end{sageverbatim}
%
\begin{sageverbatim}
harary.is_planar()
\end{sageverbatim}
%
\begin{sageverbatim}
H=harary.hamiltonian_cycle()
H.plot()
\end{sageverbatim}
%
\begin{sageverbatim}
harary.degree_sequence()
\end{sageverbatim}
%
\begin{sageverbatim}
sorted(harary.degree_sequence())
\end{sageverbatim}
%
\subsection{Graph Generators}
%
Many pre-defined graphs (digraphs, too):
%
\begin{sageverbatim}
graphs.
\end{sageverbatim}
%
\begin{sageverbatim}
trees_iterator = graphs.trees(8)
T8 = list(trees_iterator)
T8
\end{sageverbatim}
%
From a path to a star:
%
\begin{sageverbatim}
[tree.diameter() for tree in T8]
\end{sageverbatim}
%
Visually:
%
\begin{sageverbatim}
graphs_list.show_graphs(T8)
\end{sageverbatim}
%
\subsection{Tower of Hanoi Graph}
%
\begin{center}
\includegraphics{4peghanoi.jpg}\par
\end{center}
%
$\bullet$  \verb!graphs.HanoiTowerGraph(n, d)!\\
$\bullet$  Generalize to $n$ pegs and $d$ disks\\
$\bullet$  State graph: intermediate configurations, edges are ``one move''\\
$\bullet$  Labels: $d$-tuple, large disk to small disk; entries are peg numbers\\
$\bullet$  Example: $n=3$, $d=4$: $(2,0,2,1)$\\
\par
%
\begin{sageverbatim}
T = graphs.HanoiTowerGraph(3, 4, positions=True)
T.show(figsize=12)
\end{sageverbatim}
%
A solution is path between states ``all the disks on one peg'' and ``all the disks on another peg.''
%
\begin{sageverbatim}
solution=T.shortest_path((0,0,0,0), (1,1,1,1))
solution
\end{sageverbatim}
%
Minimum number of moves:
%
\begin{sageverbatim}
len(solution) - 1
\end{sageverbatim}
%
\begin{sageverbatim}
T.diameter()
\end{sageverbatim}
%
More general:
%
\begin{sageverbatim}
T = graphs.HanoiTowerGraph(4, 3, positions=True)
T.show(figsize=12)
\end{sageverbatim}
%
%
\begin{sageverbatim}
T = graphs.HanoiTowerGraph(4, 4, labels=False, positions=True)
T.show(figsize=12)
\end{sageverbatim}
%
Forget about graphics, work with graph itself.
%
\begin{sageverbatim}
T = graphs.HanoiTowerGraph(4, 8, labels=False, positions=False)
T.num_verts()
\end{sageverbatim}
%
Code vertices to integers: $d$-tuples, base $n$.  All disks on peg 0, move to all disks on peg 3.
\begin{sageverbatim}
solution = T.shortest_path(0, 4^8-1)
solution
\end{sageverbatim}
%
\begin{sageverbatim}
len(solution)-1
\end{sageverbatim}
%
Theorem: automorphisms of the state graph are just the obvious ones (renumber pegs)
%
\begin{sageverbatim}
T = graphs.HanoiTowerGraph(4, 6, labels=False, positions=False)
A = T.automorphism_group()
A.order()
\end{sageverbatim}
%
\begin{sageverbatim}
S4 = SymmetricGroup(4)
S4.is_isomorphic(A)
\end{sageverbatim}
%
Automorphisms are computed via Brendan McKay's \verb!nauty! algorithm, re-implemented as \verb!NICE!.
%
\subsection{Online Help}
%
Tab-completion, online syntax and tested examples (?), source code (??).
%
\begin{sageverbatim}
G = graphs.GrotzschGraph()
\end{sageverbatim}
%
\begin{sageverbatim}
G.degree_seq
\end{sageverbatim}
%
\subsection{Algebraic Graph Theory}
%
Comprehensive support for groups via GAP.  If need be, can ship \emph{any} GAP command directly from Sage to GAP.
%
\begin{sageverbatim}
P = graphs.PetersenGraph()
A = P.automorphism_group()
A
\end{sageverbatim}
%
\begin{sageverbatim}
A.order()
\end{sageverbatim}
%
\begin{sageverbatim}
A.is_transitive()
\end{sageverbatim}
%
\begin{sageverbatim}
S = A.stabilizer(9); S
\end{sageverbatim}
%
\begin{sageverbatim}
sg = S.subgroups(); sg
\end{sageverbatim}
%
\begin{sageverbatim}
[H.order() for H in sg]
\end{sageverbatim}
%
\begin{sageverbatim}
H.is_isomorphic(DihedralGroup(6))
\end{sageverbatim}
%%
\begin{sageverbatim}
ds = A.derived_series()
\end{sageverbatim}
%
\begin{sageverbatim}
[K.order() for K in ds]
\end{sageverbatim}
%
\begin{sageverbatim}
B = ds[1]
B.is_isomorphic(AlternatingGroup(5))
\end{sageverbatim}
%
Comprehensive support for linear algebra.
%
\begin{sageverbatim}
A = random_matrix(ZZ, 1000, 1000)
time A.determinant()
\end{sageverbatim}
%
\begin{sageverbatim}
K = graphs.KneserGraph(8,3)
\end{sageverbatim}
%
\begin{sageverbatim}
adj = K.adjacency_matrix()
adj
\end{sageverbatim}
%
\begin{sageverbatim}
K.spectrum()
\end{sageverbatim}
%
A small ``singular graph.''  (I.\ Sciriha, 2007)
%
\begin{sageverbatim}
S = graphs.CycleGraph(4)
S.add_vertices([4, 5, 6])
S.add_edges([(2,4), (2,5), (2,6)])
S.add_edges([(3,4), (3,5), (3,6)])
S.plot()
\end{sageverbatim}
%
\begin{sageverbatim}
graph_editor(S)
\end{sageverbatim}
%
\begin{sageverbatim}
adj = S.adjacency_matrix()
ker = adj.kernel()
ker
\end{sageverbatim}
%
Notice this is the kernel over the integers, and is computed as a module.  It is easy to upgrade to the rationals.
%
\begin{sageverbatim}
adjQ = adj.change_ring(QQ)
kerQ = adjQ.kernel()
kerQ
\end{sageverbatim}
%
A vector space in Sage behaves like a vector space, it has dimension, you can interect two of them, \dots
%
%
\begin{sageverbatim}
kerQ.
\end{sageverbatim}
%
\section{Combinatorics}
%
MuPAD-Combinat was sold to MathWorks (i.e.\ Matlab) in September 2008.  With foresight, MuPAD contributors had switched to Sage in June 2008.  So there is an active and experienced group contributing combinatorics to Sage as ``sage-combinat.''
\par
%
The usual functions:
%
\begin{sageverbatim}
stirling_number2(500,30)
\end{sageverbatim}
%
\begin{sageverbatim}
number_of_derangements(["a", 'b', 'c', 'd', 'e'])
\end{sageverbatim}
%
The usual objects (note multisets):
%
\begin{sageverbatim}
permutations([0,1,1,2,2,2])
\end{sageverbatim}
%
\begin{sageverbatim}
combinations([0,1,1,2,2,2], 3)
\end{sageverbatim}
%
\begin{sageverbatim}
derangements(["a", 'b', 'c', 'd', 'e'])
\end{sageverbatim}
%
Ranking and Unranking (uses classes):\\
(preserve order for \verb!.rank()!)
%
%
\begin{sageverbatim}
C = Combinations(["cat", "dog", "pig", "fox", "bug", "rat"], 3)
\end{sageverbatim}
%
\begin{sageverbatim}
C.unrank(15)
\end{sageverbatim}
%
\begin{sageverbatim}
C.rank(['cat', 'dog', 'bug'])
\end{sageverbatim}
%
Polynomials:\\
%
A symbolic object:
%
\begin{sageverbatim}
var('t')
p = bernoulli_polynomial(t, 5); p
\end{sageverbatim}
%
\begin{sageverbatim}
diff(p, t)
\end{sageverbatim}
%
\begin{sageverbatim}
p.parent()
\end{sageverbatim}
%
A function:
%
\begin{sageverbatim}
var('t')
q(t) = bernoulli_polynomial(t, 5); q
\end{sageverbatim}
%
\begin{sageverbatim}
q(10)
\end{sageverbatim}
%
\begin{sageverbatim}
q.parent()
\end{sageverbatim}
%
A polynomial over the rationals, i.e.\ in $\mathbf Q[t]$:
%
\begin{sageverbatim}
t = polygen(QQ, 't')
r = bernoulli_polynomial(t, 5); r
\end{sageverbatim}
%
\begin{sageverbatim}
r.is_irreducible()
\end{sageverbatim}
%
\begin{sageverbatim}
r.factor()
\end{sageverbatim}
%
\begin{sageverbatim}
r.parent()
\end{sageverbatim}
%
Lots of esoterica:
%
\begin{sageverbatim}
S = SkewPartitions(4)
s = S[25]; s
\end{sageverbatim}
%
\begin{sageverbatim}
print s.diagram()
\end{sageverbatim}
%
Much, much more: crystals, posets, root systems, symmetric functions, words.
%
\section{Goodies}
%
Excellent \LaTeX support.
%
\begin{sageverbatim}
latex(integrate(sec(x), x))
\end{sageverbatim}
%
\begin{sageverbatim}
A = random_matrix(QQ, 4, num_bound=9, den_bound=9)
latex(A)
\end{sageverbatim}
%
\LaTeX\ versions of graphs:
%
\begin{sageverbatim}
P = graphs.PetersenGraph()
P.set_latex_options(vertex_shape='diamond', vertex_color='red', vertex_label_color='gold', edge_color='blue')
\end{sageverbatim}
%
\begin{sageverbatim}
latex(P)
\end{sageverbatim}
%
\begin{sageverbatim}
from sage.graphs.graph_latex import setup_latex_preamble
setup_latex_preamble()
latex.jsmath_avoid_list('tikzpicture')
\end{sageverbatim}
%
\begin{center}
\includegraphics{petersen-uvic.png}
\end{center}
%
Embedded Word Processor (Shift-Click) is \TeX-aware.
%
\begin{sageverbatim}

\end{sageverbatim}
%
%
Sage\TeX\ \ allows you to call Sage to compute values, expressions or plots for automatic inclusion in your \LaTeX\  document.
%

$\bullet$ Input:  \verb!Today's date \$2011-02-12=20110212\$ factors as \sage{factor(20110212)}.!\\
$\bullet$ Output: Today's date $2011-02-12=20110212$ factors as $2^{2} \cdot 3^{2} \cdot 173 \cdot 3229$.
\par
%
Can author in \LaTeX\ and produce worksheets (this is an example).\\
%
\section{Interacts and Python Extensions}
%
$\bullet$ Full support for any Python package.\\
$\bullet$ Easy-to-construct interactive demonstrations.
\par
%
Code adapted from a demonstration by William Stein
%
\begin{sageverbatim}
import pylab
A_image = pylab.mean(pylab.imread(DATA + 'mystery.png'), 2)
@interact
def svd_image(i=(1,(1..50)), display_axes=True):
    u,s,v = pylab.linalg.svd(A_image)
    A     = sum(s[j]*pylab.outer(u[0:,j], v[j,0:]) for j in range(i))
    show(matrix_plot(A), axes=display_axes, figsize=(11,5))
    html('<h2>Compressed using %s singular values</h2>'%i)
\end{sageverbatim}
%
Code from Sage Wiki (Interacts) by Pablo Angulo.
%
\begin{sageverbatim}
%hide
def animate_contraction(g, e, frames = 12, **kwds):
    v1, v2 = e
    if not g.has_edge(v1,v2):
        raise ValueError, "Given edge not found on Graph"
    ls = []
    posd = g.get_pos()
    for j in range(frames):
        gp = Graph(g)
        posdp = dict(posd)
        p1 = posdp[v1]
        p2 = posdp[v2]
        posdp[v2] = [a*(frames-j)/frames + b*j/frames
                    for a,b in zip(p2,p1)]

        gp.set_pos(posdp)
        ls.append(plot(gp, **kwds))
    return ls

def animate_vertex_deletion(g, v, frames = 12, **kwds):
    kwds2 = dict(kwds)
    if 'vertex_colors' in kwds:
        cs = dict(kwds['vertex_colors'])
        for c, vs in kwds['vertex_colors'].items():
            if v in vs:
                vs2 = list(vs)
                vs2.remove(v)
                cs[c] = vs2
        kwds2['vertex_colors'] = cs
    else:
        kwds2 = dict(kwds)
    g2 = Graph(g)
    posd = dict(g.get_pos())
    del posd[v]
    g2.delete_vertex(v)
    g2.set_pos(posd)
    return [plot(g, **kwds),plot(g2, **kwds2)]*int(frames/2)

def animate_edge_deletion(g, e, frames = 12, **kwds):
    v1, v2 = e
    g2 = Graph(g)
    g2.delete_edge(e)
    return [plot(g, **kwds),plot(g2, **kwds)]*int(frames/2)

def animate_glide(g, pos1, pos2, frames = 12, **kwds):
    ls = []
    for j in range(frames):
        gp = Graph(g)
        pos = {}
        for v in gp.vertices():
            p1 = pos1[v]
            p2 = pos2[v]
            pos[v] = [b*j/frames + a*(frames-j)/frames
                        for a,b in zip(p1,p2)]
        gp.set_pos(pos)
        ls.append(plot(gp, **kwds))
    return ls

def medio(p1, p2):
    return tuple((a+b)/2 for a,b in zip(p1, p2))

def new_color():
    return (0.1+0.8*random(), 0.1+0.8*random(), 0.1+0.8*random())

def animate_minor(g, m, frames = 12, pause = 50, step_time = 100):
    '''Crea una animación que muestra cómo un grafo tiene un menor m
    '''
    posd = dict(g.get_pos())
    posg = posd.values()
    posm = m.get_pos().values()
    xmax = max(max(x for x,y in posg), max(x for x,y in posm))
    ymax = max(max(y for x,y in posg), max(y for x,y in posm))
    xmin = min(min(x for x,y in posg), min(x for x,y in posm))
    ymin = min(min(y for x,y in posg), min(y for x,y in posm))
    dd = g.minor(m)

    #Set colors
    m_colors = dict((v,new_color()) for v in m.vertices())
    g_colors = dict((m_colors[k],vs)
                        for k,vs in dd.items())

    extra_vs = (set(g.vertices()) -
                set(v for vs in dd.values()
                      for v in vs))
    g_colors[(0,0,0)] = list(extra_vs)

    #pics contains the frames of the animation
    #no colors at the beggining
    gg = Graph(g)
    pics = [plot(gg)]*frames

    #First: eliminate extra vertices
    for v in extra_vs:
        pics.extend(animate_vertex_deletion(gg, v, frames,
                            vertex_colors = g_colors))
        gg.delete_vertex(v)
        del posd[v]
        g_colors[(0,0,0)].remove(v)
    del g_colors[(0,0,0)]

    #Second: contract edges
    for color, vs in g_colors.items():
        while len(vs)>1:
            for j in xrange(1,len(vs)):
                if gg.has_edge(vs[0], vs[j]):
                    break
            pics.extend(animate_contraction(gg, (vs[0], vs[j]), frames,
                                vertex_colors = g_colors))
            for v in gg.neighbors(vs[j]):
                gg.add_edge(vs[0],v)
            gg.delete_vertex(vs[j])
            del posd[vs[j]]
            gg.set_pos(posd)
            posd = dict(gg.get_pos())
            del vs[j]

    #Relabel vertices of m so that they match with those of gg
    m = Graph(m)
    dd0 = dict((k, vs[0])
                  for k,vs in dd.items() )
    m.relabel(dd0)


    #Third: glide to position in m
    pos_m = m.get_pos()
    pos_g = gg.get_pos()
    pics.extend(animate_glide(gg, pos_g, pos_m, frames,
                                vertex_colors = g_colors))
    gg.set_pos(pos_m)

    #Fourth: delete redundant edges
    for e in gg.edges(labels = False):
        if not m.has_edge(e):
            pics.extend(animate_edge_deletion(gg, e, frames,
                                vertex_colors = g_colors))
            gg.delete_edge(*e)

    #And wait for a moment
    pics.extend([plot(gg, vertex_colors = g_colors)]*frames)

    return animate(pics, xmin = xmin - 0.1, xmax = xmax + 0.1,
                         ymin = ymin - 0.1, ymax = ymax + 0.1)
\end{sageverbatim}
%
\begin{sageverbatim}
graph_list = {'CompleteGraph4':graphs.CompleteGraph(4),
              'CompleteGraph5':graphs.CompleteGraph(5),
              'CompleteGraph6':graphs.CompleteGraph(6),
              'BipartiteGraph3,3':graphs.CompleteBipartiteGraph(3,3),
              'BipartiteGraph4,3':graphs.CompleteBipartiteGraph(4,3),
              'PetersenGraph':graphs.PetersenGraph(),
              'CycleGraph4':graphs.CycleGraph(4),
              'CycleGraph5':graphs.CycleGraph(5),
              'BarbellGraph(3,2)':graphs.BarbellGraph(3,2)
              }
@interact
def _(u1 = text_control(value='Does this graph'),
      g  = selector(graph_list.keys(), buttons = True),
      u2 = text_control(value='have a minor isomorphic to this other graph:'),
      m  = selector(graph_list.keys(), buttons = True),
      u3 = text_control(value='''? It is has, show it to me, 
      with an animation with the following parameters'''),
      frames = slider(1,15,1,4, label = 'Frames per second'),
      step_time = slider(1/10,2,1/10,1, label = 'Seconds per step')):
    g = graph_list[g]
    m = graph_list[m]
    if g == m:
        html('<h3>Please select two different graphs</h3>')
        return
    try:
        a = animate_minor(g, m, frames = frames)
        a.show(delay=100*step_time/frames)
    except ValueError:
        html('''<h3>The first graph have <strong>NO</strong> minor isomorphic to the second</h3>''')
\end{sageverbatim}
%
\section{Miscellaneous}
%
% List OK, section ends
%
\begin{itemize}
\item Cython: a Python-to-C translator, creates core routines in Sage
\item Optimization: CVXOPT, Exact cover with DLX, etc.
\item Active community:  $\sim$200 developers, mathematicians and programmers.
\item Online servers
\item Interfaces to GAP, PARI, Maxima, R, Singular (Python, LaTeX) --- all included packages
\item Interfaces to Mathematica, Maple, Matlab, Magma, Axiom, Kash, Macaulay2, MuPad, Octave (not included)
\item More packages: SymPy, NumPy, Linbox, MPFR, MPFI, NTL, eclib, ATLAS, FLINT, lcalc, PolyBori PyCrypto, matplotlib
\end{itemize}
%
\section{Philosophy}
%
Mathematica: ``Why You Do Not Usually Need to Know about Internals''\\
In the online Mathematica Documentation Center.
%
\begin{quote}
You should realize at the outset that while knowing about the internals of Mathematica may be of intellectual interest, it is usually much less important in practice than you might at first suppose.\par
%
\dots\par
%
Indeed, in almost all practical uses of Mathematica, issues about how Mathematica works inside turn out to be largely irrelevant.\par
%
\dots\par
%
Particularly in more advanced applications of Mathematica, it may sometimes seem worthwhile to try to analyze internal algorithms in order to predict which way of doing a given computation will be the most efficient. And there are indeed occasionally major improvements that you will be able to make in specific computations as a result of such analyses.\par
%
But most often the analyses will not be worthwhile. For the internals of Mathematica are quite complicated, and even given a basic description of the algorithm used for a particular purpose, it is usually extremely difficult to reach a reliable conclusion about how the detailed implementation of this algorithm will actually behave in particular circumstances.
%
\end{quote}
%
J. Neub\"user (Founder of GAP): ``An invitation to computational group theory.''\\
In C. M. Campbell, T. C. Hurley, E. F. Robertson, S. J. Tobin, and J. J. Ward, editors, Groups ‘93 Galway/St. Andrews, Volume 2, volume 212 of London Mathematical Society Lecture Note Series, pages 457–-475. Cambridge University Press, 1995.
%
\begin{quote}
You can read Sylow's Theorem and its proof in Huppert's book in
the library without even buying the book and then you can use
Sylow's Theorem for the rest of your life free of charge,
but\dots\ for many computer algebra systems license fees have to be
paid regularly for the total time of their use. In order to
protect what you pay for, you do not get the source, but only an
executable, i.e.\ a black box. You can press buttons and you get
answers in the same way as you get the bright pictures from your
television set but you cannot control how they were made in either
case.

With this situation two of the most basic rules of conduct in
mathematics are violated: In mathematics information is passed on
free of charge and everything is laid open for checking. Not
applying these rules to computer algebra systems that are made for
mathematical research\dots\ means moving in a most undesirable
direction. Most important: Can we expect somebody to believe a
result of a program that he is not allowed to see? Moreover: Do we
really want to charge colleagues in Moldava several years of their
salary for a computer algebra system?
\end{quote}


\end{document}