read("biquandles-maple.txt");
read("yangbaxtercohomology.txt");

gauss2vblist:=proc(g)
#
# get virtual biquandle presentation from virtual Gauss code
#
out:=[];
n:=coldim(matrix([g]));
test:=0;
for i from 1 to n do
 test:=test+g[i];
od;
if test <> 0 then 
 out:=false;
else
 p:=1;
 nxt:=[];
 z:=0;
 for i from 1 to n-1 do
  if g[i+1]=0 then
   nxt:=[op(nxt),eval(p)];
  elif g[i]=0 then
   z:=z+1;
   p:=i+1-z;
  else
   nxt:=[op(nxt),i+1-z];
  fi;
 od;
 g2:=[];
 for i from 1 to n do
  if g[i] <> 0 then
   g2:=[op(g2),g[i]];
 fi;
 od;
 ov:=[];
 for i from 1 to coldim(matrix([g2])) do
  for j from 1 to coldim(matrix([g2])) do
   if g2[i]+g2[j]=0 then
    ov:=[op(ov),j];
   fi;
  od;
 od;
 for i from 1 to coldim(matrix([g2])) do
  if Im(g2[i])=-2 then
   out:=[op(out),[S[eval(x[i])],eval(x[eval(nxt[i])])]];
  elif Im(g2[i])=2 then
   out:=[op(out),[Sinv[eval(x[i])],eval(x[eval(nxt[i])])]];
  elif Im(g2[i])=1 then
   out:=[op(out),[M[3][eval(x[i]),eval(x[eval(ov[i])])],eval(x[eval(nxt[i])])]];
  elif Im(g2[i])=-1 then
   out:=[op(out),[M[1][eval(x[i]),eval(x[eval(ov[i])])],eval(x[eval(nxt[i])])]];
  elif Im(g2[i])=0 and Re(g2[i])>0 then
   out:=[op(out),[M[4][eval(x[i]),eval(x[eval(ov[i])])],eval(x[eval(nxt[i])])]];
  elif Im(g2[i])=0 and Re(g2[i])<0 then
   out:=[op(out),[M[2][eval(x[i]),eval(x[eval(ov[i])])],eval(x[eval(nxt[i])])]];
  fi;
 od;
fi;
eval(out);
end;

reducepreslist:=proc(A)
#
# reduce biquandle presentation list
#
L:=[];
for i from 1 to coldim(matrix([A])) do 
 L:=[op(L),eval(A[i])];
od;
c:=true;
while c do
 c:=false;
 for i from 1 to coldim(matrix([L])) do
  w:=eval(L[1]);
  L:=subsop(1=NULL, L);
  if subs(eval(w[2])=z,eval(w[1]))=eval(w[1]) then 
   c:=true;
   L:=subs(eval(w[2])=eval(w[1]),L);
  else L:=[op(L),eval(w)];
  fi;
 od;
od;
eval(L);
end;

vbiqtest:=proc(A,S)
#
# Test whether (A,S) is a virtual biquandle
#
if biqtest(A) and permtest(S) and bhomtest(A,A,S) then
 eval(true);
else
 eval(false);
fi;
end;

vblist:=proc(A)
#
# lists all virtual biquandle structures on A
#
n:=coldim(A[1]);
L:=bautlist(A);
P:=permute(n);
out:=[L[1]];
for x in L do
 con:=false;
 for y in out do
  for p in P do
   tmp:=true;
   for i from 1 to n do
    if x[p[i]]<>p[y[i]] then
     tmp:=false;
    fi;
   od;
   if tmp then 
    con:=true;
   fi;
  od;
 od;
 if not con then
  out:=[op(out),eval(x)];
 fi;
od;
eval(out);
end;

listmaps := proc(a,b)
#
# list all a-component vectors with entries in 1..b
#
w:=initperm(b);
v:=[];
for i from 1 to a do
 v:=[op(v),w];
od;
ret:=[];
T:=cartprod(v);
while not T[finished] do
 t:=T[nextvalue]();
 ret:=[op(ret),t];
od;
eval(ret);
end;


