## 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()

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()

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.

## 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]

## 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.pylanguage.pylyapunov.pymatrix_cocycle.pymult_cont_frac.pyxranking_scale.pytikz_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).

**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')

Drawing the cylinders:

sage: cocycle = algo.matrix_cocycle() sage: t = cocycle.tikz_n_cylinders(3, scale=3) sage: t.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()

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)

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()

## 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.

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!?:).

634 900 493 solutions | 634 900 493 solutions |

2 days, 6:22:44.883358 | 2 days, 6:19:08.945691 |

509 560 697 solutions | 509 560 697 solutions |

2 days, 0:01:36.844612 | 2 days, 0:41:59.447773 |

628 384 422 solutions | 628 384 422 solutions |

2 days, 7:52:31.459247 | 2 days, 8:44:49.465672 |

1 212 362 145 solutions | 1 212 362 145 solutions |

3 days, 17:25:00.346627 | 3 days, 19:10:02.353063 |

197 325 298 solutions | 556 534 800 solutions |

22:51:54.439932 | 1 day, 19:05:23.908326 |

664 820 756 solutions | 468 206 736 solutions |

2 days, 8:48:54.767662 | 1 day, 20:14:56.014557 |

1 385 955 043 solutions | 1 385 955 043 solutions |

4 days, 1:40:30.270929 | 4 days, 4:44:05.399367 |

694 998 374 solutions | 694 998 374 solutions |

2 days, 11:44:29.631 | 2 days, 6:01:57.946708 |

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

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.

## 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 |

Next Page »