#
# For an irrep of W[n] that is *not* multiplicity-free with respect to
#  the action of W[n-1] (but has multiplicities no worse than two),
#  find a change of basis that optimizes the quality of the diagonal
#  entries in the orthogonal model.
# Reference: Section 4B of "Explicit matrices for irreducible..."
#
# Set infolevel[orthopt]:=2 to see info about the optimization steps
# Requires sparseops, unpack, and Rep data to be loaded

orthopt:=proc(i,n) local ch,f,rep,N,a,b,j,eq,u,r; global Rep;
  if not assigned(Rep[i,n]) then ERROR(`build the rep first`) fi;
  ch:=listchains(i,Rep[0,n],n);
  rep:=map(proc(x,y) if x[y]<0 then -x[y] fi end,ch,n-Rep[0,n]);
  rep:=[op({op(rep)})]; b:=restore(n,n,ch); N:=nops(ch);
  for r in rep do
    a:=map(`orthopt/act`,[$1..N],ch,r,n-Rep[0,n],u);
    b:=fast_prod(a,fast_prod(b,a));
    a:=zip((x,y)->coeff(x,cat('e',y)),b,[$1..N]);
    f:=`orthopt/bruteforce`(a,u);
    b:=simplify(subs({u[1]=f[1],u[2]=f[2]},b),'sqrt');
  od;
  eq:={seq(cat('e',j)=e[op(ch[j])],j=1..N)};
  for j to N do Rep[ch[j],n]:=subs(eq,b[j]) od; `Done`
end;

# Do a brute force search over positive rationals that parameterize the
# space of renormalizations (we can skip negative rationals, since that
# amounts to interchanging equivalent pairs of W[n-1]-submodules).

`orthopt/bruteforce`:=proc(f,x) local a,b,j,k,rad,ht,h,p,z,z0;
  a:=map(proc(y) if not type(y,'rational') then y fi end,f);
  rad:=convert(indets(a,'sqrt'),`*`)^2;
  a:=subs({x[1]=1/sqrt(1+rad*z^2),x[2]=z*sqrt(rad)/sqrt(1+rad*z^2)},a);
  a:=map(factor,simplify(a,'sqrt'));
  z0:=0; ht:=ilcm(op(map(denom,subs(z=0,a))));
  userinfo(2,orthopt,cat(`old denom `,ht));
  for k from 2 to 100 do for j to k-1 do 
    h:=denom(subs(z=j/(k-j),a[1]));
    for p in [2,3,5,7] do while irem(h,p,'b')=0 do h:=b od od;
    if h=1 then
      b:=subs(z=j/(k-j),a); h:=ilcm(op(map(denom,b)));
      if h<ht then ht:=h; z0:=j/(k-j);
        userinfo(2,orthopt,cat(`new denom `,ht,` at `),z0) fi
    fi
  od od;
  [1/sqrt(1+rad*z0^2),z0*sqrt(rad)/sqrt(1+rad*z0^2)];
end:

`orthopt/act`:=proc(k,ch,r,p,x) local T,l;
  if abs(ch[k][p])<>r then cat('e',k)
    else T:=subsop(p=-ch[k][p],ch[k]);
    if not member(T,ch,'l') then ERROR(`this cannot happen`) fi;
    if ch[k][p]=r then x[1]*cat('e',k)+x[2]*cat('e',l)
      else x[2]*cat('e',l)-x[1]*cat('e',k) fi
  fi
end: