# This code was written by Sara Billey at the University of
#Washington using Maple versions 8 and 9.  The package contains
#code to solve Schubert problems and to compute the unique permutation
#array associated to a list of permutations.  The algorithms are based on
#the paper "Intersections of Schubert varieties  and other permutation array schemes"
#by Sara Billey and Ravi Vakil (math.AG/0502468).  No support or guarantee is promised
#with this code.  Use at your own risk.

#### Start up

restart;
with(LinearAlgebra);
with(combinat):
with(ListTools,Reverse);


####Vector/Matrix tools

ScaleVector := proc (v)  ##forces first non-zero entry to be 1
local w;
w:=remove((x)->x=0,v);
    if nops(w)=0 then v else simplify(expand((1/(w[1])) * v)); fi;
end: # ScaleVector

Myrank := proc(m)
local b;
    b:=nops(Basis([seq(Column(m,i), i=1..ColumnDimension(m))]));
end: # Myrank


#### Permutation Array Tools
shadow := proc(pt, arr)
description `find all points in the principle array defined by pt in the array arr`;
select((x)->dominantorder(x,pt), arr);
end: # shadow


dominantorder := proc(a,b)
description `checks a in principle subarray of b`;
local i,k;
    (nops(a)=nops(b)) and (0=nops(select((i)->b[i]<a[i], [seq(k,k=1..nops(a))])));
end: # dominantorder


lexorderlists := proc(a,b)
    nops(a)=0
    or (not(nops(b)=0)
        and (a[1] < b[1]
             or (a[1] = b[1] and lexorderlists(a[2..nops(a)],b[2..nops(b)]))));
end: # lexorderlists

revlexorderlists := proc(a,b)
    nops(a)=0
    or (not(nops(b)=0)
        and (a[-1] < b[-1]
             or (a[-1] = b[-1] and revlexorderlists(a[1..-2],b[1..-2]))));
end: # revlexorderlists


flagpositions := proc(array)
boardperms(permarraytoboards(array), permrank(array), nops(array[1]) -1);
end: # flagpositions

boardperms := proc(boards,n,d)
local A,w,index,i,j,k;
    for k from 1 to d do
        A:=seq(map((x)->x[k],boards[j]), j=1..n);
        w[k]:=[seq(op({op(A[n-j])} minus {op(A[n-j-1])}), j=0..n-2),op(A[1])];
    od;
    return(convert(w,list));
end: # boardperms

dotoccupiedpositions := proc(arr,index)
   [op({op(map((x)->x[index], arr))})]
end: # dotoccupiedpositions


makeantichainsfromarray := proc(arr,n,d)
local i,x,y;
    [seq( map((y)->y[1..d-1],  select((x)->x[d]=i, arr)), i =1..n)]
end: # make_antichains_from_array

dotdownsize := proc(antichain, arr)
description `Remove antichain, include covered points, check rank condition`;
local B,coveredpts;
    if not(nops(antichain[1]) = nops(arr[1])) then print("Error in dotdownsize"); fi;
    coveredpts:={op(map(dotjoin, [seq(op(select((a)->(2=permrank(shadow(dotjoin(a),arr)))  ## assuming only rank 2 covered points will contribute.  All others are redundant.
                                                and  coveringP(a),
                                                choose(antichain,i))),i=2..3)]))};
    B:=[op(({op(arr)} minus {op(antichain)}) union coveredpts)];
    return(B);
end: # dotdownsize

permarraytoboards:= proc(arr)
description `from 4-dim array list 3-dim boards corresponding to intersection with F_i`;
local B,A,i,count,n,x,table, helper;
    n:=permrank(arr);
    A:= makeantichainsfromarray(arr,n,4);
    helper := proc(x)
        if not(assigned(table[x[1],x[2]])) or table[x[1],x[2]] > x[3]
        then table[x[1],x[2]]:=x[3]; fi;
    end: # helper
    map(helper, arr);
    table:=op(op(table));
    B[n]:= map((x)->[op(x)], table);
    for i from 0 to n-1 do
        count:=n-i;
        B[count-1]:=dotdownsize(A[count],B[count]);
    od;
    return(B);
end: # permarraytoboards

permrank := proc(arr)  # assumes arr is totally rankable
    nops({op(map((x)->x[1], arr))});
end: # permrank

DotRank := proc(arr,j)  # does not assume arr is totally rankable
    nops({op(map((x)->x[j], arr))});
end: # dotrank

TotallyRankableP := proc(arr,n,d)
local helper,result;
    result:=true;
    helper := proc(x)
        if not(RankableP(shadow(x,arr),x)) then result:=false; return(result); fi;
    end: # helper
    AllCoordinateTester(d,[seq(1,j=1..d)],[seq(n,j=1..d)],[],helper);
    return(result);
end: # TotallyRankableP

RankableP := proc(arr,x)
local s;
    s:={seq(DotRank(arr,j), j=1..nops(x))};
    evalb(nops(s)=1);
end: # RankableP


coveringP := proc (alist)
description `check if pts in alist form a covering set for their join`;
local x, i;
    x:=dotjoin(alist);
    not(member(x,alist))
    and nops(alist) = nops(select((y)->member(0,x-y), alist))
    and not(member(0,x-dotmeet(alist)))
end: # coveringP

dotmeet := proc (alist)
local j;
[seq(min(op(map((y)->y[j],alist))), j=1..nops(alist[1]))];
end: # dotmeet

dotjoin := proc (alist)
local j;
[seq(max(op(map((y)->y[j],alist))), j=1..nops(alist[1]))];
end: # dotmeet

printdotarray := proc (arr,n,d)
local A,B,i,j,x,k;
    if 3=d
    then A:=array(1..n,1..n,[seq([seq([],i=1..n)],j=1..n)]);
        for x in arr do A[x[1],x[2]]:=x[3]; od;
        print(A);
    fi;
    if 4=d
    then B:=  permarraytoboards(arr);
        for k from 1 to n do
            A[k]:=array(1..n,1..n,[seq([seq([],i=1..n)],j=1..n)]);
            for x in B[k] do A[k][x[1],x[2]]:=x[3]; od;
        od;
        print(seq(A[k],k=1..n));
    fi;
end: # printdotarray

print3dotarraylist := proc (arrlist,n)
local A,s,i,j,x,k,step;
    for step from 1 to binomial(n, 2)+1 do
        print (Step, step-1);

        for s from 1 to nops(arrlist)
        do  A[s]:=array(1..n,1..n,[seq([seq([],i=1..n)],j=1..n)]);
            for x in arrlist[s][step] do A[s][x[1],x[2]]:=x[3]; od;
        od;
        if 1 = nops({seq(A[i], i=1..nops(arrlist))})
        then          print(A[1]);
        else print(seq(A[i], i=1..nops(arrlist)));
        fi;
    od;
end: #


printvectortable := proc(vectors,n)
local eq,x,y,A,v;
    print(op(vectors));
    A:=array(1..n,1..n,[seq([seq([],i=1..n)],j=1..n)]);
    for eq in [op(op(op(vectors)))]do x:=op(eq)[1]; v:=eval(op(eq)[2]+Vector([seq(0,i=1..n)]));  A[x[1],x[2]]:= v; od;
    print("vectors:",A);
    #latex(A);
end: # printvectortable

printvectors:= proc(bb,vectors,n)   ## bb is a perm array
local solvectors,x;
    for x in bb do solvectors[x]:=eval(vectors[x]); od;
    printvectortable(solvectors,n);
end:


printforlisp := proc(alist)
description `convert maple format to lisp`;
map((x)->x+[seq(-1, i=1..nops(x))], alist);
end: # printforlisp


ListAllCoordinates := proc(minvals,maxvals)
local results,helper;
global COORDRESULTS;
    helper := proc(tail)
    global COORDRESULTS;
        COORDRESULTS:=COORDRESULTS,tail;
    end: # helper
    COORDRESULTS:=NULL;
    AllCoordinateTester(nops(minvals),minvals,maxvals,[],helper);
