#
# Procedures for testing the validity of the database of matrix models:
#
# checkreps   - verify all models in the database for W(R)
# checkone    - confirm that a list of matrices is a matrix model for W(R)
# checktraces - confirm that a rep'n of W(R) has the correct traces
# classtest   - certify a set of "discriminating" conjugacy classes

read `coxeter2.4v.txt`; withcoxeter():

read sparseops;
read unpack;

if `+`(0)=0 then # we are using Maple V.4 or later
  myspace:=proc() kernelopts(bytesalloc) end
else # we don't have to deal with the latest "improvements"
  myspace:=proc() 4*status[2] end
fi;

# checkreps - verify all models in the database for W(R)
#
# After the user loads one of the files of models via a read command,
#  checkreps(R) will unpack each of the orthogonal and rational seminormal
#  models for W(R) in that file, and use 'checkone' and 'checktraces' to
#  verify that they obey the Coxeter relations and have the right character.
# Copious information about the time cost of each test is lprinted.
# After the i-th model has been tested, a line of the form "i, N, t, sp"
#  is printed, where N is the dimension, t is the amount of time used to
#  test the i-th model, and sp is the amount of RAM allocated (in bytes).
# To save time, only a subset of the character values are checked--by
#  default, only those corresponding to para-Coxeter classes. To specify
#  an alternate set of classes, supply a sublist of [1,2,...,m] (where m is
#  the number of conjugacy classes of W(R)) as a second argument. 
# Use coxeter[class_rep](R) to see the listing of conjugacy classes.
# The first step in this program is to run 'classtest' to see if the chosen
#  list of conjugacy classes is sufficiently large to uniquely identify the
#  irreducible characters of W(R). If this test fails, no further checking
#  is done and an ERROR is generated.
# If any of the models fails one of the integrity checks, an ERROR is
#  generated and the program halts immediately.
# WARNING: testing the orthogonal models for E8 takes a long time.
#
# Uncomment the following line to omit testing of the orthogonal models:

#seminormal_only:=true;

checkreps:=proc(R) local n,st,st0,st00,a,cc,ct,i;
  interface(quiet=true);
  st00:=time(); n:=coxeter['rank'](R);
  if not assigned(Rep[n]) then
    ERROR(`load models data first`)
  elif nargs>1 then cc:=args[2] else
    cc:=coxeter['class_rep'](R);
    cc:=[seq([i,cc[i]],i=1..nops(cc))];
    cc:=map(proc(x) if nops(x[2])=nops({op(x[2])}) then x[1] fi end,cc);
  fi;
  ct:=coxeter['irr_chars'](R);
  st:=time(); if classtest(R,cc) then
    lprint(`class list is sufficient`,time()-st); printf(`\n`) else
    ERROR(`insufficent classes to detect irreps; try a different list`)
  fi;
  for i to nops(ct) do
    st0:=time();
    if not seminormal_only=true then
      st:=time(); a:=unpack(i,n);
      lprint(`orthogonal unpacking`,time()-st);
      st:=time(); checkone(a,R);
      lprint(`orthogonal checking`,time()-st);
    fi;
    st:=time(); a:=unpack(i,n,'semi');
    lprint(`seminormal unpacking`,time()-st);
    st:=time(); checkone(a,R);
    lprint(`seminormal checking`,time()-st);
    st:=time(); checktraces(a,ct[i],R,cc);
    lprint(`character checking`,time()-st);
    print(i,nops(a[1]),time()-st0,myspace());
  od;
  interface(quiet=false);
  R,time()-st00,myspace();
end;

# checkone - confirm that a list of matrices is a matrix model for W(R)
#
# If  rep  is an alleged matrix model for W(R) (orthogonal or seminormal),
#  in the format produced by the 'unpack' function, checkone(rep,R) will
#  verify that the generators satisfy the Coxeter relations of W(R).
# In the orthogonal case, also check that the matrices are orthogonal.
# In the seminormal case, also check that the matrices preserve the correct
#  diagonal form (the n+1-th matrix in the given list of matrices).

