#####################################################################################
#####################################################################################
##                                  intHaar                                        ##
##            A Maple package for integrals over compact lie groups                ##
##                   Written by Alejandro Ginory and Jongwon Kim                   ##
##                               May 29th, 2016                                    ##
#####################################################################################
#####################################################################################

withSF(); with(ListTools); with(combinat); with(group); with(Statistics); with(GraphTheory); jacks(2);

print(`intHaar`):
print(`Following functions are supported for the integrals`):
print(`IntU(polynomial function in g[i,j] and conjugate(g[i,j]))`):
print(`IntO(polynomial function in g[i,j])`):
print(`IntSp(polynomial function in g[i,j], n :: dimension of Sp)`):

#####################################################################################
##              Helper Functions with self explanatory names                       ##
#####################################################################################

nonConjugate := proc (x) not has(x, conjugate) end proc; 

isPerm := proc (list) type(list, permlist) end proc; 

PermToPart := proc (P, n) local pt, df; 
   if evalb(nops(P) <> 0) then pt := [seq(nops(P[i]), i = 1 .. nops(P))] 
    else pt := [`$`(1, n)] end if; 
    df := n-add(pt); 
    sort([op(pt), seq(1, i = 1 .. df)], `>`); 
    end proc; 

PartToPowerSum := proc (part::list) local partl, var; 
    partl := [seq(cat(p, i), i in part)]; 
    var := mul(i, i in partl) 
    end proc; 

CoefZ := proc (P, e, vlist::list) local g, res; 
    if e = 1 then res := p else g := freeze(e); res := coeff(subs(e = g, P), g) end if;     
    eval(res, [seq(v = 0, v in vlist)]) 
    end proc; 

zl := proc (lambda) 
    mul(i^numboccur(lambda, i)*factorial(numboccur(lambda, i)), i = 1 .. max(lambda)); 
    end proc; 

fLambda := proc (lambda, n) Chi(Reverse(lambda), [`$`(1, n)]) end proc; 

CosetType := proc (sigma, n) local G, perm, edges; 
    G := Graph({seq({2*i, 2*i-1}, i = 1 .. n)}); 
    perm := convert(sigma, permlist, 2*n); 
    edges := {seq({perm[2*i], perm[2*i-1]}, i = 1 .. n)}; 
    G := ConnectedComponents(AddEdge(G, edges)); 
    [(1/2)*seq(nops(G[i]), i = 1 .. nops(G))] 
    end proc; 

mInterpreter := proc (m) convert(Flatten(m), disjcyc) end proc; 


#####################################################################################
##                                Unitary Base                                     ##
#####################################################################################

Schur1 := proc (lambda) 
    mul(mul(d+j-i, j = 1 .. lambda[i]), i = 1 .. nops(lambda)) 
    end proc; 

WeingartenU := proc (sigma, n) 
    simplify(add(Chi(Q, [`$`(1, n)])*Chi(Q, Reverse(sigma))/Schur1(Reverse(Q)), Q in partition(n))/factorial(n));
    end proc; 

GetPermU := proc (I1, I2) local check, path, T, n, P; 
    check := evalb(sort(I1) = sort(I2)); 
    if check then 
    path := [seq([SearchAll(I1[i], I2)], i = 1 .. nops(I1))]; 
    T := cartprod(path); 
    n := 0; 
    while not T[finished] do n := n+1; P[n] := `~`[op](T[nextvalue]()) end do; 
    P := map(convert, select(isPerm, convert(P, list)), disjcyc); 
    else 0 end if; 
    end proc; 

IUL := proc (I1, J1, I2, J2) local m, pI, pJ, int, i, j, P; 
    m := nops(I1); 
    pI := GetPermU(I1, I2); 
    pJ := GetPermU(J1, J2); 
    int := 0; 
    if type(pI, numeric) or type(pJ, numeric) then 0; 
    else 
    for i to nops(pI) do 
      for j to nops(pJ) do 
        int := int+WeingartenU(PermToPart(mulperms(pJ[j], invperm(pI[i])), m), m);
      end do; 
    end do;
    end if; 
    simplify(int); 
    end proc; 

#####################################################################################
##                              Orthogonal Base                                    ##
#####################################################################################

