Skip to content
This repository has been archived by the owner on Jan 30, 2023. It is now read-only.

Commit

Permalink
Add new functionalities using FastAutomaton : copy, intersection, ...…
Browse files Browse the repository at this point in the history 
… and correct some bugs.
  • Loading branch information
@mercatp
mercatp committed Sep 2, 2014
1 parent 9ce37a9 commit c18febb
Show file tree
Hide file tree
Showing 5 changed files with 165 additions and 18 deletions.
35 changes: 31 additions & 4 deletions src/sage/combinat/words/automataC.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,22 @@ void ReallocAutomaton (Automaton *a, int n)
*/
}

Automaton CopyAutomaton (Automaton a)
{
Automaton r = NewAutomaton(a.n, a.na);
int i,j;
for (i=0;i<a.n;i++)
{
r.e[i].final = a.e[i].final;
for (j=0;j<a.na;j++)
{
r.e[i].f[j] = a.e[i].f[j];
}
}
r.i = a.i;
return r;
}

void init (Automaton a)
{
int i,j;
Expand Down Expand Up @@ -1289,6 +1305,9 @@ void emonde_rec3 (Automaton a, Automaton r, int *l, int etat)
*/

//retire tous les états non accessible ou non co-accessible
//
// fonction pas très éfficace : à revoir !!!!!!!!!!!!!!!!!!!!!!!!!!!!
//
Automaton emonde (Automaton a, bool verb)
{
int i,j,f;
Expand Down Expand Up @@ -1364,14 +1383,18 @@ Automaton emonde (Automaton a, bool verb)
}
printf(" ]\n");

printf("create the new automaton...\n");
printf("create the new automaton %d %d...\n", cpt, a.na);
}
//créé le nouvel automate
Automaton r = NewAutomaton(cpt, a.na);
for (i=0;i<a.n;i++)
{
if (l[i] == -1)
{
if (verb)
printf("pass %d\n", i);
continue;
}
for (j=0;j<a.na;j++)
{
f = a.e[i].f[j];
Expand All @@ -1387,7 +1410,10 @@ Automaton emonde (Automaton a, bool verb)

//remet les états finaux de a
if (verb)
{
printf("Etats supprimés : [");
fflush(stdout);
}
for (i=0;i<a.n;i++)
{
if (verb)
Expand All @@ -1403,7 +1429,8 @@ Automaton emonde (Automaton a, bool verb)
}
}
a.e[i].final &= 1;
r.e[l[i]].final = a.e[i].final;
if (l[i] != -1)
r.e[l[i]].final = a.e[i].final;
}
if (verb)
printf(" ]\n");
Expand Down Expand Up @@ -1647,7 +1674,7 @@ void PermutOP (Automaton a, int *l, int na, bool verb)
free(lf);
}

/////////////////////// the following is an implementation of Hopcroft's algorithm
/////////////////////// the following is an implementation of Hopcroft's algorithm minimization

typedef int Couple[2];

Expand Down Expand Up @@ -2019,7 +2046,7 @@ Automaton Minimise (Automaton a, bool verb)
return r;
}

///////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////

/*
int sign (int a)
Expand Down
1 change: 1 addition & 0 deletions src/sage/combinat/words/automataC.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ void printDict (Dict d);
void dictAdd (Dict *d, int e); //ajoute un élément au dictionnaire (même s'il était déjà présent)
Automaton NewAutomaton (int n, int na);
void FreeAutomaton (Automaton *a);
Automaton CopyAutomaton (Automaton a);
void init (Automaton a);
void printAutomaton (Automaton a);
void plotTikZ (Automaton a, const char **labels, const char *graph_name, double sx, double sy);
Expand Down
32 changes: 31 additions & 1 deletion src/sage/combinat/words/cautomata.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ cdef extern from "automataC.h":

Automaton NewAutomaton (int n, int na)
void FreeAutomaton (Automaton *a)
Automaton CopyAutomaton (Automaton a)
void init (Automaton a)
void printAutomaton (Automaton a)
void plotTikZ (Automaton a, const char **labels, const char *graph_name, double sx, double sy)
Expand Down Expand Up @@ -375,13 +376,28 @@ cdef class FastAutomaton:
def initial_state (self):
return self.a.i

def set_initial_state (self, int i):
self.a.i = i

def final_states (self):
l = []
for i in range(self.a.n):
if self.a.e[i].final:
l.append(i)
return l

def states (self):
return range(self.a.n)

def set_final_states (self, lf):
cdef int f
for f in range(self.a.n):
self.a.e[f].final = 0
for f in lf:
if f < 0 or f >= self.a.n:
raise ValueError("%d is not a state !"%f)
self.a.e[f].final = 1

def succ (self, int i, int j):
if i<0 or i>=self.a.n or j<0 or j>=self.a.na:
return -1
Expand Down Expand Up @@ -421,7 +437,12 @@ cdef class FastAutomaton:
return Bool(equalsAutomaton(self.a[0], b.a[0]))

#assume that the dictionnary d is injective !!!
def product (self, FastAutomaton b, dict d, verb=False):
def product (self, FastAutomaton b, dict d=None, verb=False):
if d is None:
d = {}
for la in self.A:
for lb in b.A:
d[(la,lb)] = (la,lb)
sig_on()
cdef Automaton a
cdef Dict dC
Expand Down Expand Up @@ -659,3 +680,12 @@ cdef class FastAutomaton:

def test (self):
Test()

def copy (self):
sig_on()
r = FastAutomaton(None)
r.a[0] = CopyAutomaton(self.a[0])
r.A = self.A
sig_off()
return r

73 changes: 66 additions & 7 deletions src/sage/monoids/beta_adic_monoid.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,10 +104,13 @@ cdef getElement (e, Element r, int n):
for j in range(n):
r.c[j] = p[j]

cdef InfoBetaAdic initInfoBetaAdic (self, Cd=None, verb=False):
cdef InfoBetaAdic initInfoBetaAdic (self, Cd=None, verb=False) except *:
#compute all the data in sage
b = self.b
K = b.parent()
K = NumberField((1/self.b).minpoly(), 'b')
b = K.gen()
C = [c.lift()(1/b) for c in self.C]
if verb:
print "C = %s"%C

if verb: print K

Expand All @@ -119,7 +122,9 @@ cdef InfoBetaAdic initInfoBetaAdic (self, Cd=None, verb=False):
pi = K.defining_polynomial()
from sage.rings.arith import gcd
pi = pi/gcd(pi.list()) #rend le polynôme à coefficients entiers et de contenu 1
if verb: print "pi=%s"%pi
lp = (Integer(pi.list()[0])).prime_divisors() #liste des nombres premiers concernés
if verb: print "lp=%s"%lp
pultra = [] #liste des places ultramétriques considérées
for p in lp:
#détermine toutes les places au dessus de p dans le corps de nombres K
Expand All @@ -145,7 +150,9 @@ cdef InfoBetaAdic initInfoBetaAdic (self, Cd=None, verb=False):

#calcule les bornes max pour chaque valeur absolue
if Cd is None:
Cd = Set([c-c2 for c in self.C for c2 in self.C])
Cd = Set([c-c2 for c in C for c2 in C])
else:
Cd = [c.lift()(1/b) for c in Cd]
if verb: print "Cd = %s"%Cd
m = dict([])
for p in parch:
Expand Down Expand Up @@ -174,10 +181,10 @@ cdef InfoBetaAdic initInfoBetaAdic (self, Cd=None, verb=False):
initCdInfoBetaAdic(self, &i, Cd=Cd, verb=verb)
return i

cdef initCdInfoBetaAdic (self, InfoBetaAdic *i, Cd=None, verb=False):
cdef initCdInfoBetaAdic (self, InfoBetaAdic *i, Cd, verb=False):
#recalcule les bornes max pour chaque valeur absolue
if Cd is None:
Cd = Set([c-c2 for c in self.C for c2 in self.C])
# if Cd is None:
# Cd = Set([c-c2 for c in self.C for c2 in self.C])
Cd = list(Cd)
if verb: print "Cd = %s"%Cd
m = dict([])
Expand Down Expand Up @@ -1690,6 +1697,58 @@ class BetaAdicMonoid(Monoid_class):
return ssd
return p

def intersection2 (self, ss1, ss2, verb = False): #calcule le sous-shift correspondant à l'intersection des deux monoïdes avec sous-shifts, utilise des FastAutomaton
r"""
Compute the intersection of two beta-adic monoid with subshifts given by FastAutomaton
INPUT:
- ``ss``- FastAutomaton (default: ``None``)
The first subshift to associate to the beta-adic monoid for this operation.
- ``ss2``- FastAutomaton (default: ``None``)
The second subshift to associate to the beta-adic monoid for this operation.
- ``verb``- bool (default: ``False``)
If True, print informations for debugging.
OUTPUT:
A FastAutomaton.
EXAMPLES:
#. Compute the boundary of the dragon fractal (see intersection_words for a easier way) ::
sage: m = BetaAdicMonoid(1/(1+I), {0,1})
sage: import sage.combinat.words.cautomata
sage: from sage.combinat.words.cautomata import FastAutomaton
sage: ss0 = FastAutomaton([(0,1,0)]+[(1,1,l) for l in m.C], i=0, final_states=[1])
sage: ss1 = FastAutomaton([(0,1,1)]+[(1,1,l) for l in m.C], i=0, final_states=[1])
sage: ssi = m.intersection2(ss0, ss1)
sage: m.plot2(tss = ssi) # long time
"""
a = self.relations_automaton3()
a = a.emonde_inf()
a.set_final_states(a.states())
ssp = ss1.product(ss2)
ssp = ssp.emonde()
d = {}
for (la1,la2) in ssp.Alphabet():
for lb in a.Alphabet():
if lb == la1-la2:
d[((la1,la2), lb)] = (la1,la2)
ssi = ssp.product(a, d)
ssi = ssi.emonde_inf()
ssi = ssi.emonde()
d = {}
for (l1,l2) in ssi.Alphabet():
d[(l1,l2)] = l1
ssi = ssi.determinise_proj(d)
ssi = ssi.emonde_inf()
ssi = ssi.emonde()
return ssi.minimise()

def intersection_words (self, w1, w2, ss=None, iss=None):
r"""
Compute the intersection of two beta-adic monoid with subshifts corresponding to two prefix
Expand Down
42 changes: 36 additions & 6 deletions src/sage/monoids/draw.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -378,7 +378,10 @@ void Draw (BetaAdic b, Surface s, int n, int ajust, Color col, int verb)
printf("maj = %lf\n", (1.-sqrt(cnorm(b.b))));
printf("abs(b) = %lf\n", sqrt(cnorm(b.b)));
}
n = -2.*log(max(s.sx, s.sy)*absd(1.-sqrt(cnorm(b.b))))/log(cnorm(b.b));
if (cnorm(b.b) < 1)
n = -2.*log(max(s.sx, s.sy)*absd(1.-sqrt(cnorm(b.b))))/log(cnorm(b.b));
else
n = 2.*log(max(s.sx, s.sy)*absd(1.-1./sqrt(cnorm(b.b))))/log(cnorm(b.b));
if (verb)
printf("n = %d\n", n);
}
Expand Down Expand Up @@ -407,7 +410,22 @@ void Draw (BetaAdic b, Surface s, int n, int ajust, Color col, int verb)
Fill(s, color0);
if (auto_n)
{
n = .75 - 2.*log(max(3.*s.sx/(Mx-mx), 3.*s.sy/(My-my)))/log(cnorm(b.b));
if (cnorm(b.b) < 1)
n = .75 + 2.*absd(log(max(3.*s.sx/(Mx-mx), 3.*s.sy/(My-my)))/log(cnorm(b.b)));
else
{
double ratio = pow(cnorm(b.b), n/2);
n = .75 + 2.*absd(log(ratio*max(3.*s.sx/(Mx-mx), 3.*s.sy/(My-my)))/log(cnorm(b.b)));
ratio = pow(cnorm(b.b), n/2)/ratio;
Mx = Mx*ratio;
mx = mx*ratio;
My = My*ratio;
my = my*ratio;
if (verb)
{
printf("Zone de dessin : (%lf, %lf) (%lf, %lf)\n", mx, my, Mx, My);
}
}
if (verb)
printf("n = %d\n", n);
}
Expand Down Expand Up @@ -468,7 +486,10 @@ void Draw2 (BetaAdic b, Surface s, int n, int ajust, Color col, int verb)
{
//printf("max = %d\n", max(s.sx, s.sy));
//printf("maj = %lf\n", (1.-sqrt(cnorm(b.b))));
n = -2.*log(max(s.sx, s.sy)*(1.-sqrt(cnorm(b.b))))/log(cnorm(b.b));
if (cnorm(b.b) < 1)
n = -2.*log(max(s.sx, s.sy)*(1.-sqrt(cnorm(b.b))))/log(cnorm(b.b));
else
n = 2.*log(max(s.sx, s.sy)*(1.-1./sqrt(cnorm(b.b))))/log(cnorm(b.b));
if (verb)
printf("n = %d\n", n);
}
Expand Down Expand Up @@ -497,7 +518,10 @@ void Draw2 (BetaAdic b, Surface s, int n, int ajust, Color col, int verb)
Fill(s, color0);
if (auto_n)
{
n = .75 - 2.*log(max(3.*s.sx/(Mx-mx), 3.*s.sy/(My-my)))/log(cnorm(b.b));
if (cnorm(b.b) < 1)
n = .75 + 2.*absd(log(max(3.*s.sx/(Mx-mx), 3.*s.sy/(My-my)))/log(cnorm(b.b)));
else
n = .75 + 2.*absd(log(max(.3*s.sx/(Mx-mx), .3*s.sy/(My-my)))/log(cnorm(b.b)));
if (verb)
printf("n = %d\n", n);
}
Expand Down Expand Up @@ -559,7 +583,10 @@ void DrawList (BetaAdic2 b, Surface s, int n, int ajust, ColorList cl, double al
{
//printf("max = %d\n", max(s.sx, s.sy));
//printf("maj = %lf\n", (1.-sqrt(cnorm(b.b))));
n = -2.*log(max(s.sx, s.sy)*(1.-sqrt(cnorm(b.b))))/log(cnorm(b.b));
if (cnorm(b.b) < 1)
n = -2.*log(max(s.sx, s.sy)*(1.-sqrt(cnorm(b.b))))/log(cnorm(b.b));
else
n = 2.*log(max(s.sx, s.sy)*(1.-1./sqrt(cnorm(b.b))))/log(cnorm(b.b));
if (verb)
printf("n = %d\n", n);
}
Expand Down Expand Up @@ -593,7 +620,10 @@ void DrawList (BetaAdic2 b, Surface s, int n, int ajust, ColorList cl, double al
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 (cnorm(b.b) < 1)
n = .5 + 2.*absd(log(max(3.*s.sx/(Mx-mx), 3.*s.sy/(My-my)))/log(cnorm(b.b)));
else
n = .5 + 2.*absd(log(max(.3*s.sx/(Mx-mx), .3*s.sy/(My-my)))/log(cnorm(b.b)));
if (verb)
printf("n = %d\n", n);
}
Expand Down

0 comments on commit c18febb

Please sign in to comment.