vbhomlist:= proc(g,A,s)
#
# find virtual biquandle homomorphisms
#
sinv:=[];
n:=coldim(A[1]);
n2:=coldim(matrix([g]));
for i from 1 to n do
 for j from 1 to n do
  if s[j]= i then
   sinv:=[op(sinv),j];
  fi;
 od;
od;
L3:=gauss2vblist(g);
L2:=reducepreslist(L3);
L:=eval(L2);
for i from 1 to coldim(matrix([L])) do
 L:=subs(eval(L[i][2])=f[i],L);
od;
out:=[];
F:=listmaps(coldim(matrix([L])),coldim(A[1]));
for i from 1 to coldim(A[1])^coldim(matrix([L])) do
 w:=eval(F[i]);
 c:=true;
 for j from 1 to coldim(matrix([L])) do
  X:=L[j][1];
  Y:=L[j][2];
  if eval(subs(f=w,M=A,S=s,Sinv=sinv,eval(X))) <>
     eval(subs(f=w,M=A,S=s,Sinv=sinv,eval(Y))) then
   c:=false;
  fi;
 od; 
 if c then 
  out:=[op(out),eval(w)];
 fi;
od;
p:=1;
nxt:=[];
z:=0;
for i from 1 to n2-1 do
 if g[i+1]=0 then
  nxt:=[op(nxt),eval(p)];
 elif g[i]=0 then
  z:=z+1;
  p:=i+1-z;
 else
  nxt:=[op(nxt),i+1-z];
 fi;
od;
g2:=[];
for i from 1 to coldim(matrix([g])) do
 if g[i] <> 0 then
  g2:=[op(g2),g[i]];
 fi;
od;
ov:=[];
for i from 1 to coldim(matrix([g2])) do
 for j from 1 to coldim(matrix([g2])) do
  if g2[i]+g2[j]=0 then
   ov:=[op(ov),j];
  fi;
 od;
od;
out2:=[];
for w in out do
 w2:=[];
 for i from 1 to coldim(matrix([g2])) do
  w2:=[op(w2),x[i]];
 od;
 for i from 1 to coldim(matrix([L2])) do
  w2:=eval(subs(L2[i][2]=w[i],w2));
 od;
 c:=true;
 while c do
  c:=false;
  for i from 1 to coldim(matrix([w2])) do
   if Im(g2[i])=-2 and type(w2[i],integer) then
    w2[eval(nxt[i])]:=s[eval(w2[i])]; 
   elif Im(g2[i])=2 and type(w2[i],integer) then
     w2[eval(nxt[i])]:=sinv[eval(w2[i])];
   elif type(w2[i],integer) and type(w2[ov[i]],integer) then
    if Im(g2[i])=1 then
     w2[eval(nxt[i])]:=A[3][eval(w2[i]),eval(w2[ov[i]])];
    elif Im(g2[i])=-1 then
     w2[eval(nxt[i])]:=A[1][eval(w2[i]),eval(w2[ov[i]])];
    elif Im(g2[i])=0 and Re(g2[i])>0 then
     w2[eval(nxt[i])]:=A[4][eval(w2[i]),eval(w2[ov[i]])];
    elif Im(g2[i])=0 and Re(g2[i])<0 then
     w2[eval(nxt[i])]:=A[2][eval(w2[i]),eval(w2[ov[i]])];
    fi;
   else
    c:=true;
   fi;
  od;
 od;
 w3:=[];
 for i from 1 to coldim(matrix([g2])) do
  if g[i]=0 then
   w3:=[op(w3),0,w2[i]];
  else
   w3:=[op(w3),w2[i]];
  fi;
 od;
 w3:=[op(w3),0]; 
 out2:=[op(out2),w3];
od;
eval(out2);
end;

strongcompatcheck:=proc(A,S,phi)
#
# check for strong compatibility given compatibility
#
n:=coldim(A[1]);
c:=true;
for i from 1 to n do
 for j from 1 to n do
  if phi[n*(i-1)+j]<>phi[n*(S[i]-1)+S[j]] then
   c:=false;
  fi;
 od;
od;
eval(c);
end;

