From d7d44558e95be693516baa92b9455a650ae78fb8 Mon Sep 17 00:00:00 2001 From: Paul Mercat Date: Tue, 5 Aug 2014 18:46:49 +0200 Subject: [PATCH] Add a better way to draw a beta-adic monoid (new function plot2) using extern C files. I had to change the automaton associated with the beta-adic monoid for WordMorphism. --- src/sage/combinat/words/automata.pyx | 26 +- src/sage/combinat/words/morphism.py | 25 +- src/sage/monoids/beta_adic_monoid.pyx | 236 ++++++++++++++++ src/sage/monoids/complex.h | 57 ++++ src/sage/monoids/draw.c | 388 ++++++++++++++++++++++++++ src/sage/monoids/draw.h | 41 +++ 6 files changed, 759 insertions(+), 14 deletions(-) create mode 100644 src/sage/monoids/complex.h create mode 100644 src/sage/monoids/draw.c create mode 100644 src/sage/monoids/draw.h diff --git a/src/sage/combinat/words/automata.pyx b/src/sage/combinat/words/automata.pyx index bd977e5810f..c0038c4c4e8 100644 --- a/src/sage/combinat/words/automata.pyx +++ b/src/sage/combinat/words/automata.pyx @@ -67,14 +67,32 @@ class Automaton (DiGraph): """ return DiGraph.plot(self, edge_labels=True) - def plot2 (self, name="Automaton"): + def plot2 (self, name="Automaton", verb=False): r""" Plot the automaton. """ + I = [] + if hasattr(self, 'I'): + I = self.I + F = [] + if hasattr(self, 'F'): + F = self.F + if verb: + print I,F txt = '{\n' - #for e in self.vertices(): - # txt += ' "'+str(e)+'"\n' - #txt += '\n' + for e in self.vertices(): + txt += ' "'+str(e)+'" [shape=' + if e in F: + txt += 'doublecircle' + else: + txt += 'circle' + txt += ', style=' + if e in I: + txt += 'bold' + else: + txt += 'solid' + txt += ', fontsize=20, margin=0]\n' + txt += '\n' for e,f,l in self.edges(): txt += ' "'+str(e)+'" -> "'+str(f)+'" [label="'+str(l)+'"];\n' txt += '}\n' diff --git a/src/sage/combinat/words/morphism.py b/src/sage/combinat/words/morphism.py index 393623f0b40..91a4450088a 100644 --- a/src/sage/combinat/words/morphism.py +++ b/src/sage/combinat/words/morphism.py @@ -3088,7 +3088,7 @@ def abelian_rotation_subspace(self): return M._column_ambient_module().change_ring(QQ).subspace(basis) - def rauzy_fractal_beta_adic_monoid (self, I=None, eig=None, invss=False, verb=False): + def rauzy_fractal_beta_adic_monoid (self, I=None, eig=None, invss=True, det=True, verb=False): r""" Returns the beta-adic monoid associated to the substitution. @@ -3116,26 +3116,31 @@ def rauzy_fractal_beta_adic_monoid (self, I=None, eig=None, invss=False, verb=Fa for a2 in self._morph[a]: C[i] = s i=i+1 - ss.add_edge(a2, a, s) + ss.add_edge(a, a2, s) s += C0[a2] if verb: print "C=%s"%C - C = Set(C.values()) #ensemble des polynômes en beta qui sont des chiffres + C = set(C.values()) #ensemble des polynômes en beta qui sont des chiffres from sage.monoids.beta_adic_monoid import BetaAdicMonoid res = BetaAdicMonoid(K.gen(), C) - ss.A = C - ss.I = ss.vertices() + ss.F = ss.vertices() if I is not None: - ss.F = I + ss.I = I else: - ss.F = ss.vertices() + ss.I = [ss.vertices()[0]] + ss.A = C if invss: - res.ss = ss.transpose() - else: - res.ss = ss.determinize(noempty=True) + ss = ss.transpose() + + if verb: + print ss + + if det: + ss = ss.determinize(noempty=True) + res.ss = ss return res diff --git a/src/sage/monoids/beta_adic_monoid.pyx b/src/sage/monoids/beta_adic_monoid.pyx index 5e51f9120f8..8a65daa42db 100644 --- a/src/sage/monoids/beta_adic_monoid.pyx +++ b/src/sage/monoids/beta_adic_monoid.pyx @@ -54,6 +54,136 @@ def emonde (a, K): if K.zero() in s: return s +cdef extern from "draw.c": + cdef cppclass Etat: + int* f + int final + cdef cppclass Automate: + int n + int na + Etat* e + int i + ctypedef unsigned char uint8 + cdef cppclass Color: + uint8 r + uint8 g + uint8 b + uint8 a + cdef cppclass Surface: + Color **pix + int sx, sy + cdef cppclass Complexe: + double x + double y + cdef cppclass BetaAdic: + Complexe b + Complexe *t #liste des translations + int n #nombre de translations + Automate a + + Surface NewSurface (int sx, int sy) + void FreeSurface (Surface s) + Automate NewAutomate (int n, int na) + void FreeAutomate(Automate a) + BetaAdic NewBetaAdic (int n) + void FreeBetaAdic (BetaAdic b) + void Draw (BetaAdic b, Surface s, int n, int ajust, Color col, int verb) + void Draw2 (BetaAdic b, Surface s, int n, int ajust, Color col, int verb) + void print_word (BetaAdic b, int n, int etat) + +cdef Complexe complex (c): + cdef Complexe r + r.x = c.real() + r.y = c.imag() + return r + +cdef surface_to_img (Surface s): + import numpy as np + from PIL import Image + arr = np.zeros([s.sy, s.sx], dtype = [('r', 'uint8'), ('g', 'uint8'), ('b', 'uint8'), ('a', 'uint8')]) + cdef int x, y + cdef Color c + for x in range(s.sx): + for y in range(s.sy): + c = s.pix[x][s.sy -y-1] + arr[y,x][0] = c.r + arr[y,x][1] = c.g + arr[y,x][2] = c.b + arr[y,x][3] = c.a + img = Image.fromarray(arr, 'RGBA') + img.save("/Users/mercat/Desktop/output.png") + img.save("output.png") + +cdef Automate getAutomate (a, d, iss=None, verb=False): + if verb: + print "getAutomate %s..."%a + lv = a.vertices() + if hasattr(a, 'F'): + F = a.F + else: + F = lv + #alloue l'automate + cdef Automate r = NewAutomate(a.num_verts(), len(a.Alphabet())) + #réindice les sommets + dv = {} + cdef int i + for u,i in zip(lv, range(len(lv))): + dv[u] = i + if u in F: + r.e[i].final = 1 + if verb: + print len(lv) + #copie l'automate en C + le = a.edges() + if verb: + print "len(le)=%s"%len(le) + for u,v,l in le: + if d.has_key(l): + #if dv.has_key(u) and dv.has_key(v): + r.e[dv[u]].f[d[l]] = dv[v] + #else: + # print "Erreur : pas de clef %s ou %s !"%(u,v) + else: + print "Erreur : pas de clef %s !"%l + if verb: + print "I..." + if iss is not None: + r.i = iss + else: + if hasattr(a, 'I') and len(a.I) > 0: + r.i = dv[list(a.I)[0]] + else: + raise ValueError("The initial state must be defined !") + if verb: + print "...getAutomate" + return r + +cdef BetaAdic getBetaAdic (self, prec=53, ss=None, iss=None, add_letters=True, verb=False): + from sage.rings.complex_field import ComplexField + CC = ComplexField(prec) + cdef BetaAdic b + if ss is None: + if hasattr(self, 'ss'): + ss = self.ss + else: + ss = self.default_ss() + C = set(self.C) + if add_letters: + C.update(ss.Alphabet()) + b = NewBetaAdic(len(C)) + b.b = complex(CC(self.b)) + d = {} + for i,c in zip(range(b.n), C): + b.t[i] = complex(CC(c)) + d[c] = i + #automaton + b.a = getAutomate(ss, d, iss=iss, verb=verb) + return b + +def PrintWord (m, n): + b = getBetaAdic(m) + print_word(b, n, b.a.i) + class BetaAdicMonoid(Monoid_class): r""" ``b``-adic monoid with numerals set ``C``. @@ -333,6 +463,112 @@ class BetaAdicMonoid(Monoid_class): # orbit_points = orbit_points.union(Set([place(b)*p+place(c) for p in orbit_points0])) # return orbit_points + def plot2 (self, n=None, ss=None, iss=None, sx=800, sy=600, ajust=True, prec=53, color=(0,0,0,255), method=0, add_letters=True, verb=False): + r""" + Draw the limit set of the beta-adic monoid (with or without subshift). + + INPUT: + + - ``n`` - integer (default: ``None``) + The number of iterations used to plot the fractal. + Default values: between ``5`` and ``16`` depending on the number of generators. + + - ``place`` - place of the number field of beta (default: ``None``) + The place we should use to evaluate elements of the number field. + + - ``ss`` - Automaton (default: ``None``) + The subshift to associate to the beta-adic monoid for this drawing. + + - ``iss`` - set of initial states of the automaton ss (default: ``None``) + + - ``sx, sy`` - dimensions of the resulting image (default : ``800, 600``) + + - ``ajust`` - adapt the drawing to fill all the image, with ratio 1 (default: ``True``) + + - ``prec`` - precision of returned values (default: ``53``) + + - ``color`` - list of three integer between 0 and 255 (default: ``(0,0,255,255)``) + Color of the drawing. + + - ``verb`` - bool (default: ``False``) + Print informations for debugging. + + OUTPUT: + + A Graphics object. + + EXAMPLES: + + #. The dragon fractal:: + + sage: m=BetaAdicMonoid(1/(1+I), {0,1}) + sage: m.plot2() # long time + + #. The Rauzy fractal of the Tribonacci substitution:: + + sage: s = WordMorphism('1->12,2->13,3->1') + sage: m = s.rauzy_fractal_beta_adic_monoid() + sage: m.plot2() # long time + + #. A non-Pisot Rauzy fractal:: + + sage: s = WordMorphism({1:[3,2], 2:[3,3], 3:[4], 4:[1]}) + sage: m = s.rauzy_fractal_beta_adic_monoid() + sage: m.b = 1/m.b + sage: m.ss = m.ss.transpose().determinize().minimize() + sage: m.plot2() # long time + + #. The dragon fractal and its boundary:: + + sage: m = BetaAdicMonoid(1/(1+I), {0,1}) + sage: p1 = m.plot() # long time + sage: ssi = m.intersection_words(w1=[0], w2=[1]) # long time + sage: p2 = m.plot2(ss = ssi, n=18) # long time + sage: p1+p2 # long time + + #. The "Hokkaido" fractal and its boundary:: + + sage: s = WordMorphism('a->ab,b->c,c->d,d->e,e->a') + sage: m = s.rauzy_fractal_beta_adic_monoid() + sage: p1 = m.plot() # long time + sage: ssi = m.intersection_words(w1=[0], w2=[1]) # long time + sage: p2 = m.plot2(ss=ssi, n=40) # long time + sage: p1+p2 # long time + + #. A limit set that look like a tiling:: + + sage: P=x^4 + x^3 - x + 1 + sage: b = P.roots(ring=QQbar)[2][0] + sage: m = BetaAdicMonoid(b, {0,1}) + sage: m.plot2(18) # long time + + """ + cdef Surface s = NewSurface (sx, sy) + cdef BetaAdic b + b = getBetaAdic(self, prec=prec, ss=ss, iss=iss, add_letters=add_letters, verb=verb) + #dessin + cdef Color col + col.r = color[0] + col.g = color[1] + col.b = color[2] + col.a = color[3] + if n is None: + n = -1 + if method == 0: + Draw(b, s, n, ajust, col, verb) + elif method == 1: + print "Not implemented !" + return + #lv = s.rauzy_fractal_projection_exact().values() + #for i,v in zip(range(len(lv)),lv): + # b.t[i] = complex(CC(v)) + #Draw2(b, s, n, ajust, col, verb) + #enregistrement du résultat + surface_to_img(s) + FreeSurface(s) + FreeAutomate(b.a) + FreeBetaAdic(b) + def plot (self, n=None, place=None, ss=None, iss=None, prec=53, point_size=None, color='blue', verb=False): r""" Draw the limit set of the beta-adic monoid (with or without subshift). diff --git a/src/sage/monoids/complex.h b/src/sage/monoids/complex.h new file mode 100644 index 00000000000..3623879d00f --- /dev/null +++ b/src/sage/monoids/complex.h @@ -0,0 +1,57 @@ + +struct Complexe +{ + double x,y; +}; +typedef struct Complexe Complexe; + +Complexe prod (Complexe a, Complexe b) +{ + Complexe r; + r.x = a.x*b.x - a.y*b.y; + r.y = a.x*b.y + a.y*b.x; + return r; +} + +Complexe zero () +{ + Complexe r; + r.x = 0; + r.y = 0; + return r; +} + +Complexe un () +{ + Complexe r; + r.x = 1; + r.y = 0; + return r; +} + +Complexe add (Complexe a, Complexe b) +{ + Complexe r; + r.x = a.x + b.x; + r.y = a.y + b.y; + return r; +} + +inline double carre (double x) +{ + return x*x; +} + +double cnorm (Complexe c) +{ + return carre(c.x) + carre(c.y); +} + +Complexe inv (Complexe c) +{ + Complexe r; + double cn = cnorm(c); + r.x = c.x/cn; + r.y = -c.y/cn; + return r; +} diff --git a/src/sage/monoids/draw.c b/src/sage/monoids/draw.c new file mode 100644 index 00000000000..48a04133c21 --- /dev/null +++ b/src/sage/monoids/draw.c @@ -0,0 +1,388 @@ +#include +#include "complex.h" +#include "draw.h" + +double mx = -2, my = -2, Mx = 2, My = 2; //zone de dessin + +double mx2 = 1000000, my2 = 1000000, Mx2 = -1000000, My2 = -1000000; //extremum observés + +Color color0; //couleur du fond +Color color; //couleur de dessin +Color* colors; //liste de couleurs de dessin + +Surface NewSurface (int sx, int sy) +{ + Surface s; + s.sx = sx; + s.sy = sy; + s.pix = (Color **)malloc(sizeof(Color *)*sx); + int x; + for (x=0;x= Mx || p.y < my || p.y >= My) + return; + double fx = (p.x - mx)*s.sx/(Mx-mx); + double fy = (p.y - my)*s.sy/(My-my); + unsigned int x = fx; + unsigned int y = fy; + if (x < s.sx && y < s.sy) + { + if (x+1 < s.sx && y+1 < s.sy) + { + s.pix[x][y] = moy(s.pix[x][y], color, (1.-fx+x)*(1.-fy+y)); + s.pix[x+1][y] = moy(s.pix[x+1][y], color, (fx-x)*(1.-fy+y)); + s.pix[x][y+1] = moy(s.pix[x][y+1], color, (1.-fx+x)*(fy-y)); + s.pix[x+1][y+1] = moy(s.pix[x+1][y+1], color, (fx-x)*(fy-y)); + }else + s.pix[x][y] = color; + } +} + +void print_word (BetaAdic b, int n, int etat) +{ + if (etat < 0) + return; + if (n == 0) + printf("(%d)", etat); + else + { + int i; + for (i=0;i= 0 && etat < b.a.n) + { + if (!b.a.e[etat].final) + { + //printf("%d pas final !", etat); + return; + } + }else + { + printf("état %d !\n", etat); + return; + } + if (pos.x < mx2) + mx2 = pos.x; + if (pos.x > Mx2) + Mx2 = pos.x; + if (pos.y < my2) + my2 = pos.y; + if (pos.y > My2) + My2 = pos.y; + pos = add(pos, b.t[etat]); + set_pix (s, pos); + }else + { + int i; + for (i=0;i= 0 && etat < b.a.n) + { + if (!b.a.e[etat].final) + { + //printf("%d pas final !", etat); + return; + } + }else + { + printf("état %d !\n", etat); + return; + } + if (p.x < mx2) + mx2 = p.x; + if (p.x > Mx2) + Mx2 = p.x; + if (p.y < my2) + my2 = p.y; + if (p.y > My2) + My2 = p.y; + set_pix (s, p); + }else + { + int i; + for (i=0;i 0) + { + My = My + delta/(2*s.sx); + my = my - delta/(2*s.sx); + }else + { + Mx = Mx - delta/(2*s.sy); + mx = mx + delta/(2*s.sy); + } + } + if (verb) + { + printf("Zone de dessin : (%lf, %lf) (%lf, %lf)\n", mx, my, Mx, My); + } + Fill(s, color0); + if (auto_n) + { + n = .5 - 2.*log(max(3.*s.sx/(Mx-mx), 3.*s.sy/(My-my)))/log(cnorm(b.b)); + if (verb) + printf("n = %d\n", n); + } + //niter = n; + Draw_rec (b, s, n, zero(), un(), b.a.i); +} + +void Draw2 (BetaAdic b, Surface s, int n, int ajust, Color col, int verb) +{ + int auto_n = (n < 0); + //set global variables + mx2 = 1000000, my2 = 1000000, Mx2 = -1000000, My2 = -1000000; //extremum observés + color0.r = color0.g = color0.b = 255; + color0.a = 0; + color = col; + colors = (Color *)malloc(sizeof(Color)*b.a.n); + int i; + for (i=0;i 0) + { + My = My + delta/(2*s.sx); + my = my - delta/(2*s.sx); + }else + { + Mx = Mx - delta/(2*s.sy); + mx = mx + delta/(2*s.sy); + } + } + if (verb) + { + printf("Zone de dessin : (%lf, %lf) (%lf, %lf)\n", mx, my, Mx, My); + } + Fill(s, color0); + if (auto_n) + { + n = .5 - 2.*log(max(3.*s.sx/(Mx-mx), 3.*s.sy/(My-my)))/log(cnorm(b.b)); + if (verb) + printf("n = %d\n", n); + } + pos = zero(); + Draw_rec2 (b, s, n, b.a.i); +} diff --git a/src/sage/monoids/draw.h b/src/sage/monoids/draw.h new file mode 100644 index 00000000000..2e4e34c5490 --- /dev/null +++ b/src/sage/monoids/draw.h @@ -0,0 +1,41 @@ + +//typedef int *Etat; +struct Etat +{ + int* f; //liste des na fils indexés par les lettres + int final; //l'état est final ? +}; +typedef struct Etat Etat; + +struct Automate +{ + Etat* e; //liste des états (l'état initial est le premier) + int n; //nombre d'états + int na; //nombre de lettres différentes + int i; //état initial +}; +typedef struct Automate Automate; + +typedef unsigned char uint8; + +struct Color +{ + uint8 r,g,b,a; +}; +typedef struct Color Color; + +struct Surface +{ + Color** pix; + int sx, sy; //dimensions de l'image +}; +typedef struct Surface Surface; + +struct BetaAdic +{ + Complexe b; // + Complexe *t; //liste des translations + int n; //nombre de translations + Automate a; +}; +typedef struct BetaAdic BetaAdic;