Since Sage Days 20 at Marseille held in January 2010, I have been doing the same example over and over again each time I showed someone else how object oriented coding works in Python: using fruits, strawberry, oranges and bananas.

Here is my file: fruit.py. I use to build it from scratch by adding one line at a time using attach command to see what has changed starting with Banana class, then Strawberry class then Fruit class which gathers all common methods.

This time, I wrote the complete documentation (all tests pass, coverage is 100%) and followed Sage coding convention as far as I know them. Thus, I hope this file can be useful as an example to explain those coding convention to newcomers.

One may check that all tests pass using:

$ sage -t fruit.py
sage -t  "fruit.py"
[3.7 s]

----------------------------------------------------------------------
All tests passed!
Total time for all tests: 3.8 seconds

One may check that documentation and doctest coverage is 100%:

$ sage -coverage fruit.py
----------------------------------------------------------------------
fruit.py
SCORE fruit.py: 100% (10 of 10)
----------------------------------------------------------------------

This message is about differences between a Python function and a symbolic function. This is also explained in the Some Common Issues with Functions page in the Sage Tutorial.

In Sage, one may define a symbolic function like:

sage: f(x) = x^2-1

And draw it using one the following way (both works):

sage: plot(f, (x,-10,10))
sage: plot(f(x), (x,-10,10))

Here both f and f(x) are symbolic expressions:

sage: type(f)
<type 'sage.symbolic.expression.Expression'>
sage: type(f(x))
<type 'sage.symbolic.expression.Expression'>

although there are different:

sage: f
x |--> x^2 - 1
sage: f(x)
x^2 - 1

Now if f is a Python function defined with a def statement:

sage: def f(x):
....:     if x>0:
....:         return x
....:     else:
....:         return 0

It is really a Python function:

sage: f
<function f at 0xb933470>
sage: type(f)
<type 'function'>

As above, one can draw the function f:

sage: plot(f, (x,-10,10))

But be carefull, drawing f(x) will not work as expected:

sage: plot(f(x), (x,-10,10))

Why? Because, the python function f gets evaluated on the variable x and this may either raise an exception depending on the definition of f or return some result which might not be a symbolic expression. Here f(x) gets always evaluated to zero because in the definition of f, bool(x > 0) returns False:

sage: x
x
sage: bool(x > 0)
False
sage: f(x)
0

Hence the following constant function is drawn:

sage: plot(0, (x,-10,10))

which is not what we want.

Today, I am presenting the Chapter 3 of the book Probability on Graphs of Geoffrey Grimmett during a monthly reading seminar at LIAFA. The title of the chapter is Percolation and self-avoiding walks. I did some computations to improve my intuition on the question. My code is in the following file : bond_percolation.sage. This post is about some of my computations. You might want to test them yourself online using the Sage Cell Server.

Some samples when p=0.5

A bond percolation sample inside the box \(\Lambda(m)=[-m,m]^d\) when \(p=0.5\) and \(d=2\):

sage: S = BondPercolationSample(p=0.5, d=2)
sage: S.plot(m=40, pointsize=10, thickness=1)
Graphics object consisting of 7993 graphics primitives
sage: _.show()

Another time gives something different:

sage: S = BondPercolationSample(p=0.5, d=2)
sage: S.plot(m=40, pointsize=10, thickness=1)
Graphics object consisting of 10176 graphics primitives
sage: _.show()

Some samples for ranges of values of p

From p=0.1 to p=0.9:

sage: percolation_graphics_array(srange(0.1,1,0.1), d=2, m=5)

From p=0.41 to p=0.49:

sage: percolation_graphics_array(srange(0.41,0.50,0.01), d=2, m=5)

From p=0.51 to p=0.59:

sage: percolation_graphics_array(srange(0.51,0.60,0.01), d=2, m=5)

Upper bound and lower bound for percolation probability \(\Theta(p)\)