checkone:=proc(rep,R) local n,N,m,j,k,v,fl,f;
  n:=coxeter['rank'](R);
  m:=coxeter['cox_matrix'](R);
  if nops(rep)=n then fl:=NULL else fl:='semi' fi;
  N:=nops(rep[1]); v:=[seq(cat('e',j),j=1..N)];
  for j to n do
    f:=braid_eqns(1,fast_prod(rep[j],rep[j]),v,fl);
    if simplify(f,'sqrt')<>{0} then ERROR(cat(`not involution `,j)) fi;
    if fl<>NULL then f:=fast_prod(rep[j],rep[n+1]) else f:=rep[j] fi;
    if not is_sym(f) then ERROR(cat(`does not preserve form `,j)) fi;
    for k to j-1 do
      f:=braid_eqns(m[j,k],rep[j],rep[k],fl);
      if simplify(f,'sqrt')<>{0} then ERROR(cat(`bad braid `,j,` `,k)) fi;
    od
  od; true
end;

# checktraces - confirm that a rep'n of W(R) has the correct traces
#
# If  rep  is a rational matrix model for W(R), given as a list of matrices
#  representing the generators (in the format produced by 'unpack'), and f
#  is a character of W(R) (a list of traces), checktraces(rep,f,R) will
#  verify that the character of 'rep' is f.
# If a fourth argument L is supplied, where L is a sublist of [1,2,...,m]
#  (m = number of conjugacy classes of W(R)), then only the traces on
#  conjugacy classes indexed by L are tested.
# Use coxeter[class_rep](R) to see the listing of conjugacy classes.
# We do not assume that the model is hereditary, nor that the ordering of
#  simple reflections is graceful: traces are computed naively.

checktraces:=proc(rep,f,R) local cr,cc,j,w,k,tr;
  cr:=coxeter['class_rep'](R);
  if nargs>3 then cc:=args[4] else cc:=[$1..nops(f)] fi;
  for j in cc do w:=cr[j];
    if nops(w)=0 then tr:=nops(rep[1]) else tr:=rep[w[1]];
      for k from 2 to nops(w) do tr:=fast_prod(tr,rep[w[k]]) od;
      tr:=convert([seq(coeff(tr[k],cat('e',k)),k=1..nops(tr))],`+`);
    fi;
    if f[j]<>tr then ERROR(cat(`wrong trace at class `,j)) fi
  od; true
end;

# classtest - certify a set of "discriminating" conjugacy classes
#
# classtest(R,cc) uses linear programming methods to test whether a
#  subset of the conjugacy classes of W(R) are sufficient to distinguish
#  every irreducible character of W(R) from among all characters of W(R).
# In other words, suppose that f and g are known to be characters of
#  representations of W(R), and that f is irreducible. If the values of
#  of f and g agree on a particular subset of conjugacy classes, must f=g?
# The set of conjugacy classes is specified as a sublist cc of [1,...,m],
#  where m is the number of conjugacy classes of W(R).
# Use coxeter[class_rep](R) to see the listing of conjugacy classes.
# If cc is omitted, the default is to use the list of para-Coxeter classes.
# Add the optional argument 'verbose' to see the list of irreps of W(R)
#  that fail to be uniquely determined by the chosen conjugacy classes.
#
# Note: for A_n, B_n (n odd), D_n, E6, E7, E8, the para-Coxeter classes
#  are sufficient, but there can be much smaller discriminating sets.

classtest:=proc(R) local N,ct,cc,i,j,f,c,lp,dims,flag;
  flag:=-1; N:=nargs;
  if args[N]='verbose' then flag:=0; N:=N-1 fi;
  ct:=coxeter['irr_chars'](R);
  if N>1 then cc:=args[2] else
    cc:=coxeter['class_rep'](R);
    cc:=[seq([j,cc[j]],j=1..nops(cc))];
    cc:=map(proc(x) if nops(x[2])=nops({op(x[2])}) then x[1] fi end,cc);
  fi;
  dims:=[seq([j,ct[j][1]],j=1..nops(ct))];
  f:={seq(convert([seq(c[j]*ct[j][i],j=1..nops(ct))],`+`)=0,i=cc)};
  for i to nops(ct) do
    lp:=map(proc(x,y) if x[2]>y then x[1] fi end,dims,ct[i][1]);
    lp:=subs({seq(c[j]=0,j=lp),c[i]=-1},f);
    lp:=traperror(simplex['feasible'](lp,NONNEGATIVE)); # trap a bug in V.3
    if lp=lasterror and lp<>'infeasible' then ERROR(lasterror)
      elif lp<>lasterror and lp then
        if flag<0 then RETURN(false) else flag:=flag+1;
          lprint(cat(`irrep `,i,` fails the classtest`)) fi;
    fi;
  od;
  evalb(flag<1)
end;