vbinv:=proc(g,A,S,phi,v)
#
# compute virtual 2-cocycle invariant
#
H:=vbhomlist(g,A,S);
Sinv:=[];
n:=coldim(A[1]);
n2:=coldim(matrix([g]));
for i from 1 to n do
 for j from 1 to n do
  if S[j]= i then
   sinv:=[op(Sinv),j];
  fi;
 od;
od;
p:=1;
nxt:=[];
z:=0;
for i from 1 to n2-1 do
 if g[i+1]=0 then
  nxt:=[op(nxt),eval(p)];
 elif g[i]=0 then
  z:=z+1;
  p:=i+1-z;
 else
  nxt:=[op(nxt),i+1-z];
 fi;
od;
g2:=[];
for i from 1 to coldim(matrix([g])) do
 if g[i] <> 0 then
  g2:=[op(g2),g[i]];
 fi;
od;
ov:=[];
for i from 1 to coldim(matrix([g2])) do
 for j from 1 to coldim(matrix([g2])) do
  if g2[i]+g2[j]=0 then
   ov:=[op(ov),j];
  fi;
 od;
od;
str:=strongcompatcheck(A,S,phi);
f:=0;
for w in H do
 vc:=0;
 pc:=0;
 w2:=[];
 for i from 1 to coldim(matrix([w])) do
  if w[i]<>0 then
   w2:=[op(w2),w[i]];
  fi;
 od;
 for i from 1 to coldim(matrix([g2])) do
  if Im(g2[i])=-2 then
   vc:=vc+v[n*(w2[i]-1)+w2[ov[i]]]-v[n*(w2[nxt[ov[i]]]-1)+w2[nxt[i]]];
  elif Im(g2[i])=-1 then
   pc:=pc-phi[n*(w2[nxt[i]]-1)+w2[nxt[ov[i]]]];
  elif Im(g2[i])=0 and Re(g2[i])<0 then
   pc:=pc+phi[n*(w2[i]-1)+w2[ov[i]]];
  fi;
 od;
 if str then
  f:=f+t^pc*s^vc;
 else 
  f:=f+T^(pc+vc);
 fi;
od;
eval(f);
end;


getkernel:=proc(N)
#
# read off a basis for the kernel of N
#
M:=gaussjord(N);
out:=[];
lp:=[];
lead:={}; # set of column numbers of leading ones
for i from 1 to rowdim(matrix(M)) do
 j:=1;
 while M[i,j]=0 and j < coldim(M) do
  j:=j+1;
 od;
 if M[i,j]=1 then
  lead:=lead union {eval(j)};
  lp:=[op(lp),eval(j)];
 fi;
od;
for i from 1 to coldim(M) do
 if not eval(i) in lead then 
  v:=[];
  for l from 1 to coldim(M) do
   v:=[op(v),0];
  od;
 for z from 1 to coldim(matrix([lp])) do
  v[lp[z]]:=eval(M[z,i]);
  od;
  v[i]:=-1;
  out:=[op(out),eval(v)];
 fi;
od;
eval(out);
end;


rmzeroes:=proc(M);
#
# remove zero rows
#
out:=[];
for i from 1 to rowdim(M) do
 c:=false;
 for j from 1 to coldim(M) do
  if M[i,j]<>0 then
   c:=true;
  fi;
 od;
 if c then
  m:=[];
  for j from 1 to coldim(M) do
   m:=[op(m),eval(M[i,j])];
  od;
  out:=[op(out),eval(m)];
 fi;
od;
eval(matrix(out));
end;


bautlist:=proc(A)
#
# list automorphisms of A
#
out:=[];
L:=bhomlist(A,A);
for x in L do
 if permtest(x) then
  out:=[op(out),x];
 fi;
od;
eval(out);
end;


ybcompattest:=proc(A,S,phi,v)
#
# test whether a virtual cocycle v is compatible with phi
#
n:=coldim(A[1]);
Sinv:=[];
for i from 1 to n do
 for j from 1 to n do
  if S[j]= i then
   Sinv:=[op(Sinv),j];
  fi;
 od;
od;
c:=true;
for i from 1 to n do
 for j from 1 to n do
  for k from 1 to n do
   if phi[n*(i-1)+j] 
     +v[n*(A[2][i,j]-1)+k]
     -v[n*(Sinv[k]-1)+S[A[2][i,j]]]
     +v[n*(A[4][j,i]-1)+Sinv[k]]
     -v[n*(Sinv[Sinv[k]]-1)+S[A[4][j,i]]] 
     <> phi[n*(S[i]-1)+S[j]] 
     +v[n*(j-1)+k]
     -v[n*(Sinv[k]-1)+S[j]]
     +v[n*(i-1)+Sinv[k]]
     -v[n*(Sinv[Sinv[k]]-1)+S[i]] then
    c:=false;
   fi;
  od;
 od;
od;
eval(c);
end;


sinv:=proc(S)
#
# Inverse of S
#
n:=coldim(matrix([S]));
Sinv:=[];
for i from 1 to n do
 for j from 1 to n do
  if S[j]= i then
   Sinv:=[op(Sinv),j];
  fi;
 od;
od;
eval(Sinv);
end;


vcocycles:=proc(S)
#
# find a basis for the space of rational S 2-cocycles
#
n:=coldim(matrix([S]));
Sinv:=[];
for i from 1 to n do
 for j from 1 to n do
  if S[j]= i then
   Sinv:=[op(Sinv),j];
  fi;
 od;
od;
r:=[];
for i from 1 to n^2 do
 r:=[op(r),0];
od;
M:=[];
for i from 1 to n do
 for j from 1 to n do
  for k from 1 to n do
   rij:=eval(r);
   rij[n*(i-1)+j]:=rij[n*(i-1)+j]-1;
   rij[n*(Sinv[j]-1)+S[i]]:=rij[n*(Sinv[j]-1)+S[i]]+1;
   rij[n*(S[i]-1)+S[j]]:=rij[n*(S[i]-1)+S[j]]+1;
   rij[n*(j-1)+S[S[i]]]:=rij[n*(j-1)+S[S[i]]]-1;
   rij[n*(S[i]-1)+k]:=rij[n*(S[i]-1)+k]-1;
   rij[n*(Sinv[k]-1)+S[S[i]]]:=rij[n*(Sinv[k]-1)+S[S[i]]]+1;
   rij[n*(i-1)+Sinv[k]]:=rij[n*(i-1)+Sinv[k]]+1;
   rij[n*(Sinv[Sinv[k]]-1)+S[i]]:=rij[n*(Sinv[Sinv[k]]-1)+S[i]]-1;
   rij[n*(Sinv[j]-1)+Sinv[k]]:=rij[n*(Sinv[j]-1)+Sinv[k]]-1;
   rij[n*(Sinv[Sinv[k]]-1)+j]:=rij[n*(Sinv[Sinv[k]]-1)+j]+1;
   rij[n*(j-1)+k]:=rij[n*(j-1)+k]+1;
   rij[n*(Sinv[k]-1)+S[j]]:=rij[n*(Sinv[k]-1)+S[j]]-1;
   M:=[op(M),eval(rij)];
  od;
 od;
od;
eval(getkernel(matrix(M)));
end;

vcob:=proc(S)
#
# Find S-coboundaries
#
out:=[];
n:=coldim(matrix([S]));
Sinv:=[];
for i from 1 to n do
 for j from 1 to n do
  if S[j]= i then
   Sinv:=[op(Sinv),j];
  fi;
 od;
od;
for i from 1 to n do
 cob:=[];
 for j from 1 to n do
  for k from 1 to n do
   temp:=0;
   if k=i then 
    temp:=temp+1;
   fi;
   if j=i then
    temp:=temp+1;
   fi;
   if S[j]=i then
    temp:=temp-1;
   fi;
   if Sinv[k]=i then
    temp:=temp-1;
   fi;
   cob:=[op(cob),eval(temp)];
  od;
 od;
 out:=[op(out),eval(cob)];
od;
M:=gaussjord(out);
out2:=[];
for i from 1 to rowdim(M) do
 j:=getrow(i,M);
 if not zerow(j) then
  out2:=[op(out2),eval(j)];
 fi;
od;
eval(out2);
end;

