Blogue

A time evolution picture of packages built in parallel by Sage

16 décembre 2016 | Catégories: sage | View Comments

Compiling sage takes a while and does a lot of stuff. Each time I am wondering which components takes so much time and which are fast. I wrote a module in my slabbe version 0.3b2 package available on PyPI to figure this out.

This is after compiling 7.5.beta6 after an upgrade from 7.5.beta4:

sage: from slabbe.analyze_sage_build import draw_sage_build
sage: draw_sage_build().pdf()
/Files/2016/sage_build.png

From scratch from a fresh git clone of 7.5.beta6, after running MAKE='make -j4' make ptestlong, I get:

sage: from slabbe.analyze_sage_build import draw_sage_build
sage: draw_sage_build().pdf()
/Files/2016/sage_build_from_scratch.png

The picture does not include the start and ptestlong because there was an error compiling the documentation.

By default, draw_sage_build considers all of the logs files in logs/pkgs but options are available to consider only log files created in a given interval of time. See draw_sage_build? for more info.

Read and Post Comments

unsupported operand parent for *, Matrix over number field, vector over symbolic ring

18 février 2016 | Catégories: sage | View Comments

Yesterday I received this email (in french):

Salut,
avec Thomas on a une question bête:

K.<x>=NumberField(x*x-x-1)

J'aimerais multiplier une matrice avec des coefficients en x par un vecteur
contenant des variables a et b.  Il dit "unsupported operand parent for *,
Matrix over number field, vector over symbolic ring"

Est ce grave ?

Here is my answer. Indeed, in Sage, symbolic variables can't multiply with elements in an Number Field in x:

sage: x = var('x')
sage: K.<x> = NumberField(x*x-x-1)
sage: a = var('a')
sage: a*x
Traceback (most recent call last)
...
TypeError: unsupported operand parent(s) for '*': 'Symbolic Ring' and
'Number Field in x with defining polynomial x^2 - x - 1'

But, we can define a polynomial ring with variables in a,b and coefficients in the NumberField. Then, we are able to multiply a with x:

sage: x = var('x')
sage: K.<x> = NumberField(x*x-x-1)
sage: K
Number Field in x with defining polynomial x^2 - x - 1
sage: R.<a,b> = K['a','b']
sage: R
Multivariate Polynomial Ring in a, b over Number Field in x with
defining polynomial x^2 - x - 1
sage: a*x
(x)*a

With two square brackets, we obtain powers series:

sage: R.<a,b> = K[['a','b']]
sage: R
Multivariate Power Series Ring in a, b over Number Field in x with
defining polynomial x^2 - x - 1
sage: a*x*b
(x)*a*b

It works with matrices:

sage: MS = MatrixSpace(R,2,2)
sage: MS
Full MatrixSpace of 2 by 2 dense matrices over Multivariate Power
Series Ring in a, b over Number Field in x with defining polynomial
x^2 - x - 1
sage: MS([0,a,b,x])
[  0   a]
[  b (x)]
sage: m1 = MS([0,a,b,x])
sage: m2 = MS([0,a+x,b*b+x,x*x])
sage: m1 + m2 * m1
[              (x)*b + a*b       (x + 1) + (x + 1)*a]
[                (x + 2)*b (3*x + 1) + (x)*a + a*b^2]
Read and Post Comments

slabbe-0.2.spkg released

30 novembre 2015 | Catégories: sage, slabbe spkg | View Comments

These is a summary of the functionalities present in slabbe-0.2.spkg optional Sage package. It works on version 6.8 of Sage but will work best with sage-6.10 (it is using the new code for cartesian_product merged the the betas of sage-6.10). It contains 7 new modules:

  • finite_word.py
  • language.py
  • lyapunov.py
  • matrix_cocycle.py
  • mult_cont_frac.pyx
  • ranking_scale.py
  • tikz_picture.py

Cheat Sheets

The best way to have a quick look at what can be computed with the optional Sage package slabbe-0.2.spkg is to look at the 3-dimensional Continued Fraction Algorithms Cheat Sheets available on the arXiv since today. It gathers a handful of informations on different 3-dimensional Continued Fraction Algorithms including well-known and old ones (Poincaré, Brun, Selmer, Fully Subtractive) and new ones (Arnoux-Rauzy-Poincaré, Reverse, Cassaigne).