Zonal1 := proc (lambda) mul(mul(d+2*j-i-1, j = 1 .. lambda[i]), i = 1 .. nops(lambda)) 
    end proc; 

CoefNormalizerO := proc (rho, n) zl(2*rho)/(2^n*factorial(n)) end proc; 

OLS := proc (lambda, sigma) local n, j, P, vb, cf; 
    n := add(lambda); 
    j := J[seq(i, i in lambda)]; 
    P := subs(a = 2, top(j)); 
    vb := PartToPowerSum(sigma); 
    cf := CoefNormalizerO(sigma, n)*CoefZ(P, vb, [seq(cat(p, i), i = 1 .. n)]); 
    end proc; 

WeingartenO := proc (p1, p2, n) 2^n*factorial(n)*add(fLambda(2*lambda, 2*n)*OLS(lambda, CosetType(mulperms(p1, p2), n))/Zonal1(lambda), lambda in map(Reverse, partition(n)))/factorial(2*n); 
    end proc; 

GetPermO := proc (lst) local de, df, dp, T, n, P; 
    de := MakeUnique(lst); 
    df := [seq([SearchAll(i, lst)], i in de)]; 
    if convert(map(type, map(nops, df), even), `and`) then 
      dp := map(setpartition, df, 2); 
      T := combinat[cartprod](dp); 
      n := 0; 
      while not T[finished] do n := n+1; P[n] := `~`[op](T[nextvalue]()) end do; 
      P := map(`@`(mInterpreter, sort), convert(P, list)); else 0 
    end if; 
    end proc; 

IOL := proc (Il, Jl) local m, pI, pJ, int, i, j; 
    m := nops(Il); 
    pI := GetPermO(Il); 
    pJ := GetPermO(Jl); 
    int := 0; 
    if type(pI, numeric) or type(pJ, numeric) then 0 
    else 
      for i to nops(pI) do 
        for j to nops(pJ) do 
          int := int+WeingartenO(pJ[j], invperm(pI[i]), (1/2)*m);
        end do;
      end do;
    end if;
    simplify(int); 
    end proc;

#####################################################################################
##                              Symplectic Base                                    ##
#####################################################################################

Bijections := proc (A, B) 
    map[3](zip, `[]`, A, combinat:-permute(B)) 
    end proc; 

GetPermSp := proc (lst, m) local k, de, dI, dJ, check, dp, T, n, P; 
    k := nops(lst); 
    de := remove(proc (x) options operator, arrow; m < x end proc, MakeUnique(lst)); 
    dI := [seq([SearchAll(i, lst)], i = 1 .. m)]; 
    dJ := [seq([SearchAll(i+m, lst)], i = 1 .. m)]; 
    check := convert([seq(evalb(nops(dI[i]) = nops(dJ[i])), i = 1 .. m)], `and`); 
    if check then 
      dp := [seq(Bijections(dI[i], dJ[i]), i = 1 .. nops(dI))]; 
      T := combinat[cartprod](dp); 
      n := 0; 
      while not T[finished] do n := n+1; P[n] := `~`[op](T[nextvalue]()) end do; 
      P := map(`@`(mInterpreter, sort), convert(P, list)); 
    else 0 end if; 
    end proc; 

WeingartenSp := proc (p1, p2, n) 
    simplify(subs(d = -2*d, WeingartenO(p1, p2, n)*(-1)^n*parity(mulperms(p1, p2)))); 
    end proc; 

WgSign := proc (p1, p2, n) 
    sign(mul([seq((p1[2*i]-p1[2*i-1])*(p2[2*i]-p2[2*i-1]), i = 1 .. n)])); 
    end proc; 

ISpL := proc (Il, Jl, n) local m, pI, pJ, int, i, j; 
    m := (1/2)*nops(Il); 
    pI := GetPermSp(Il, n); 
    pJ := GetPermSp(Jl, n); 
    int := 0; 
    if type(pI, numeric) or type(pJ, numeric) then 0 
    else 
    for i to nops(pI) do 
      for j to nops(pJ) do 
      int := int+WgSign(convert(pJ[j], permlist, 2*m), convert(pI[i], permlist, 2*m), m)*WeingartenSp(pJ[j], invperm(pI[i]), m); 
      end do; 
    end do; 
    end if; 
    subs(d = n, simplify(int)); 
    end proc; 