COORDRESULTS:=[COORDRESULTS];
return(COORDRESULTS);
end: # ListAllCoordinates note, I am using this ordering in ranktable


AllCoordinateTester := proc(j,minvals,maxvals,tail,helper)
local a;
    if j=0 then apply(helper,tail);
    else for a from minvals[j] to maxvals[j]
         do AllCoordinateTester(j-1,minvals, maxvals,  [a,op(tail)], helper);
         od;
    fi;
end: # AllCoordinateTester, note, I am using this ordering in ranktable

ranktable := proc(perms)
description `recursively compute the unique perm array for given perms`;
local table,i,j,l,n,d,mat,xi,index;
global RANKTABLEVAR;
    n:=nops(perms[1]);
    d:=nops(perms)+1;
    RANKTABLEVAR:=table();
    RANKTABLEVAR[1]:=n;
    ### establish the base case
    for i from 1 to nops(perms) do
        perms[i];
        mat:=permmatrix(Reverse(perms[i]));
        for xi from 1 to n do
            for j from 1 to n do
                index:=subsop(i=xi,[seq(n,i=1..d-1),j]);
                RANKTABLEVAR[index]:=rank(submatrix(mat, [$1..xi], [$1..j]));
            od;    od;  od;
    ### recursively set all other values
    AllCoordinateTester(d,[seq(1,j=1..d)],[seq(n,j=1..d)],[],recranktable);
    RANKTABLEVAR;
end: # permarray


recranktable := proc(index)
local x;
global RANKTABLEVAR;
    if assigned(RANKTABLEVAR[index])
    then return(RANKTABLEVAR[index]);
    elif member(0,index) then return(0);
    else x:=recranktable(takefirst(index)) + recranktable(dropfirst(index)) - index[-1];
        RANKTABLEVAR[index]:=max(0,x,recranktable(subsop(-1=index[-1]-1,index)));
        return(RANKTABLEVAR[index]);
    fi;
end: # recranktable

takefirst := proc (alist)
description `keep list up to first entry not equal to n and append on all n's after that.`;
    local n,i,count;
    count:=0;
    n:=RANKTABLEVAR[1];
    for i from 1 to n
    while alist[i]=n do count:=count+1; od;
    if count>nops(alist)-1 then alist;
    else [seq(n, i=1..count),alist[count+1],seq(n,i=count+2..nops(alist)-1),alist[-1]]; fi;
end: # takefirst

dropfirst := proc (alist)
description `keep list except first entry not equal to n and append on all n's after that.`;
    local n,i,count;
    count:=0;
    n:=RANKTABLEVAR[1];
    for i from 1 to n
    while alist[i]=n do count:=count+1; od;
    if count>nops(alist)-1 then alist;
    else [seq(n, i=1..count),n,op(alist[count+2..nops(alist)])]; fi;
end: # dropfirst

printranktable := proc(table)
local i,j,k,l,n;
    n:=table[1];
    for k from 1 to n do
        printf("\n");
        for i from 1 to n do
            printf("\n");
            for l from 1 to n do
                printf("    ");
                for j from 1 to n do
                    printf("%a", table[[i,j,k,l]]);
                od;
            od;
        od;
        printf("\n");
    od;
end: # printranktable

permarray := proc(permlist)
description `convert the recursively constructed ranktable to a permutation array`;
local result,helper,table,n,d;
    table:=ranktable(permlist);
    n:=nops(permlist[1]);
    d:=nops(permlist)+1;
    result:=NULL;
    helper := proc(tail)
        if not(table[tail]=0)
        and table[tail]>max(seq(table[subsop(k=tail[k]-1,tail)],
                           k=select((i)->not(tail[i]=1), [$1..d])))
        and not(redundantP(tail,[result]))
        then result:=result,tail; fi;
    end: # helper
    AllCoordinateTester(d,[seq(1,j=1..d)],[seq(n,j=1..d)],[],helper);
    return([result]);
end: # permarray



redundantP := proc(x,arr)
local shad;
    shad:=select((y)->not(y=x), shadow(x,arr));
    not(nops(shad)=0) and permrank(shad)=permrank([op(shad),x]);
end: # redundantP


TripleIdentityArrayP := proc(arr)
description `Check arr has 3 flags in transverse position and top board is the generic triple array`;
local B,A,i,count,n,x,table, helper;
    n:=permrank(arr);
    helper := proc(x)
        if not(assigned(table[x[1],x[2]])) or table[x[1],x[2]] > x[3]
        then table[x[1],x[2]]:=x[3]; fi;
    end: # helper
    map(helper, arr);
    table:=op(op(table));
    B[n]:= sort(map((x)->[op(x)], table),lexorderlists);
    evalb(B[n] = TripleIdentityArray(n,3));
end: # TripleIdenityArrayP

TripleIdentityArray := proc(n)
description `top board for 3 transverse flags`;
    sort([seq(seq( [(n - c),( (n-k+1)+ c), k],
             c=0..k-1),
        k=1..n)],  lexorderlists);
end: # TripleIdenityArray

#### Equation Tools
RemoveTrivialEqns := proc(eqns)
remove((x)->op(x)[1] = op(x)[2], eqns);
end: # RemoveTrivialEqns

Minors := proc(M,r)
description `list of all rxr minors in Matrix M`;
    map((y)->op(map((x)->det(submatrix(M, x,y)),
                    choose([$1..RowDimension(M)],r))),
        choose([$1..ColumnDimension(M)],r));
end: # Minors

#### Permutation Tools
permmatrix := proc(perm)
local helper,i,j;
helper:=proc(k,l) if perm[k]=l then 1 else 0; fi; end:
array([seq([seq(helper(i,j), i=1..nops(perm))], j=1..nops(perm))]);
end: # permmatrix

PermsFixedLength := proc(k,n)
description `find all perms of length k in S_n`;
local results,helper;
    helper := proc(tail)
        results:=results,tail;
    end: # helper
    results:=NULL;
    AllPermFixedLengthTester(k,n,[],[$1..n],helper);
results:=[results];
return(results);
end: # permslength


AllPermFixedLengthTester := proc(k,n,tail,stack,helper)
local a;
    if k=0 then apply(helper,[op(tail),op(stack)]);
    else for a from max(0,k - binomial(nops(stack)-1, 2))
               to min((nops(stack)-1), k)
         do AllPermFixedLengthTester(k-a,n,[op(tail),stack[a+1]], subsop((a+1)=NULL,stack),helper);
         od;
    fi;
end: # AllCoordinateTester, note, I am using this ordering in ranktable

W0 := proc(n)
    [seq(n-i,i=0..n-1)];
end: # W0
MultPerms := proc(u,v)
description `multiply perms in 1-line notation`;
    map((x)->u[x], v);
end: # MultPerms

PermLength := proc(perm)
local count,j;
    count:=0;
    for j from 1 to nops(perm)
    do count:=count+nops(select((x)->perm[j]>x, perm[(j+1)..nops(perm)])); od;
    count;
end: # PermLength

##### Specific Flags

RandomFlag := proc(n)
description `Choose a random flag in R^n`;
local v;
    v:=randmatrix(n,n,unimodular);
    [seq(Vector(n,(x)->v[i,x]), i=1..n)];
end: # RandomFlag


FixedFlag := proc(n)
[seq(UnitVector(j,n), j=1..n)];
end: # FixedFlag

BinomialFlag := proc (n)
local vectors, pt;
    [seq(Vector([seq(binomial(n-j,k), k=0..n-j),seq(0,l=(n-j+2)..n)]), j=1..n)];
end: # binomialflag

TransverseFlags := proc(n,d)
description `Give a list of d flags hopefully in transverse position`;
    if d=0 then NULL;
    elif d=1 then [FixedFlag(n)];
    elif d=2 then [FixedFlag(n), Reverse(FixedFlag(n))];
    elif d>2 then [FixedFlag(n), Reverse(FixedFlag(n)), seq(RandomFlag(n), i=1..d-2)];
    fi;
