#############################################################################
##
#W  foldings.gi      GAP library     Manuel Delgado <[email protected]>
#W                                   Jose Morais    <[email protected]>
##
#H  @(#)$Id: foldings.gi,v 1.13 $
##
#Y  Copyright (C)  2004,  CMUP, Universidade do Porto, Portugal
##
#############################################################################
## A finitely generated subgroup of a free group of finite rank rk can be given
## as a list [rk, gen1, gen2,...]. The generators can be given as strings
## on the generators of the free group (and its inverses which are 
## represented by the corresponding capital letters) or as lists of integers
## where if i<=rk then  i represents the ith generator; if i>rk, then i
## represents the inverse of the rk-ith generator. The generators of the 
## free group are assumed to be a, b, c, ...
##
## Example: [2,"abA","bbabAB"] means the subgroup of the free group on 2 
## generators generated by aba^{-1} ...
##
## Another representation could be [2,[1,2,3],[2,2,1,2,3,4]].
##
#############################################################################
##
#F IsGenRep(L)
##
##
InstallGlobalFunction(IsGenRep, function(L)
    local   abc,  ABC,  alph;
    
    if IsPosInt(L[1]) then 
        abc := "abcdefg";
        ABC := "ABCDEFG";
        alph := Concatenation(abc{[1..L[1]]},ABC{[1..L[1]]}); 
    else 
        return false;
    fi;
    
    return ForAll(L{[2..Length(L)]}, x-> IsString(x) and IsSubset(alph,x));
end);
#############################################################################
##
#F IsListRep(L)
##
##
InstallGlobalFunction(IsListRep, function(L)
    return IsPosInt(L[1])  and ForAll(L{[2..Length(L)]}, x-> IsList(x) and ForAll(x, y->IsPosInt(y) and y <= 2 * L[1]));
end);
##
#############################################################################
## The following functions allow us to pass from one representation to 
## another
##
#############################################################################
##
#F GeneratorsToListRepresentation(L)
##
## L is a list whose first element is the number of generators of the 
## free group. The remaining elements are the generators of the subgroup. 
##
## Example: when the input is [2,"abA","bbabAB"], the output will be
## [2,[1,2,3],[2,2,1,2,3,4]]
##
## Warning: Alphabets with more than 7 letters must not be used 
##
InstallGlobalFunction(GeneratorsToListRepresentation, function(L)
    local   K,  abc,  ABC,  alph,  g,  T;
    
    if not IsPosInt(L[1]) or L[1] > 7 then
        Error("The rank in IsGeneratorsToListRepresentation must be as an integer not greater that 7");
    fi;
    
    if not IsGenRep(L) then  
        Error("The generators in GeneratorsToListRepresentation must be given as strings");
    fi;
    
    K := [L[1]];
    
    abc := "abcdefg";
    ABC := "ABCDEFG";
    alph := Concatenation(abc{[1..L[1]]},ABC{[1..L[1]]});        
    
    for g in [2.. Length(L)] do
        T := List(L[g], i -> Position(alph,i));
        Add(K,T);
    od;
    return K;
end);
#############################################################################
#F  ListToGeneratorsRepresentation(K) 
##
## is the inverse of GeneratorsToListRepresentation
##
InstallGlobalFunction(ListToGeneratorsRepresentation, function(K)
    local   L,  abc,  ABC,  alph,  g;
    
    if not IsPosInt(K[1]) or K[1] > 7 then
        Error("The rank in IsListToGeneratorsRepresentation must be as an integer not greater that 7");
    fi;
    
    if not IsListRep(K) then  
        Error("The generators in ToListGeneratorsRepresentation must be given as lists of integers");
    fi;
    
    L := [K[1]];
    
    abc := "abcdefg";
    ABC := "ABCDEFG";
    alph := Concatenation(abc{[1..L[1]]},ABC{[1..L[1]]});        
    
    for g in [2.. Length(K)] do
        Add(L,alph{K[g]});
    od;
    return L;
    
end);
#############################################################################
##
#F FlowerAutomaton(L)
##
## Given a finitely generated subgroup of a free group (by any of the two 
## means indicated above) the flower automaton is constructed.
##
InstallGlobalFunction(FlowerAutomaton, function(L)
    local   n,  abc,  ABC,  alph,  states,  i,  q,  T,  j,  g,  p,  a;
    
    if IsListRep(L) then
        L := ListToGeneratorsRepresentation(L);
    elif not IsGenRep(L) then
        Error("The argument of FlowerAutomaton must be a representation of a subgroup of the free group");
    fi;
    
    n := L[1];
    
    abc := "abcdefg";
    ABC := "ABCDEFG";
    alph := Concatenation(abc{[1..L[1]]},ABC{[1..L[1]]});        

    states := 1;
    for i in [2..Length(L)] do
        states := states + Length(L[i]) - 1;
    od;
    q := 1;
    T := NullMat(n,states);
    for i in [1..n] do
        for j in [1..states] do
            T[i][j] := [];
        od;
    od;
    for i in [2..Length(L)] do
        g := L[i];
        if Length(g) > 1 then
            q := q+1;                #add a new state
            p := Position(alph, g[1]);
            if p <= n then
                AddSet(T[p][1],q);
            else
                p := p - n;
                AddSet(T[p][q], 1);
            fi;
        else
            p := Position(alph, g[1]);
            if p <= n then
                AddSet(T[p][1],1);
            else
                p := p - n;
                AddSet(T[p][1], 1);
            fi;
        fi;
        
        for a in [2..Length(g)-1] do 
            q := q+1;
            p := Position(alph, g[a]);
            if p <= n then
                AddSet(T[p][q-1],q);
            else
                p := p - n;
                AddSet(T[p][q], q-1);
            fi;
        od;
        p := Position(alph, g[Length(g)]);
        if p <= n then
            AddSet(T[p][q],1);
        else
            p := p - n;
            AddSet(T[p][1], q);
        fi;
    od;
    
    return Automaton("nondet",states, n, T, [1],[1]); #n is to be replaced by alph
end);

