#
# Unpack orthogonal and seminormal data for matrix reps of Weyl groups.
#
# Requires a file of representation data to be loaded first.
#
# With two arguments, unpack(i,n) returns a list [O[1],...,O[n]] of
# orthogonal matrices for the i-th irrep of the Weyl group of rank n in
# the (current) database. The j-th matrix represents the action of s[j].
#
# With three arguments, unpack(i,n,<anything>) returns a list
# [A[1],...,A[n],D] of (rational) seminormal matrices for the i-th irrep,
# along with the (positive, rational) diagonal matrix D that provides
# an invariant quadratic form (i.e., A[j] * D * A[j]^T = D).
# Note that the positive square root of D converts between the orthogonal
# and seminormal forms:  D^(1/2) * O[j] * D^(-1/2) = A[j].
#
# For more information, see the accompanying READ_ME,
# or visit <www.math.lsa.umich.edu/~jrs/data/models>
 
unpack:=proc(i,n) local ch,k,T,di;
  if not assigned(Rep[i,n]) or (nargs>2 and not (assigned(Semi[n])
    and member(i,Semi[n]))) then ERROR(`not available`) fi;
  ch:=listchains(i,0,n); di:=NULL;
  if nargs>2 then
    di:=[seq(cat('e',k),k=1..nops(ch))];
    for k from 2 to n do
      di:=zip((x,y)->x*y,di,[seq(Semi[T[k],T[k+1],k]^2,T=ch)])
    od
  fi;
  [seq(restore(k,n,ch,args[3..nargs]),k=1..n),di];
end;

# NOTE: if a new file of models is read during the same Maple session,
#  it will overwrite the previous data, but it will not unpack correctly
#  until a remember table is erased.
# Use this command (with zero arguments) to erase the remember table:

`unpack/erase`:=proc() global `restore/block`;
  `restore/block`:=subsop(4=NULL,op(`restore/block`)); NULL
end:

############################ Internal Stuff #################################
#
# List all chains in the Bratteli diagram that start at level k and
# end at vertex i of level n. These index a basis for a representation
# of the parabolic subgroup of W[n] generated by those simple reflections 
# that centralize W[k].

listchains:=proc(i,k,n) local ch,j,T,r;
  ch:=[[i]];
  for j from n by -1 to k+1 do
    ch:=[seq(seq([r,op(T)],r=Branch[T[1],j]),T=ch)]
  od; ch
end:

# Given a list ch of partial chains of fixed length ending at rank n that
# is stable under the action of s[k], restore the matrix of this action.
# Note: this requires that the chains have >= n+1-Rep[0,k] terms each.

# With a fourth argument, restore the seminormal version.

restore:=proc(k,n,ch) local a,N,r,i,k0,k1,T0,T1,v,vars;
  a:=table(); r:=nops(ch[1]); N:=nops(ch);
  k0:=Rep[0,k]-n+r; k1:=k-n+r;
  for i to N do
    T0:=op(1..k0-1,ch[i]); T1:=op(k1+1..r,ch[i]);
    vars:=Rep[ch[i][k1],k][ch[i][k0]];
    if nargs=3 then a[i]:=Rep[[op(k0..k1,ch[i])],k] else
      member([op(k0..k1,ch[i])],vars,'v');
      a[i]:=`restore/block`(ch[i][k0],ch[i][k1],k)[v];
    fi;
    a[i]:=subs({seq(e[op(v)]=e[T0,op(v),T1],v=vars)},a[i]);
  od;
  a:=[seq(a[i],i=1..N)];
  subs({seq(e[op(ch[i])]=cat('e',i),i=1..N)},a);
end:

# Restore the j-th seminormal block for the action of s[n] on
# the i-th irrep of W[n], and remember it.

`restore/block`:=proc(j,i,n) local ch,sc,rep,T,k;
  option remember;
  ch:=Rep[i,n][j]; rep:=[seq(Rep[T,n],T=ch)];
  sc:=[seq(convert([seq(Semi[T[k-1],T[k],n-nops(T)+k],
    k=2..nops(T))],`*`),T=ch)];
  rep:=subs(zip((x,y)->e[op(x)]=e[op(x)]/y,ch,sc),rep);
  rep:=zip((x,y)->simplify(x*y,'sqrt'),sc,rep);
  if type(rep,'list(polynom(rational))') then rep
    else ERROR(`restoration failed -- corrupt data`) fi;
end: