#
# Write a collection of reports about the models in the database:
#
#  canonical_names
#  clone_list
#  complexity
#  quality
#
# Simply read this file into a Maple session and prepare to wait.
# Running time on a 2.8GHz Pentium IV: about 8 minutes.
#
# WARNING: this will overwrite existing files in the 'stats' directory,
# and you will lose the scintillating comments that were in the originals.
#
# WARNING^2: If you are using a Maple version older than Maple 6, make sure
# you have first saved a copy of the E8_data file in the `.m` format that
# your Maple version can read:
#
#  read `../data/E8_data`; save Rep,Semi,Branch,`../data/E8_data.m`; quit;
#

read `coxeter2.4v.txt`; withcoxeter():
read orthogonal;
read clone_test;
read sparseops;
read unpack;

#####################
#
# canonical_names
#
# Print the "canonical" name of each irrep of W(R).
# Format: [N,t,e] or [N,t1,t2,e], where N=dimension,
#  t = trace of a reflection; t1,t2=traces of either conjugacy class
#  of reflections, if there is more than one,
#  and e = sign of the trace of w0.
# Note that this fails to distinguish two irreps of W(F4),
#  and fails sometimes in the classical cases (e.g., W(A13)).
# Also print the number of eigenvalues of s[1] that equal 1,
#  as a percentage. (A useful parameter for calculations involving the
#  spherical unitary dual.)

foo:=proc(a) if a<0 then `-` elif a=0 then `O` else `+` fi end;
score:=proc(f,i1) .01*round(50*(f[1]+f[i1])/f[1]) end;

repnames:=proc(R) local S,n,ct,cr,w0,i,j,i1,i2,formatstring;
  interface(quiet=true);
  S:=coxeter['base'](R); n:=nops(S);
  ct:=coxeter['irr_chars'](S);
  cr:=coxeter['class_rep'](S);
  w0:=coxeter['longest_elt'](S);
  member(coxeter['class_rep'](w0,S),cr,'j');
  member(coxeter['class_rep']([1],S),cr,'i1');
  member(coxeter['class_rep']([n],S),cr,'i2');
  if i1=i2 then
    formatstring:=`%3d  [%4d,%4d, %1s] %4.2f\n`;
    for i to nops(ct) do
      printf(formatstring,i,ct[i][1],ct[i][i1],
        foo(ct[i][j]),score(ct[i],i1))
    od;
  else
    formatstring:=`%3d  [%4d,%4d,%4d, %1s] %4.2f\n`;
    for i to nops(ct) do
      printf(formatstring,i,ct[i][1],ct[i][i1],ct[i][i2],
        foo(ct[i][j]),score(ct[i],i1))
    od;
  fi;
  interface(quiet=false); printf(`\n`);
end;

writeto(`../stats/canonical_names`);
for R in [G2,F4,E5,E6,E7,E8] do R; repnames(R) od;
writeto(terminal);

#####################
#
# clonelist
#
# Generate a list of irreps of W(R) that have nontrivial clones.
# In each case, return a list of clones; each clone is a list of
#  integers indicating a set of W[n]-irreps whose direct sum, when
#  restricted to W[n-1], is isomorphic to the W[n-1] restriction of
#  the given irrep.

listclones:=proc(R) local S,n,S0,i; global Rep,Branch;
  interface(quiet=true);
  S:=coxeter['base'](R); n:=nops(S);
  S0:=subsop(n=NULL,S);
  Rep[n]:=coxeter['irr_chars'](S);
  Rep[n-1]:=coxeter['irr_chars'](S0):
  Branch:=table();
  for i to nops(Rep[n]) do clone_test(i,S0,S) od;
  interface(quiet=false); printf(`\n`);
end;

writeto(`../stats/clonelist`);
for R in [F4, E4, E5, E6, E7, E8] do R; listclones(R) od;
writeto(terminal);

#####################
#
# complexity
#
# Report on the number of equations and variables involved in the system
#  that determines each hereditary model of W(R), assuming that s[n]
#  centralizes W[n-2], s[n-1] centralizes W[n-3], and that s[n]s[n-1] and
#  s[n-1]s[n-2] both have order 3.
# Also, mark with *'s the cases that are not totally free.

complexity:=proc(R) local ct,S,n,S1,S2,S3,f1,f2,f3,i,a,b,q;
  interface(quiet=true);
  S:=coxeter['base'](R); n:=nops(S);
  S1:=[op(1..n-1,S)];
  S2:=[op(1..n-2,S)];
  S3:=[op(1..n-3,S)];
  ct:=coxeter['irr_chars'](S);
  f1:=[seq(cat('e',i),i=1..nops(ct))];
  f3:=map(coxeter['cprod'],coxeter['irr_chars'](S3),
    coxeter['restrict'](f1,S3,S),S3);
  f2:=map(coxeter['cprod'],coxeter['irr_chars'](S2),
    coxeter['restrict'](f1,S2,S),S2);
  f1:=map(coxeter['cprod'],coxeter['irr_chars'](S1),
    coxeter['restrict'](f1,S1,S),S1);
  for i to nops(ct) do
    q:=max(op(fast_prod(f1,ct[i])));
    if q>1 then q:=`*` else q:=` ` fi;
    a:=convert(map(x->x*(x+1)/2,fast_prod(f3,ct[i])),`+`);
    b:=convert(map(x->x*(x+1)/2,fast_prod(f2,ct[i])),`+`);
    printf(`%3d %1s%4d  %5d  %3d\n`,i,q,ct[i][1],a+b,b);
  od;
  interface(quiet=false);
end;

writeto(`../stats/complexity`);
for R in [E6, E7, E8] do R; complexity(R) od;
writeto(terminal);

#####################
#
# quality
#
# Measure the sparsity of a matrix (average number of nonzero's per row)

sparsity:=proc(a) local sp;
  sp:=convert(map(x->nops(indets(x)),a),`+`);
  .01*round(100*sp/nops(a));
end;

# Report MN,MD,LCD for the entries of a matrix
#   MN=maximum (absolute) numerator
#   MD=maximum denominator
#  LCD=least common denominator

quality:=proc(a) local lcd,tm;
  tm:=map(x->x^2,map(coeffs,{op(a)}));
  lcd:=sqrt(ilcm(op(map(denom,tm))));
  sqrt(max(op(map(numer,tm)))),sqrt(max(op(map(denom,tm)))),lcd;
end;

# Report the dimension, sparsity, and quality (both orthogonal and
# seminormal forms) of s[n] in every irrep of W[n].

report:=proc(R) local n,i,sp,ch,q,a,j,b;
  n:=coxeter['rank'](R);
  if not assigned(Rep[n]) then ERROR(`no model data loaded`) fi;  
  for i to nops(Rep[n]) do
    ch:=listchains(i,0,n); b:=restore(n,n,ch,0);
    sp:=sparsity(b); b:=[quality(b)];
    a:=[quality(restore(n,n,ch))];
    if min(op(Branch[i,n]))<0 then q:=`*` else q:=` ` fi;
    a:=subs({seq(j=cat('r',j^2),j=indets(a,'sqrt'))},a);
    a:=map(proc(x) if type(x,`*`) then cat(op(1,x),`*`,op(2,x))
      else x fi end,a);
    a:=cat(`[`,a[1],`, `,a[2],`, `,a[3],`]`);
    b:=cat(`[`,b[1],`, `,b[2],`, `,b[3],`]`);
    printf(`%3d %1s%4d %5.2f %25s %25s\n`,i,q,nops(ch),sp,
      convert(a,'string'),convert(b,'string'))
  od;
  printf(`\n`);
end;

writeto(`../stats/quality`);
interface(quiet=true);
read `../data/F4_data`;
F4; report(F4);

`unpack/flush`();    # clear out the remember table
read `../data/E8_data.m`;
for R in [E6,E7,E8] do R; report(R) od;
interface(quiet=false);

writeto(terminal);

quit;