#############################################################################
##
#F FoldFlowerAutomaton(arg)
##
## The first (and usually also the last) argument must be a flower automaton.
## (The first state must be the initial and final state; all vertices, except 
## the initial state, must be of degree 2.)
##
## The second argument, when present, only has effect when it is <true>. 
## WARNING: It should only be used when facilities to draw automata are 
## avaiable. In that case, one may visualize the identifications that 
## are taking place.
## 
## Makes Stallings foldings on the flower automaton <A>
##
##
InstallGlobalFunction(FoldFlowerAutomaton, function(arg)
    local   bool,  A,  ug,  n,  na,  ns,  T,  changes1,  changes2,  identify,  
            deleteAndRename,  a,  q,  p,  c1,  c2,  c,  newtable,  b,  aut,  
            s,  r;
    
    bool := false;
    A := arg[1];
    if IsBound(arg[2]) and arg[2] = true then
        bool := true;
    fi;
    
    
    if not A!.type = "nondet" then
        Error(" A must be non deterministic");
    fi;
    if not (A!.initial = [1] and A!.accepting = [1]) then
        Error(" 1 must be initial and accepting state");
    fi;
    ug := UnderlyingMultiGraphOfAutomaton(A);
    if not ForAll([2..A!.states], q -> AutoVertexDegree(ug,q)=2) then
        Error(" A must be a flower automaton");
    fi;
        
    n := 1;
    
    na := A!.alphabet;
    ns := A!.states;
    T := StructuralCopy(A!.transitions);
    changes1 := true;
    changes2 := true;
    
    
    ####################################
    identify := function(p1,p2)
        local   a,  q;
        
        if p2 = 1 then  # let the initial state never be removed
            p2 := p1;
            p1 := 1;
        fi;
        if bool then
            Print("I am identifying states ",p1, " and ",p2, "\n");
        fi;
        
        for a in [1..na] do # all occurrences of p2 in the transition 
                             # matrix are substituted by p1.
            for q in [1..ns] do
                if p2 in T[a][q] then
                    T[a][q] := Union(T[a][q],[p1]);
                    T[a][q] := Difference(T[a][q],[p2]);
                fi;
            od;
        od;
        
        for a in [1..na] do 
            T[a][p1] := Union(T[a][p1], T[a][p2]);
            if not p1 = p2 then
                 T[a][p2] := [];
            fi;
        od;
        
        T[a][p1] := Set(Flat(T[a][p1]));
        SubtractSet(T[a][p1],[0]);
        if bool then
            n := n+1;
            DrawAutomaton(Automaton("nondet",ns,na,T,[1],[1]),String(n));
#            Error("...");
        fi;
        
    end;
    ######################
    deleteAndRename := function(T,c)# delete a list c of vertices
        local   TR,  acc,  nt,  newtable,  n1,  n2,  newnewtable,  r,  s;
            
        TR := TransposedMat(T);
        acc := Difference([1..Length(T[1])],c);
    
        nt := TR{acc};
        newtable := TransposedMat(nt);
        
        n1 := Length(newtable);
        n2 := Length(newtable[1]);
        newnewtable := NullMat(n1,n2);
        for r in [1 .. n1] do
            for s in [1 .. n2] do
                if newtable[r][s] <> 0 then
                    if Position(acc, newtable[r][s]) <> fail then
                        newnewtable[r][s] := Position(acc, newtable[r][s]);
                    fi;
                else
                    newnewtable[r][s] := 0;
                fi;
                
            od;
        od;
        return newnewtable;
    end;
    ###########################
    while changes1 or changes2 do 
       while changes1 do
            changes1 := false;
            for a in [1..na] do
                for q in [1..ns] do
                    if Length(T[a][q]) > 1 then
                        changes1 := true;
                        changes2 := true;
                        identify(T[a][q][1],T[a][q][2]);
                    fi;
                od;
            od;
        od;
        while changes2 do
            changes2 := false;
            for a in [1..na] do
                for p in [1..ns] do
                    for q in [1..ns] do
                        if p <> q and Intersection(T[a][p],T[a][q]) <> [] then
                            changes1 := true;
                            changes2 := true;
                            identify(p,q);
                        fi;
                    od;
                od;
            od;
        od;
    od;
    for a in [1..na] do
        for q in [1..ns] do
            if T[a][q] <> [] then
                T[a][q] := T[a][q][1];
            else
                T[a][q] := 0;
            fi;
        od;
    od;
    ### computes the inaccessible states
    c1 := Filtered([1..ns], q -> ForAll([1..na],a -> T[a][q] = 0));
    c2 := Difference([1..ns],Set(Flat(T)));
    c := Intersection(c2,c1);
    
    newtable := deleteAndRename(T,c); ## removes the inaccessible states
    
    ## remove states of degree 1
    b := true;
    while b do
        b := false;
        aut := Automaton("det", Length(newtable[1]), na, newtable,[1],[1]);
        ug := UnderlyingMultiGraphOfAutomaton(aut);
#        ug := UnderlyingGraphOfAutomaton(aut);
        T := aut!.transitions;
        s := []; #list of vertices of degree 1
        for r in [2..aut!.states] do
            if AutoVertexDegree(ug,r) = 1 then
                Add(s,r);
            fi;
        od;
        if s <> [] then
            b := true;
            newtable := deleteAndRename(T,s); ## removes states of degree 1
        fi;
    od;
   
    aut := Automaton("det", Length(newtable[1]), na, newtable,[1],[1]);
    if bool then
        DrawAutomaton(aut,"aut");
    fi;
    return aut;

end);

#############################################################################
##
#F SubgroupGenToInvAut(L)
##
## Returns the inverse automaton corresponding to the subgroup given by 
## <A>L</A>.
InstallGlobalFunction(SubgroupGenToInvAut,function(L)
    return FoldFlowerAutomaton(FlowerAutomaton(L));
end);

##########################################################################
##
#F AddInverseEdgesToInverseAutomaton(aut)
##
## Given an inverse automaton, adds the edges labeled by the inverses
##
InstallGlobalFunction(AddInverseEdgesToInverseAutomaton,function(aut)
    local   T,  q,  L,  i,  a,  ai,  alph;

    if not IsInverseAutomaton(aut) then
        Error("The argument must be an inverse automaton");
    fi;
    if not IsInt(AlphabetOfAutomaton(aut)) then
        Error("The automaton must be defined over the alphabet abc...");
    fi;
    T := StructuralCopy(aut!.transitions); 
    q := aut!.states;
    for L in T do 
        for i in [1..Length(L)] do
            if IsBound(L[i]) and L[i] = 0 then
                Unbind(L[i]);
            fi;
        od;
    od;
    for a in aut!.transitions do
        ai := [];
        for i in [1..q] do
            if i in a then
                Add(ai, Position(a, i));
            else
                Add(ai,0);
            fi;
        od;
        Append(T,[ai]);
    od;
    for L in T do 
        for i in [1..q] do
            if not IsBound(L[i]) then
                L[i] := 0;
            fi;
        od;
    od;
    alph := "";
    for i in [1 .. aut!.alphabet] do
        alph := Concatenation(alph, [jascii[68+i]]);
    od;
    for i in [1 .. aut!.alphabet] do
        alph := Concatenation(alph, [jascii[68+i-32]]);
    od;
    FamilyObj(aut)!.alphabet := alph;
    aut!.alphabet := Length(alph);
    aut!.transitions := T;
#    return(aut);
#    return Automaton(aut!.type,aut!.states,alph,T,aut!.initial,aut!.accepting);
end);