In every case, we have the following upper bound for the percolation probability: \[ \Theta(p) = \mathbb{P}_p(\vert C_0\vert=\infty) \leq \mathbb{P}_p(\vert C_0\vert > 1) = 1 - \mathbb{P}_p(\vert C_0\vert = 1) = 1 - (1-p)^{2d}. \] In particular, if \(p\neq 1\), then \(\Theta(p)<1\). In Sage, define the upper bound:

sage: p,n = var('p,n')
sage: d = var('d')
sage: upper_bound = 1 - (1-p)^(2*d)

Also, from Equation (3.8), we have the following lower bound: \[ \Theta(p) \geq 1 - \sum_{n=4}^{\infty} n (4(1-p))^n. \]

In Sage, define the lower bound:

sage: p,n = var('p,n')
sage: lower_bound = 1 - sum(n*(4*(1-p))^n,n,4,oo)
sage: lower_bound.factor()
-(3072*p^5 - 14336*p^4 + 26624*p^3 - 24592*p^2 + 11288*p - 2057)/(4*p - 3)^2

This is not defined when \(p=3/4\), but we are interested in the values in the interval \(]3/4,1]\). In particular, for which value of \(p\) is this lower bound strictly larger than zero:

sage: root = lower_bound.find_root(0.76, 0.99); root
0.8639366490304586

Let's now draw a graph of the lower and upper bound:

sage: U = plot(upper_bound(d=2),(0,1),color='red', thickness=3)
sage: L = plot(lower_bound,(0.86,1),color='green', thickness=3)
sage: G = U + L
sage: G += point((root, 0), color='red', size=20)
sage: lower = r"$1-\sum_{n=4}^{\infty} n4^n(1-p)^n$"
sage: upper = r"$1 -(1-p)^{4}$"
sage: title = r"$1-\sum_{n=4}^{\infty} n4^n(1-p)^n\leq\Theta(p)\leq 1 -(1-p)^{2d}$"
sage: G += text(title, (.5, 1.05), color='black', fontsize=15)
sage: G += text(upper, (.3, 0.5), color='red', fontsize=20)
sage: G += text(lower, (.7, 0.5), color='green', fontsize=20)
sage: G += text("%.5f"%root,(0.88, .03), color='green', horizontal_alignment='left')
sage: G.show()

Thus we conclude that \(\Theta(p) >0\) for \(p>0.8639\) and thus \(p_c \leq 0.8639\).

Percolation probability - dimension 2

The code allows to define the percolation probability function for a given dimension d. It generates n samples and consider the cluster to be infinite if its cardinality is larger than the given stop value.

Here we use Sage adaptative recursion algorithm for drawing the plot of the percolation probability which finds the particular important intervals to ask for more values of the function. See help section of plot function for details. Because T might be long to compute we start with only 4 points.

When stop=100:

sage: T = PercolationProbability(d=2, n=10, stop=100)
sage: T.return_plot((0,1),adaptive_recursion=4,plot_points=4).show()

When stop=1000:

sage: T = PercolationProbability(d=2, n=10, stop=1000)
sage: T.return_plot((0,1),adaptive_recursion=4,plot_points=4).show()

When stop=2000:

sage: T = PercolationProbability(d=2, n=10, stop=2000)
sage: T.return_plot((0,1),adaptive_recursion=4,plot_points=4).show()

Percolation probability - dimension 3

When stop=100:

sage: T = PercolationProbability(d=3, n=10, stop=100)
sage: T.return_plot((0,1),adaptive_recursion=4,plot_points=4).show()

When stop=1000:

sage: T = PercolationProbability(d=3, n=10, stop=1000)
sage: T.return_plot((0,1),adaptive_recursion=4,plot_points=4).show()

Percolation probability - dimension 4

When stop=100:

sage: T = PercolationProbability(d=4, n=10, stop=100)
sage: T.return_plot((0,1),adaptive_recursion=4,plot_points=4).show()

