A time evolution picture of packages built in parallel by Sage
16 décembre 2016 | Catégories: sage | View CommentsCompiling 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 CommentsYesterday 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 CommentsThese 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).
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 succesful 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 useful 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 CommentsSome 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 CommentsIn 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 |
« Previous Page -- Next Page »
Propulsé par Blogofile
S'abonner au Flux RSS
et aux Commentaires.
This work by Sébastien Labbé is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.