/Files/2015/arp_cheat_sheet.png

Installation

sage -i http://www.slabbe.org/Sage/slabbe-0.2.spkg    # on sage 6.8
sage -p http://www.slabbe.org/Sage/slabbe-0.2.spkg    # on sage 6.9 or beyond

Examples

Computing the orbit of Brun algorithm on some input in \(\mathbb{R}^3_+\) including dual coordinates:

sage: from slabbe.mult_cont_frac import Brun
sage: algo = Brun()
sage: algo.cone_orbit_list((100, 87, 15), 4)
[(13.0, 87.0, 15.0, 1.0, 2.0, 1.0, 321),
 (13.0, 72.0, 15.0, 1.0, 2.0, 3.0, 132),
 (13.0, 57.0, 15.0, 1.0, 2.0, 5.0, 132),
 (13.0, 42.0, 15.0, 1.0, 2.0, 7.0, 132)]

Computing the invariant measure:

sage: fig = algo.invariant_measure_wireframe_plot(n_iterations=10^6, ndivs=30)
sage: fig.savefig('a.png')
/Files/2015/brun_invm_wireframe_plot.png

Drawing the cylinders:

sage: cocycle = algo.matrix_cocycle()
sage: t = cocycle.tikz_n_cylinders(3, scale=3)
sage: t.png()
/Files/2015/brun_cylinders_3.png

Computing the Lyapunov exponents of the 3-dimensional Brun algorithm:

sage: from slabbe.lyapunov import lyapunov_table
sage: lyapunov_table(algo, n_orbits=30, n_iterations=10^7)
  30 succesfull orbits    min       mean      max       std
+-----------------------+---------+---------+---------+---------+
  $\theta_1$              0.3026    0.3045    0.3051    0.00046
  $\theta_2$              -0.1125   -0.1122   -0.1115   0.00020
  $1-\theta_2/\theta_1$   1.3680    1.3684    1.3689    0.00024

Dealing with tikzpictures

Since I create lots of tikzpictures in my code and also because I was unhappy at how the view command of Sage handles them (a tikzpicture is not a math expression to put inside dollar signs), I decided to create a class for tikzpictures. I think this module could be usefull in Sage so I will propose its inclusion soon.

I am using the standalone document class which allows some configurations like the border:

sage: from slabbe import TikzPicture
sage: g = graphs.PetersenGraph()
sage: s = latex(g)
sage: t = TikzPicture(s, standalone_configs=["border=4mm"], packages=['tkz-graph'])

The repr method does not print all of the string since it is often very long. Though it shows how many lines are not printed:

sage: t
\documentclass[tikz]{standalone}
\standaloneconfig{border=4mm}
\usepackage{tkz-graph}
\begin{document}
\begin{tikzpicture}
%
\useasboundingbox (0,0) rectangle (5.0cm,5.0cm);
%
\definecolor{cv0}{rgb}{0.0,0.0,0.0}
...
... 68 lines not printed (3748 characters in total) ...
...
\Edge[lw=0.1cm,style={color=cv6v8,},](v6)(v8)
\Edge[lw=0.1cm,style={color=cv6v9,},](v6)(v9)
\Edge[lw=0.1cm,style={color=cv7v9,},](v7)(v9)
%
\end{tikzpicture}
\end{document}

There is a method to generates a pdf and another for generating a png. Both opens the file in a viewer by default unless view=False:

sage: pathtofile = t.png(density=60, view=False)
sage: pathtofile = t.pdf()
/Files/2015/petersen_graph.png

Compare this with the output of view(s, tightpage=True) which does not allow to control the border and also creates a second empty page on some operating system (osx, only one page on ubuntu):

sage: view(s, tightpage=True)
/Files/2015/petersen_graph_view.png

One can also provide the filename where to save the file in which case the file is not open in a viewer:

sage: _ = t.pdf('petersen_graph.pdf')

Another example with polyhedron code taken from this Sage thematic tutorial Draw polytopes in LateX using TikZ:

sage: V = [[1,0,1],[1,0,0],[1,1,0],[0,0,-1],[0,1,0],[-1,0,0],[0,1,1],[0,0,1],[0,-1,0]]
sage: P = Polyhedron(vertices=V).polar()
sage: s = P.projection().tikz([674,108,-731],112)
sage: t = TikzPicture(s)
sage: t
\documentclass[tikz]{standalone}
\begin{document}
\begin{tikzpicture}%
        [x={(0.249656cm, -0.577639cm)},
        y={(0.777700cm, -0.358578cm)},
        z={(-0.576936cm, -0.733318cm)},
        scale=2.000000,
...
... 80 lines not printed (4889 characters in total) ...
...
\node[vertex] at (1.00000, 1.00000, -1.00000)     {};
\node[vertex] at (1.00000, 1.00000, 1.00000)     {};
%%
%%
\end{tikzpicture}
\end{document}
sage: _ = t.pdf()
/Files/2015/polyhedron.png
Read and Post Comments

There are 13.366.431.646 solutions to the Quantumino game

21 septembre 2015 | Catégories: sage | View Comments

Some years ago, I wrote code in Sage to solve the Quantumino puzzle. I also used it to make a one-minute video illustrating the Dancing links algorithm which I am proud to say it is now part of the Dancing links wikipedia page.

/Files/2015/Quantumino.png

Let me recall that the goal of the Quantumino puzzle is to fill a \(2\times 5\times 8\) box with 16 out of 17 three-dimensional pentaminos. After writing the sage code to solve the puzzle, one question was left: how many solutions are there? Is the official website realist or very prudent when they say that there are over 10.000 potential solutions? Can it be computed in hours? days? months? years? The only thing I knew was that the following computation (letting the 0-th pentamino aside) never finished on my machine:

sage: from sage.games.quantumino import QuantuminoSolver
sage: QuantuminoSolver(0).number_of_solutions()   # long time :)

Since I spent already too much time on this side-project, I decided in 2012 to stop investing any more time on it and to really focus on finishing writing my thesis.

So before I finish writing my thesis, I knew that the computation was not going to take a light-year, since I was able to finish the computation of the number of solutions when the 0-th pentamino is put aside and when one pentamino is pre-positioned somewhere in the box. That computation completed in 4 hours on my old laptop and gave about 5 millions solutions. There are 17 choices of pentatminos to put aside, there are 360 distinct positions of that pentamino in the box, so I estimated the number of solution to be something like \(17\times 360\times 5000000 = 30 \times 10^9\). Most importantly, I estimated the computation to take \(17\times 360\times 4= 24480\) hours or 1020 days. Therefore, I knew I could not do it on my laptop.

But last year, I received an email from the designer of the Quantumino puzzle:

-------- Message transféré --------
Sujet : quantumino
Date : Tue, 09 Dec 2014 13:22:30 +0100
De : Nicolaas Neuwahl
Pour : Sebastien Labbe

hi sébastien labbé,

i'm the designer of the quantumino puzzle.
i'm not a mathematician, i'm an architect. i like mathematics.
i'm quite impressed to see the sage work on quantumino, also i have not the
knowledge for full understanding.

i have a question for you - can you tell me HOW MANY different quantumino-
solutions exist?

ty and bye

nicolaas neuwahl

This summer was a good timing to launch the computation on my beautiful Intel® Core™ i5-4590 CPU @ 3.30GHz × 4 at Université de Liège. First, I improved the Sage code to allow a parallel computation of number of solutions in the dancing links code (#18987, merged in a Sage 6.9.beta6). Secondly, we may remark that each tiling of the \(2\times 5\times 8\) box can be rotated in order to find 3 other solutions. It is possible to gain a factor 4 by avoiding to count 4 times the same solution up to rotations (#19107, still needs work from myself). Thanks to Vincent Delecroix for doing the review on both ticket. Dividing the estimated 1024 days of computation needed by a factor \(4\times 4=16\) gives an approximation of 64 days to complete the computation. Two months, just enough to be tractable!

With those two tickets (some previous version to be honest) on top of sage-6.8, I started the computation on August 4th and the computation finished last week on September 18th for a total of 45 days. The computation was stopped only once on September 8th (I forgot to close firefox and thunderbird that night...).

The number of solutions and computation time for each pentamino put aside together with the first solution found is shown in the table below. We remark that some values are equal when the aside pentaminoes are miror images (why!?:).

/Files/2015/b0.png /Files/2015/b1.png
634 900 493 solutions 634 900 493 solutions
2 days, 6:22:44.883358 2 days, 6:19:08.945691
/Files/2015/b2.png /Files/2015/b3.png
509 560 697 solutions 509 560 697 solutions
2 days, 0:01:36.844612 2 days, 0:41:59.447773
/Files/2015/b4.png /Files/2015/b5.png
628 384 422 solutions 628 384 422 solutions
2 days, 7:52:31.459247 2 days, 8:44:49.465672
/Files/2015/b6.png /Files/2015/b7.png
1 212 362 145 solutions 1 212 362 145 solutions
3 days, 17:25:00.346627 3 days, 19:10:02.353063
/Files/2015/b8.png /Files/2015/b9.png
197 325 298 solutions 556 534 800 solutions
22:51:54.439932 1 day, 19:05:23.908326
/Files/2015/b10.png /Files/2015/b11.png
664 820 756 solutions 468 206 736 solutions
2 days, 8:48:54.767662 1 day, 20:14:56.014557
/Files/2015/b12.png /Files/2015/b13.png
1 385 955 043 solutions 1 385 955 043 solutions
4 days, 1:40:30.270929 4 days, 4:44:05.399367
/Files/2015/b14.png /Files/2015/b15.png
694 998 374 solutions 694 998 374 solutions
2 days, 11:44:29.631 2 days, 6:01:57.946708
/Files/2015/b16.png  
1 347 221 708 solutions  
3 days, 21:51:29.043459  

Therefore the total number of solutions up to rotations is 13 366 431 646 which is indeed more than 10000:)

sage: L = [634900493, 634900493, 509560697, 509560697, 628384422,
628384422, 1212362145, 1212362145, 197325298, 556534800, 664820756,
468206736, 1385955043, 1385955043, 694998374, 694998374, 1347221708]
sage: sum(L)
13366431646
sage: factor(_)
2 * 23 * 271 * 1072231
Summary
The machine (4 cores) Intel® Core™ i5-4590 CPU @ 3.30GHz × 4 (Université de Liège)
Computation Time 45 days, (Aug 4th -- Sep 18th, 2015)
Number of solutions (up to rotations) 13 366 431 646
Number of solutions / cpu / second 859

My code will be available on github.

About the video on wikipedia.

I must say that the video is not perfect. On wikipedia, the file talk page of the video says that the Jerky camera movement is distracting. That is because I managed to make the video out of images created by .show(viewer='tachyon') which changes the coordinate system, hardcodes a lot of parameters, zoom properly, simplifies stuff to make sure the user don't see just a blank image. But, for making a movie, we need access to more parameters especially the placement of the camera (to avoid the jerky movement). I know that Tachyon allows all of that. It is still a project that I have to create a more versatile Graphics3D -> Tachyon conversion allowing to construct nice videos of evolving mathematical objects. That's another story.

Read and Post Comments

Arnoux-Rauzy-Poincaré sequences

26 février 2015 | Catégories: sage | View Comments

In a recent article with Valérie Berthé [BL15], we provided a multidimensional continued fraction algorithm called Arnoux-Rauzy-Poincaré (ARP) to construct, given any vector \(v\in\mathbb{R}_+^3\), an infinite word \(w\in\{1,2,3\}^\mathbb{N}\) over a three-letter alphabet such that the frequencies of letters in \(w\) exists and are equal to \(v\) and such that the number of factors (i.e. finite block of consecutive letters) of length \(n\) appearing in \(w\) is linear and less than \(\frac{5}{2}n+1\). We also conjecture that for almost all \(v\) the contructed word describes a discrete path in the positive octant staying at a bounded distance from the euclidean line of direction \(v\).

In Sage, you can construct this word using the next version of my package slabbe-0.2 (not released yet, email me to press me to finish it). The one with frequencies of letters proportionnal to \((1, e, \pi)\) is:

sage: from slabbe.mcf import algo
sage: D = algo.arp.substitutions()
sage: it = algo.arp.coding_iterator((1,e,pi))
sage: w = words.s_adic(it, repeat(1), D)
word: 1232323123233231232332312323123232312323...

The factor complexity is close to 2n+1 and the balance is often less or equal to three:

sage: w[:10000].number_of_factors(100)
202
sage: w[:100000].number_of_factors(1000)
2002
sage: w[:1000].balance()
3
sage: w[:2000].balance()
3

Note that bounded distance from the euclidean line almost surely was proven in [DHS2013] for Brun algorithm, another MCF algorithm.

Other approaches: Standard model and billiard sequences

Other approaches have been proposed to construct such discrete lines.

One of them is the standard model of Eric Andres [A03]. It is also equivalent to billiard sequences in the cube. It is well known that the factor complexity of billiard sequences is quadratic \(p(n)=n^2+n+1\) [AMST94]. Experimentally, we can verify this. We first create a billiard word of some given direction:

sage: from slabbe import BilliardCube
sage: v = vector(RR, (1, e, pi))
sage: b = BilliardCube(v)
sage: b
Cubic billiard of direction (1.00000000000000, 2.71828182845905, 3.14159265358979)
sage: w = b.to_word()
sage: w
word: 3231232323123233213232321323231233232132...

We create some prefixes of \(w\) that we represent internally as char*. The creation is slow because the implementation of billiard words in my optional package is in Python and is not that efficient:

sage: p3 = Word(w[:10^3], alphabet=[1,2,3], datatype='char')
sage: p4 = Word(w[:10^4], alphabet=[1,2,3], datatype='char') # takes 3s
sage: p5 = Word(w[:10^5], alphabet=[1,2,3], datatype='char') # takes 32s
sage: p6 = Word(w[:10^6], alphabet=[1,2,3], datatype='char') # takes 5min 20s

We see below that exactly \(n^2+n+1\) factors of length \(n<20\) appears in the prefix of length 1000000 of \(w\):

sage: A = ['n'] + range(30)
sage: c3 = ['p_(w[:10^3])(n)'] + map(p3.number_of_factors, range(30))
sage: c4 = ['p_(w[:10^4])(n)'] + map(p4.number_of_factors, range(30))
sage: c5 = ['p_(w[:10^5])(n)'] + map(p5.number_of_factors, range(30)) # takes 4s
sage: c6 = ['p_(w[:10^6])(n)'] + map(p6.number_of_factors, range(30)) # takes 49s
sage: ref = ['n^2+n+1'] + [n^2+n+1 for n in range(30)]
sage: T = table(columns=[A,c3,c4,c5,c6,ref])
sage: T
  n    p_(w[:10^3])(n)   p_(w[:10^4])(n)   p_(w[:10^5])(n)   p_(w[:10^6])(n)   n^2+n+1
+----+-----------------+-----------------+-----------------+-----------------+---------+
  0    1                 1                 1                 1                 1
  1    3                 3                 3                 3                 3
  2    7                 7                 7                 7                 7
  3    13                13                13                13                13
  4    21                21                21                21                21
  5    31                31                31                31                31
  6    43                43                43                43                43
  7    52                55                56                57                57
  8    63                69                71                73                73
  9    74                85                88                91                91
  10   87                103               107               111               111
  11   100               123               128               133               133
  12   115               145               151               157               157
  13   130               169               176               183               183
  14   144               195               203               211               211
  15   160               223               232               241               241
  16   176               253               263               273               273
  17   192               285               296               307               307
  18   208               319               331               343               343
  19   224               355               368               381               381
  20   239               392               407               421               421
  21   254               430               448               463               463
  22   268               470               491               507               507
  23   282               510               536               553               553
  24   296               552               583               601               601
  25   310               596               632               651               651
  26   324               642               683               703               703
  27   335               687               734               757               757
  28   345               734               787               813               813
  29   355               783               842               871               871

Billiard sequences generate paths that are at a bounded distance from an euclidean line. This is equivalent to say that the balance is finite. The balance is defined as the supremum value of difference of the number of apparition of a letter in two factors of the same length. For billiard sequences, the balance is 2:

sage: p3.balance()
2
sage: p4.balance() # takes 2min 37s
2

Other approaches: Melançon and Reutenauer

Melançon and Reutenauer [MR13] also suggested a method that generalizes Christoffel words in higher dimension. The construction is based on the application of two substitutions generalizing the construction of sturmian sequences. Below we compute the factor complexity and the balance of some of their words over a three-letter alphabet.

On a three-letter alphabet, the two morphisms are:

sage: L = WordMorphism('1->1,2->13,3->2')
sage: R = WordMorphism('1->13,2->2,3->3')
sage: L
WordMorphism: 1->1, 2->13, 3->2
sage: R
WordMorphism: 1->13, 2->2, 3->3

Example 1: periodic case \(LRLRLRLRLR\dots\). In this example, the factor complexity seems to be around \(p(n)=2.76n\) and the balance is at least 28:

sage: from itertools import repeat, cycle
sage: W = words.s_adic(cycle((L,R)),repeat('1'))
sage: W
word: 1213122121313121312212212131221213131213...
sage: map(W[:10000].number_of_factors, [10,20,40,80])
[27, 54, 110, 221]
sage: [27/10., 54/20., 110/40., 221/80.]
[2.70000000000000, 2.70000000000000, 2.75000000000000, 2.76250000000000]
sage: W[:1000].balance()  # takes 1.6s
21
sage: W[:2000].balance()  # takes 6.4s
28

Example 2: \(RLR^2LR^4LR^8LR^{16}LR^{32}LR^{64}LR^{128}\dots\) taken from the conclusion of their article. In this example, the factor complexity seems to be \(p(n)=3n\) and balance at least as high (=bad) as \(122\):

sage: W = words.s_adic([R,L,R,R,L,R,R,R,R,L]+[R]*8+[L]+[R]*16+[L]+[R]*32+[L]+[R]*64+[L]+[R]*128,'1')
sage: W.length()
330312
sage: map(W.number_of_factors, [10, 20, 100, 200, 300, 1000])
[29, 57, 295, 595, 895, 2981]
sage: [29/10., 57/20., 295/100., 595/200., 895/300., 2981/1000.]
[2.90000000000000,
 2.85000000000000,
 2.95000000000000,
 2.97500000000000,
 2.98333333333333,
 2.98100000000000]
sage: W[:1000].balance()  # takes 1.6s
122
sage: W[:2000].balance()  # takes 6s
122

Example 3: some random ones. The complexity \(p(n)/n\) occillates between 2 and 3 for factors of length \(n=1000\) in prefixes of length 100000:

sage: for _ in range(10):
....:     W = words.s_adic([choice((L,R)) for _ in range(50)],'1')
....:     print W[:100000].number_of_factors(1000)/1000.
2.02700000000000
2.23600000000000
2.74000000000000
2.21500000000000
2.78700000000000
2.52700000000000
2.85700000000000
2.33300000000000
2.65500000000000
2.51800000000000

For ten randomly generated words, the balance goes from 6 to 27 which is much more than what is obtained for billiard words or by our approach:

sage: for _ in range(10):
....:     W = words.s_adic([choice((L,R)) for _ in range(50)],'1')
....:     print W[:1000].balance(), W[:2000].balance()
12 15
8 24
14 14
5 11
17 17
14 14
6 6
19 27
9 16
12 12

References

[BL15]V. Berthé, S. Labbé, Factor Complexity of S-adic words generated by the Arnoux-Rauzy-Poincaré Algorithm, Advances in Applied Mathematics 63 (2015) 90-130. http://dx.doi.org/10.1016/j.aam.2014.11.001
[DHS2013]Delecroix, Vincent, Tomás Hejda, and Wolfgang Steiner. “Balancedness of Arnoux-Rauzy and Brun Words.” In Combinatorics on Words, 119–31. Springer, 2013. http://link.springer.com/chapter/10.1007/978-3-642-40579-2_14.
[A03]E. Andres, Discrete linear objects in dimension n: the standard model, Graphical Models 65 (2003) 92-111.
[AMST94]P. Arnoux, C. Mauduit, I. Shiokawa, J. I. Tamura, Complexity of sequences defined by billiards in the cube, Bull. Soc. Math. France 122 (1994) 1-12.
[MR13]G. Melançon, C. Reutenauer, On a class of Lyndon words extending Christoffel words and related to a multidimensional continued fraction algorithm. J. Integer Seq. 16, No. 9, Article 13.9.7, 30 p., electronic only (2013). https://cs.uwaterloo.ca/journals/JIS/VOL16/Reutenauer/reut3.html
Read and Post Comments

Next Page »