vcohom:=proc(S)
#
# Computes S-cohomology with Q coefficients
#
Z:=vcocycles(S);
Bn:=vcob(S);
out:=[];
while coldim(matrix([Z]))<>0 do
 x:=eval(Z[1]);
 Z:=subsop(1=NULL,Z);
 if not inspan(eval(x),Bn) then
  out:=[op(out),eval(x)];
  Bn:=[op(Bn),eval(x)]
 fi;
od;
eval(out);
end;

ybcompatfind:=proc(A,S)
#
# find linear combos of compatible S-cocycles and Yang-Baxter cocycles
#
n:=coldim(A[1]);
L:=vcocycles(S);
P1:=ybcocy(A);
P:=[];
for x in P1 do
 red:=true;
 for i from 1 to n do
  for j from 1 to n do
   if j=A[1][i,j] and i=A[3][j,i] then 
    if x[n*(j-1)+i]<>0 then
     red:=false;
    fi;
   fi;
   if i=A[2][j,i] and j=A[4][i,j] then 
    if x[n*(i-1)+j]<>0 then
     red:=false;
    fi;
   fi;
  od;
 od;
 if red then 
  P:=[op(P),x];
 fi;
od;
Sinv:=[];
for i from 1 to n do
 for j from 1 to n do
  if S[j]= i then
   Sinv:=[op(Sinv),j];
  fi;
 od;
od;
r:=[];
for i from 1 to coldim(matrix([L]))+coldim(matrix([P])) do
 r:=[op(r),0];
od;
vz:=[];
for i from 1 to coldim(matrix([L[1]])) do
 vz:=[op(vz),0];
od;
m:=[];
for i from 1 to n do
 for j from 1 to n do
  for k from 1 to n do
   rij1:=eval(r);
   for x from 1 to coldim(matrix([L])) do
    rij1[x]:= rij1[x]+L[x][n*(A[2][i,j]-1)+k]
                     -L[x][n*(Sinv[k]-1)+S[A[2][i,j]]]
                     +L[x][n*(A[4][j,i]-1)+Sinv[k]]
                     -L[x][n*(Sinv[Sinv[k]]-1)+S[A[4][j,i]]]
                     -L[x][n*(j-1)+k]
                     +L[x][n*(Sinv[k]-1)+S[j]]
                     -L[x][n*(i-1)+Sinv[k]]
                     +L[x][n*(Sinv[Sinv[k]]-1)+S[i]];
   od;
   for x from 1 to coldim(matrix([P])) do
    rij1[coldim(matrix([L]))+x]:=rij1[coldim(matrix([L]))+x]
                                -P[x][n*(i-1)+j]+P[x][n*(S[i]-1)+S[j]];
   od; 
   m:=[op(m),rij1];
  od;
 od;
od;
M:=rmzeroes(gaussjord(matrix(m)));
if coldim(M)=0 then 
 out:=false;
else
 out:=[];
  SolSet:=getkernel(M);
  for x in SolSet do
   v:=eval(vz);
   for j from 1 to coldim(matrix([L])) do
    v:=v+x[j]*L[j]; 
   od;
   phi:=eval(vz);
   for j from 1 to coldim(matrix([P])) do
    phi:=phi-x[j+coldim(matrix([L]))]*P[j];
   od; 
   dc:=true;
   for j from 1 to coldim(matrix([phi])) do
    if phi[j]<>0 then
     dc:=false;
    fi;
   od;
  if dc then
   if not degcheck(S,v) then
    out:=[op(out),[phi,v]];
   fi;
  else
   out:=[op(out),[phi,v]];
  fi;
 od;
fi;
eval(out);
end;


degcheck:=proc(S,v)
#
# Check for degeneracy
#
n:=coldim(matrix([S]));
Sinv:=[];
for i from 1 to n do
 for j from 1 to n do
  if S[j]= i then
   Sinv:=[op(Sinv),j];
  fi;
 od;
od;
c:=true;
for i from 1 to n do
 for j from 1 to n do
  if v[n*(i-1)+j]<>v[n*(Sinv[j]-1)+S[i]] then
   c:=false;
  fi;
 od;
od;
eval(c);
end;