#####################################################################################
##                                  IntHaar                                        ##
#####################################################################################

GetCoeff := proc (exp) local parts; 
    parts := [op(exp)]; 
    if has(exp, g) then 
    if type(exp, `^`) then 1 
    else mul(select(not has, parts, g)); 
    end if; 
    else exp 
    end if;
    end proc; 

OneVar := proc (mn) local parts, IJ, i; 
    parts := [op(mn)]; 
    IJ := []; 
    if type(mn, 'indexed') then parts 
    else 
    for i to parts[2] do 
      IJ := [op(IJ), op(parts[1])]; 
    end do; 
    end if; 
    end proc; 

MonomialToList := proc (mn) 
    if type(mn, `*`) then Flatten(map(OneVar, [op(mn)])) else OneVar(mn) 
    end if; 
    end proc; 

IOM := proc (mn) local cf, vr, IJ, Il, Jl; 
    cf := GetCoeff(mn); 
    vr := mn/cf; 
    if type(vr, 'numeric') then cf*vr 
    else 
    IJ := MonomialToList(vr); 
    Il := [seq(IJ[2*i-1], i = 1 .. (1/2)*nops(IJ))]; 
    Jl := [seq(IJ[2*i], i = 1 .. (1/2)*nops(IJ))]; 
    cf*IOL(Il, Jl); 
    end if; 
    end proc; 

IntO := proc (pn) if type(pn, `+`) then add(map(IOM, [op(pn)])) else IOM(pn) end if;
    end proc; 

ISpM := proc (mn, n) local cf, vr, IJ, Il, Jl; 
    cf := GetCoeff(mn); 
    vr := mn/cf; 
    if type(vr, 'numeric') then cf*vr 
    else 
    IJ := MonomialToList(vr); 
    Il := [seq(IJ[2*i-1], i = 1 .. (1/2)*nops(IJ))]; 
    Jl := [seq(IJ[2*i], i = 1 .. (1/2)*nops(IJ))]; 
    cf*ISpL(Il, Jl, n) 
    end if; 
    end proc; 

IntSp := proc (pn, n) 
    if type(pn, `+`) then add(map(ISpM, [op(pn)], n)) else ISpM(pn, n) end if; 
    end proc; 

hasAbs := proc (x) has(x, abs) end proc; 

absSelector := proc (x) 
    if type(x, `^`) and has(x, abs) then x else select(hasAbs, x) end if; 
    end proc; 

breaker := proc (x) 
    if type(x, `^`) then breakerH(x) else mul(map(breakerH, [op(x)])) end if; 
    end proc; 

breakerH := proc (x) 
    if type(op(x)[2], even) then (op(op(x)[1])*conjugate(op(op(x)[1])))^((1/2)*op(x)[2]) else 0 
    end if; 
    end proc; 

IUM := proc (mn) local cf, vr, as, nc, cc, IJ, Il, Jl, KL, Kl, Ll; 
    cf := GetCoeff(mn); 
    vr := mn/cf; 
    as := absSelector(vr); 
    if type(as, numeric) then 0 else vr := vr*breaker(as)/as end if; 
    if type(vr, 'numeric') then cf*vr 
    else 
    nc := select(nonConjugate, vr); 
    cc := conjugate(vr/nc); 
    IJ := MonomialToList(nc); 
    Il := [seq(IJ[2*i-1], i = 1 .. (1/2)*nops(IJ))]; 
    Jl := [seq(IJ[2*i], i = 1 .. (1/2)*nops(IJ))]; 
    KL := MonomialToList(cc); 
    Kl := [seq(KL[2*i-1], i = 1 .. (1/2)*nops(KL))]; 
    Ll := [seq(KL[2*i], i = 1 .. (1/2)*nops(KL))]; 
    cf*IUL(Il, Jl, Kl, Ll) 
    end if;
    end proc; 

IntU := proc (pn, n) 
    if type(pn, `+`) then add(map(IUM, [op(pn)])) else IUM(pn) end if;
    end proc;