end: # TransverseFlags


GenericFlag := proc (n)
local i,j;
    [seq(Vector([seq(0,j=1..i-1),1,seq(XX[i,j], j=i+1..n)]), i=1..n)];
end: # GenericFlag



RankToPermArray := proc(atable,n,d)
description `convert atable to a permutation array`;
local a, helper;
    a:=NULL;
    helper :=proc(tail)
        ## add a dot if necessary
        if not({atable[tail]} = {seq(DotRank(shadow(tail, [a]),j), j=1..d)})
        then
            a:=tail,a;
        fi;
        ## verify rankable at tail and correct rank
        if not(RankableP(shadow(tail,[a]), d))
        then print("dot rank Problem", tail);
        fi;
    end: #helper
    ###initial table setup
    AllCoordinateTester(d, [seq(1, r=1..d)], [seq(n, r=1..d)], [],helper);
    [a];
end: # RankToPermArray

FlagPermArray := proc (flaglist)
local n,d;
    d:=nops(flaglist);
    n:=Dimension(flaglist[1][1]);
    RankToPermArray(FlagRankTable(flaglist), n,d);
end: # FlagPermArray


GenericPermArray := proc(n,d)
    description `Compute vectors in 1-dim subspaces `;
local results,helper;
    results:=NULL;
    helper :=proc(tail)
        if add(j, j=tail) = (d-1)*n+1
        then results:=results,tail;
        fi;
    end: #helper
    ###initial table setup
    AllCoordinateTester(d, [seq(1, r=1..d)], [seq(n, r=1..d)], [],helper);
    [results];
end: # GenericPermArray


ProjArray := proc(P,i)
description `compute P_i from P`;
local Q;
Q:=select((x)->x[-1] <= i, P);
Q:=map((x)->x[1..-2], Q);
select((x)->not(redundantP(x,Q)), Q);
end: # ProjArray

AntichainArray := proc (P,i)
description `Find the antichain A_i in the EL algorithm`;
local d;
    d:=nops(P[1]);
    select((x)->x[d]=i, P);
end: # AntichainArray

##### Flag Tools


FlagIntersection := proc(E,x)
description `E=flaglist, x=coordinate in [n]^d. Give a basis for  E_x_1, ...E_x_d  `;
    IntersectionBasis([seq(E[i][1..x[i]], i=1..nops(x))]);
end: # FlagIntersection

FlagVectorTable := proc(flaglist)
description `Compute vectors in 1-dim subspaces `;
local results,helper, base, vtable,x,n,d,i,j,index,A,invA,k;
global FVTABLE, FRTABLE, PERMARRAY;
    unassign('FVTABLE');
    unassign('FRTABLE');
    PERMARRAY:=NULL;
    d:=nops(flaglist);
    n:=Dimension(flaglist[1][1]);
    vprint(n,d);
    helper :=proc(tail)
        #if add(j, j=tail) = (d-1)*n+1
        base:=FlagIntersection(flaglist, tail);
        if 1 = nops(base) then
            FVTABLE[tail]:= eval(base[1]);
            PERMARRAY:=PERMARRAY,tail;
        elif nops(base) > max(seq(DotRank([PERMARRAY], j), j=1..d))
        then PERMARRAY:=PERMARRAY,tail;
        fi;
    end: #helper
    ###initial table setup
    AllCoordinateTester(d, [seq(1, r=1..d)], [seq(n, r=1..d)], [],helper);
    PERMARRAY:= [PERMARRAY];
    if PERMARRAY = GenericPermArray(n,d)
    then vprint("generic");
    else print("not generic");
    fi;
    FVTABLE;
end: # FlagVectorTable

FlagRankTable := proc(flaglist)
description `Computes dimensions of all principle subspaces `;
local results,helper, base, vtable,x,n,d,i,j,index,A,invA,k;
global FRTABLE;
    unassign('FRTABLE');
    d:=nops(flaglist);
    n:=Dimension(flaglist[1][1]);
    helper :=proc(tail)
        FRTABLE[tail]:=nops(FlagIntersection(flaglist, tail));
    end: #helper
    AllCoordinateTester(d, [seq(1, r=1..d)], [seq(n, r=1..d)], [],helper);
    FRTABLE;
end: # FlagRankTable

######  Main Function: Solving For Flags in Generalized Schubert Cell

FlagSolver := proc(P,flaglist)
description `P=perm array, flaglist=list of flags in general position.`;
local vectors, boards, solution, V, pt,i,j,k,l,n,d, perms;
global FLAGSOLS, FLAGEQNS,FLAGVARIABLES;
    n:=permrank(P);
    d:=nops(P[1]);
    perms:=flagpositions(P);
    print("Testing:", perms);
    V:=FlagVectorTable(flaglist);
    #vectors:=eval(V);  #copy V
    FLAGSOLS:=[];
    FLAGEQNS:=NULL;
    FLAGVARIABLES:=NULL;
    RecFlagSolver(P,V, {op(ProjArray(P,n))}, n-1, n, d);
    if nops(FLAGSOLS)>0 then print("Solution Found.  See FLAGSOLS");
        else print("No Solution Found.");
    fi;
    FLAGSOLS;
end: # FlagSolver

RecFlagSolver := proc(P,V,B, i, n, d)
description `recursively compute solutions to the board/vector equations, A=Pi, B=P_i+1`;
local NV,x,y,bb,Vsol,SolVectors,sols, sol,pt,A,R,Q,restorevariables, dfreedom:
global voo,FLAGSOLS,FIXEDVECTORS,FLAGVARIABLES:
    vprint(ELlevel,i);
    restorevariables:=NULL:
    A:={op(ProjArray(P,i))};
    B;
    Q:= A intersect B;
    R:= sort([op(A minus B)], lexorderlists);
    bb:=map((y)->V[y], Q);
    FIXEDVECTORS:=Basis([op(bb)]):
   #1.verify these vectors have rank <=i or solve
    if nops(FIXEDVECTORS) > i
    then
        vprint(fixed);
        sols:=SolveMatrixRankEqns(i,convert(FIXEDVECTORS,Matrix),{FLAGVARIABLES}):
        vprint(nops(sols)):
        for sol in sols do
            vprint(solution, sol):
            NV:=map(Rescale, SubsTable(sol, V)):
            vprint(rescaled);
            FIXEDVECTORS:=Basis(map((y)->NV[y], Q)):
            if nops(FIXEDVECTORS) <= i
            then RecFlagSolver(P,NV,B,i,n,d):
            fi:
        od:
    else
  #2.setup new vectors
        if nops(FIXEDVECTORS)=i then dfreedom:=0:
        else dfreedom:=i - nops(FIXEDVECTORS): fi:
        for x in R do
            V[x]:= VectorSet(x, Basis(map((y)->V[y], shadow(x,B)))):
            ### if extra degrees of freedom then add to FIXEDVECTORS
            if dfreedom > 0
            and Myrank(convert([op(FIXEDVECTORS),V[x]], Matrix)) >nops(FIXEDVECTORS)
            then FIXEDVECTORS:=[op(FIXEDVECTORS),V[x]]:
                FLAGVARIABLES:=op(remove((y)->y=CC[x][1], [FLAGVARIABLES])):
                restorevariables:=restorevariables, CC[x][1]:
                dfreedom:=dfreedom - 1:
            fi:
        od:
  #3.Find all viable vector tables of this index and proceed to i-1 or stop if i=1
        vprint(solving, i):
        SolVectors:=SolveAllRankEqns(i,A,V):
        for Vsol in SolVectors do
            if (i =1) then FLAGSOLS:= [op(FLAGSOLS),Vsol]:
            else FLAGVARIABLES:=FLAGVARIABLES, restorevariables:
                vprint(map((y)->V[y], A));
                RecFlagSolver(P,Vsol,A, i-1, n,d):
            fi:
        od:
    fi:

end: # RecFlagSolver

SolutionFlagBasis := proc (arr)
description ` For each solution vector table in FLAGSOLS, choose a basis from the vectors in the vector table.`;
local bases,dots,a,k,n,d,i,sol;
    bases:=NULL;
    n:=DotRank(arr,1);
    d:=nops(arr[1]);
    for i from 1 to n do
        a:=AntichainArray(arr, i);
        k:=min(op(map((x)->x[3], a)));
        a:=select((x)->x[3]=k, a);
        a:=sort(a, lexorderlists);
        dots[i]:=a[-1];
        od;
    dots:=map((x)->x[1..3],convert(dots,list));
    for sol in FLAGSOLS do
        bases:=bases,map((x)->Rescale(sol[x]),dots);
    od;
    [bases];
end: # SolutionFlagBasis


VectorSet := proc(x, abasis)
description `define vectors[pt] in the span of coving points`;
local v,w,i;
global FLAGVARIABLES;
    v:=abasis[1];
    for i from 2 to nops(abasis) do
        w:=abasis[i];
        v:=Add(v, w, 1, CC[x][i-1]);
        FLAGVARIABLES:=FLAGVARIABLES,CC[x][i-1];
    od;
    v;
end: # VectorSet



SubsTable := proc(eqs,A)    ## this seems to make a copy of the table A.
description `A=table, eqs= equations to sub into the elements of the table`;
        map((x)->simplify(subs(eqs, x)), A);
end: # SubsTable

SolveAllRankEqns := proc (i,bb,vectors)
description `Verify every principle subarray of bb has permrank equal to spanning-dim`;
local n,d,currentvectors,helper,j;
global TEMP;
    n:=Dimension(vectors[bb[1]]);
    d:=nops(bb[1]);
    currentvectors:=[vectors];
    #### run through all x in bb checking the points in shad(x)
    ### correspoind with vectors a vector set with same rank.
    helper := proc(tail)
    local shad;
        if ((d-1)*n < add(j, j=tail) and nops(currentvectors)>0 ) then
            shad:=shadow(tail,bb);
            ##insure all fixed 1-dim spaces along with v[tail] have at most rank i.
            if  (nops(shad) = 1) and member(tail,bb)
            then
                currentvectors:=FindSolutionVectorTables(eval(currentvectors),FIXEDVECTORS,tail, 1, i);
            fi;

            ##insure all shadows have given rank
            if  (nops(shad) > 1)
            then currentvectors:=FindSolutionVectorTables(currentvectors,
                                                     shad, tail, permrank(shad), permrank(shad)); fi;

        fi;
    end: # helper
    AllCoordinateTester(d,[seq(1,j=1..d)],[seq(n,j=1..d)],[],helper);
    #for cv in currentvectors do print("level",i); printvectors(bb, cv, n);  od;
    return(UniqueElems(currentvectors,EqualVectorTables));
end: # SolveRankEqns

UniqueElems := proc(alist,EqualityTest)
if nops(alist)>1 then
    if MyMember(alist[1],alist[2..-1],EqualityTest)
    then UniqueElems(alist[2..-1],EqualityTest)
    else [alist[1],op(UniqueElems(alist[2..-1],EqualityTest))]
    fi;
else alist; fi;
end: # UniqueElems

MyMember := proc(a,B,atest)
description `test if a is a member of a list B`;
local result,b;
    result:=false;
    for b in B do
        if apply(atest,a,b) then result:=true; break; fi;
    od;
    result;
end: # MyMember


EqualVectorTables := proc(VT1,VT2)   ###assumes equal set of indices
local result,a;
result:=true;
for a in indices(VT1) do
    if not(Equal(VT1[op(a)],VT2[op(a)])) then result:=false; return(false); fi;
od;
result;
end: # EqualVectorTables

FindSolutionVectorTables := proc (Vtables,B,tail,minrank,maxrank)
description `check given rank conditions on all vector tables in Vtables, solve if necessary`;
local M,eqns, sol,sols,V,newVtables,NV,rkM,n,y;
global FLAGEQNS,FLAGVARIABLES,MFOO;
    newVtables:=NULL;
    #for each currently valid set of vectors, verify the rank of the matrix of vectors corresponding to the elements in B is in the given interval
    if nops(B)=0 then     newVtables:=op(Vtables);
    else for V in Vtables do
            if type(B[1],Vector) then
                M:=convert([op(B),V[tail]], Matrix);
            else
                M:=convert(map((y)->V[y], [op(B)]), Matrix);            #matrix of corresponding vectors
            fi;
             if  ColumnDimension(M) > 5 then MFOO:=eval(M); fi;
             rkM:=Myrank(M);
#check minimum rank
             if (minrank > rkM) then print("Unexpected Dependencies Exist"); fi;
#check maximum rank
             if (rkM <= maxrank)  then newVtables:=newVtables,V;
             else
                 sols:=SolveMatrixRankEqns(maxrank,M,{FLAGVARIABLES});
                 for sol in sols do
                     NV:=map(Rescale, SubsTable(sol,V));
#check ranks again, if correct then add NV to newVtables
                     if type(B[1],Vector) then
                         M:=convert([op(B),NV[tail]], Matrix);
                     else
                         M:=convert(map((y)->NV[y], [op(B)]), Matrix);
                     fi;
                     rkM:=Myrank(M);
                     if (minrank > rkM) then print("Unexpected Dependencies Exist");
                     elif (rkM <= maxrank)  then newVtables:=newVtables,NV; fi;
                 od;
             fi;
         od;
    fi;
# return all vector table solutions whose rank is correct at tail.
    [newVtables];
end: # FindSolutionVectorTables

SolveMatrixRankEqns := proc(i,M,vars)
description `Solve equations which force the rank of M <= i`;
local eqns,sols;
global FLAGEQNS;
    eqns:={op(Minors(M,i+1))};
    sols:= map(RemoveTrivialEqns, [solve(eqns,vars)]);
    if nops(sols)>0 then FLAGEQNS:=sols,FLAGEQNS; fi;
    #sols:=map(allvalues, sols);
    #if nops(sols)>1 then print("multiple solutions:", sols);  fi;
    sols;
end: # SolveMatrixRankEqns

SelectTable := proc(sel,A)
local i,B,a;
i:=select(sel, [indices(A)]);
for a in i do B[a]:=eval(A[a]); od;
op(eval(B));
end: # SelectTable



###testing
foo:=[op(TransverseFlags(4,2)), BinomialFlag(4)];
doo:=[[4,4,2,2],[2,4,4,2],[4,4,1,3],[4,3,2,3],[3,4,2,3],[3,3,4,3],[4,2,3,4],[3,3,3,4],[2,4,3,4],[3,2,4,4],[2,3,4,4],[1,4,4,4], [4,1,4,1]];
#FlagVectorTable(foo);
#FlagSolver(doo,FVTABLE);
#printvectors(ProjArray(doo,1),FLAGSOLS[1],4);

davis:=[[2, 4, 4, 2], [4,1,4,2],[4, 2, 3, 3], [3, 4, 3, 3], [3, 2, 4, 3],
        [4, 3, 2, 4], [3, 3, 3, 4], [2, 4, 3, 4], [3, 4, 2, 4],
        [4,4,1,4], [2, 3, 4, 4], [1, 4, 4, 4], [4, 4, 2, 1]];


##########   Example:  uncomment to run.


#doo:=permarray([[1, 3, 2, 4], [2, 3, 4, 1], [3, 1, 2, 4]]);
#printdotarray(doo, 4,4);
#GenericFlag(4);
#foo:=[op(TransverseFlags(4,2)), GenericFlag(4)];
##foo:=[op(TransverseFlags(4,2)), RandomFlag(4)];  #uncomment for another choice of flag
#FlagSolver(doo, foo);
#SolutionFlagBasis(doo);
#printvectors(ProjArray(doo,1),FLAGSOLS[1],4);
#printvectors(ProjArray(doo,2),FLAGSOLS[1],4);
#printvectors(ProjArray(doo,3),FLAGSOLS[1],4);
#printvectors(ProjArray(doo,4),FLAGSOLS[1],4);