#############################################################################
##
#F GeodesicTreeOfInverseAutomatonWithInformation
##
## Is an auxiliar function to the following functions
## InverseAutomatonToGenerators and GeodesicTreeOfInverseAutomaton
##
InstallGlobalFunction(GeodesicTreeOfInverseAutomatonWithInformation, function(A)
    local   Ainv,  T,  visited,  bool,  NEW,  lista,  u,  new,  a,  ai,  
            tree;
    
    if not IsInverseAutomaton(A) or A!.accepting <> A!.initial 
       or Length(A!.initial) <> 1 then
        Error("<A> must be an inverse automaton");
    fi;
    Ainv := Automaton(A!.type,A!.states,A!.alphabet,StructuralCopy(A!.transitions),A!.initial,A!.accepting);
    AddInverseEdgesToInverseAutomaton(Ainv);
    T := StructuralCopy(Ainv!.transitions);
    visited := [Ainv!.initial[1]];
    bool := true;
    NEW := [Ainv!.initial[1]];
    lista := [];
    for u in [1..A!.states] do
        Add(lista,[]);
    od;
    
    while bool do
        new := ShallowCopy(NEW);
        NEW := [];
        bool := false;
        for u in new do 
            for a in [1..Ainv!.alphabet] do
                if not IsBound(T[a][u]) or T[a][u] in visited or T[a][u] = 0 then
                    T[a][u] := 0;
                else
                    bool := true;
                    Add(visited, T[a][u]);
                    Add(NEW, T[a][u]);
                    lista[T[a][u]] := Concatenation(lista[u],[a]);
                fi;
            od;
        od;
    od;
    for ai in [A!.alphabet+1..Ainv!.alphabet] do
        a := ai - A!.alphabet;
        for u in [1..A!.states]do
            if T[ai][u] <> 0 then
                T[a][T[ai][u]] := u;
            fi;
        od;
    od;
    T := T{[1..A!.alphabet]};
    
    tree := Automaton(A!.type,A!.states,A!.alphabet,T,A!.initial,A!.accepting);

    return [tree,lista];
end);

#############################################################################
##
#F GeodesicTreeOfInverseAutomaton
##
## Returns an automaton whose underlying graph is a geodesic tree of the 
## underlying graph of the automaton given.
##
InstallGlobalFunction(GeodesicTreeOfInverseAutomaton, function(A)
    return GeodesicTreeOfInverseAutomatonWithInformation(A)[1];
end);

#############################################################################
##
#F InverseAutomatonToGenerators
##
## returns a set of generators (given trough the representation above) of the 
## subgroup of the free group corresponding to the automaton given. 
##
InstallGlobalFunction(InverseAutomatonToGenerators, function(A)
    local a, ll, i, GEN, gen, generator, g, e, abc, ABC, alph, PO, T, tree, 
          lista, u, 
          posedges; #positive edges that are not part of the geodesic tree.
    
    if A!.alphabet > 7 then
        Error("The alphabet in GeodesicTreeOfInverseAutomaton must be given as an integer not greater that 7");
    fi;
    abc := "abcdefg";
    ABC := "ABCDEFG";
    alph := Concatenation(abc{[1..A!.alphabet]},ABC{[1..A!.alphabet]});        
    
    lista := GeodesicTreeOfInverseAutomatonWithInformation(A)[2];
    T := GeodesicTreeOfInverseAutomatonWithInformation(A)[1]!.transitions; 
    
    PO := StructuralCopy(A!.transitions);
    
    for a in [1..A!.alphabet] do
        for u in [1..A!.states]do
            if T[a][u] <> 0 then
                PO[a][u] := 0;
            fi;
        od;
    od;    
    posedges := [];
    for a in [1..A!.alphabet] do
        for u in [1..A!.states]do
            if PO[a][u] <> 0 then
                Add(posedges, [u,a,PO[a][u]]);
            fi;
        od;
    od;    
    gen := [];
    for e in posedges do
        generator := ShallowCopy(lista[e[1]]);
        Add(generator, e[2]);
        
        ll := List(lista[e[3]], i -> (i + A!.alphabet) mod (2* A!.alphabet));
        for i in [1..Length(ll)] do
            if ll[i] = 0 then
                ll[i] := 2* A!.alphabet;
            fi;
        od;
        
        generator := Concatenation(generator, Reversed(ll));
        Add(gen, generator);
    od;
    GEN := [];
    for g in gen do
        Add(GEN, alph{g});
    od;
    return Concatenation([A!.alphabet],GEN);
end);


#E