Percolation probability - dimension 13

When stop=100:

sage: T = PercolationProbability(d=13, n=10, stop=100)
sage: T.return_plot((0,1),adaptive_recursion=4,plot_points=4).show()

sage: L = [1, 4, 12, 36, 100, 284, 780, 2172, 5916, 16268, 44100, 120292,
324932, 881500, 2374444, 6416596, 17245332, 46466676, 124658732, 335116620,
897697164, 2408806028, 6444560484, 17266613812, 46146397316, 123481354908,
329712786220, 881317491628]
sage: for k in range(1, len(L)): print numerical_approx(L[k]^(1/k))
4.00000000000000
3.46410161513775
3.30192724889463
3.16227766016838
3.09502148400370
3.03400133198980
2.99705187539871
2.96144397263395
2.93714926770637
2.91369345857619
2.89627439045790
2.87949308754677
2.86632078916860
2.85362749495679
2.84328447096562
2.83329615650289
2.82493415671599
2.81684125361654
2.80992368218258
2.80321554383456
2.79738645741910
2.79172363211806
2.78673687369245
2.78188437392354
2.77756387722633
2.77335345579129
2.76956977331575

sage: L = [1, 6, 30, 150, 726, 3534, 16926, 81390, 387966, 1853886,
8809878, 41934150, 198842742, 943974510, 4468911678, 21175146054,
100121875974, 473730252102, 2237723684094, 10576033219614, 49917327838734,
235710090502158, 1111781983442406, 5245988215191414, 24730180885580790,
116618841700433358, 549493796867100942,2589874864863200574,
12198184788179866902, 57466913094951837030, 270569905525454674614]
sage: for k in range(1, len(L)): print numerical_approx(L[k]^(1/k))
6.00000000000000
5.47722557505166
5.31329284591305
5.19079831727404
5.12452137580198
5.06709510955294
5.02933019629493
4.99573287588832
4.97111339009676
4.94876680377358
4.93129192790635
4.91521453865211
4.90209314463520
4.88990167518413
4.87964724632057
4.87004597517131
4.86178722582108
4.85400655861169
4.84719703702142
4.84074902256992
4.83502763526502
4.82958688248615
4.82470487210973
4.82004549244633
4.81582557693112
4.81178552451599
4.80809774735294
4.80455755518719
4.80130435575213
4.79817388859565

sage: k < 2.76956977331575
k < 2.76956977331575
sage: _ / (2.76956977331575  * k)
0.361066909970928 < (1/k)
sage: 1 - 0.361066909970928
0.638933090029072

sage: k < 4.79817388859565
k < 4.79817388859565
sage: _ / (4.79817388859565 * k)
0.208412621805310 < (1/k)

Below are some Sage tricks that I gathered from other users of Sage, from sage-devel and other places since one year.

Stop the focus in the Notebook

This Notebook hack of the day was published on sage-devel by William Stein during May 2012 to fix that focus movement in the Notebook:

html('<script>cell_focus=function(x,y){} </script>')

Consult the documentation of a function in the browser

Open the documentation of a particular function in a web browser, from either the command-line or the notebook:

sage: browse_sage_doc(factor)

It works also for methods of an object:

sage: m = matrix(2, range(4))
sage: browse_sage_doc(m.inverse)

I found this command when I consulted the help():

sage: help()

Implicit multiplication

Some behavior in Sage can be made implicit like multiplication and variable definition. This might be good for new users coming from Maple for instance.

Normally, this syntax raises an error:

sage: 34x
------------------------------------------------------------
   File "<ipython console>", line 1
     34x
       ^
SyntaxError: invalid syntax

It is possible to make it work:

sage: implicit_multiplication(True)
sage: 34x
34*x

Implicit variable definition

The following works only in the Sage Notebook. It allows to turn on automatic definition of variables:

sage: automatic_names(True)
sage: x + y + z
x + y + z

Turn out automatic show when using plot

Set the default to False for showing plots using any plot commands:

sage: show_default(False)

I prefer False but you may not.

Rerun the patchbot

One can re-run the tests of a particular ticket by adding ?kick to the url of the ticket on the patchbot Sage server. For example, to rerun the tests on the ticket #13461, one can load the following page:

http://patchbot.sagemath.org/ticket/13461/?kick

This trick was shared on sage-devel during August 2012.

Python debugger

The majority of the Best Practices for Scientific Computing are followed by the Sage Development model. But, there is one principle that at least I do not use enough : a debugger. However, a Python debugger exists and a tutorial for using the Python debugger is available on onlamp.com.

I remember that discussion on sage-devel Poll: which debugger do you use? where developers were sharing their debugger tricks. It seems that using print statements is unavoidable.

Computing basic statistics with R

I always knew R is in Sage but never used it even if sometimes I need to compute some statistics of list of numbers. In Sage, one can print statistics of a list of numbers using R:

sage: L = [randint(0,100) for _ in range(20)]
sage: r.summary(L)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  14.00   27.75   57.50   53.45   69.00   98.00

Creation of the number field in \(\sqrt{5}\)

Before learning it from this video of William Stein, I did not know that the square brackets could be used to create such number field:

sage: A = ZZ[sqrt(5)]
sage: a,sqrt5 = A.gens()
sage: A
Order in Number Field in sqrt5 with defining polynomial x^2 - 5
sage: sqrt5
sqrt5

Using Sage locally in the notebook from a server

First, log into the server using the following port setup:

ssh -L 8389:localhoat:8389 [SERVER]

Start the notebook with the given port:

sage -notebook port=8389

This ask for a password, if you forget it, you may reset it by first opening Sage, and by starting the notebook with the option reset=True:

sage: notebook(port=8389, reset=True)

Then, by opening the browser at the following adress, I can log in to the notebook from the server:

http://localhost:8389/home/admin/

Going from Mercurial to Git

Sage will soon move from Mercurial to Git. In the past, I tried to understand the difference between Mercurial and Git. I was never able to find a simple text or blog post on the web explaining in a simple way the differences. One diffence would be about branching. But I was not using branches with Mercurial... As I see it, differences between Mercurial and Git are hard to explain or at least hard to understand.

Anyway, I once found this tutorial on the Git version control system: Understanding Git Conceptually where the approach is "conceptual" and maybe easier for the mathematician. For instance, in Section 1 it is written that repositories can be "visualized as a directed acyclic graph of commit objects". I still haven't go through this tutorial but I consider to start with this one.

Inheritance tree

One may draw the inheritance tree of a class with the following command:

sage: class_graph(Integer).plot()

La toute dernière version de Blogofile 0.8.b vient de sortir plus d'un an après la version 0.7.1 que j'utilise pour ce site web. Blogofile est un générateur de site web statique, c'est-à-dire qui dépend uniquement de fichiers textes. Il existe une multitude de projets similaires dont Pelican, nanoc, Octopress, Cytoplasm, rstblog et plusieurs autres. Le site web mathematism.com a fait le choix de nanoc et mickgardner.com a choisi Jekyll. Un autre a choisi de coder son propre ReStructuredText-Writer en le basant sur docutils.

Pour ma part, j'ai choisi Blogofile, car il est basé sur le langage Python, comprend la syntaxe ReStructuredText et aussi parce que j'ai été capable de l'utiliser dès la première journée contrairement au projet Hyde aussi basé sur Python, mais que je n'ai jamais réussi à comprendre.

D'autres personnes ont fait le même choix que moi tel que le créateur de Mako Templates qui explique son choix dans son texte how coders blog ainsi que ses blogofile hacks. Finalement, voici trois autres sites qui expliquent leur passage à blogofile: asktherelic.com, paolocorti.net et lukeplant.me.uk.