(**This file computes our bound on the variation of the hessian.**)

(*First we load in the files having the relevant derivatives*)
<<derivatives1
<<derivatives2

(*This coordinate list specifies which variables are involved with which index, as
   discussed in the paper.  The mth item on the list gives the variables that are
   involved with the index m.  *)

coords={{1},{2,3},{4,5},{6,7},{1,2,3},{2,3,4,5},{2,3,6,7},{1,4,5},{1,6,7},{4,5,6,7}};

(**This routine returns 1 if all the variables of L={i,j,k} are involved with m.*)

involve0[L_,m_]:=(
a1=Table[Position[coords[[m]],L[[i]]],{i,1,Length[L]}];
Product[Length[a1[[i]]],{i,1,Length[L]}])

(**This routine relables the variable indices.  The energy function involves the
      indices 1,...,7, whereas each of the component summands only involves the
      indices 1,2 or 1,2,3,4.   Let's just write out an example to see what is going on.
      We have 

      E(x1,x2,x3,x4,x5,x6,x7) = ....  + U10(x4,x5,x6,x7).

      The variable list (4,5,6)  corresponds to the first three variables of U10,
      so the function involve1[{4,5,6},10] would return {1,2,3}*)


involve1[L_,m_]:=(
a1=involve0[L,m];
len=Length[L];
zero=Table[0,{len}];
a2=If[a1==0,zero,coords[[m]]];
a3=If[a1==0,zero,Table[Position[a2,L[[i]]],{i,1,len}]];
a4=If[a1==0,zero,Table[a3[[i,1,1]],{i,1,len}]])


(*This computes the absolute value of the derivative at the TBP.  These
    are the numbers a(i,m) and a(i,j,m) and a(i,j,k,m) in the paper.
    Here we call the function aa*)

aa[L_,m_]:=(
a1=involve1[L,m];
If[a1[[1]]==0,0,Abs[ExactBound[m,a1]]])


(*This is the function b from the paper.  Basically, it is the function a, plus a
    small fudge factor that depends on the choice of m*)

bb[L_,m_]:=(
a1=aa[L,m];
vec=If[m<5,{1/1000,0,0},{1/200,1/50,1/10}];
a2=a1+vec[[Length[L]]])


(*This computes the function Phi hat*)

PhiHat[{i_,j_,k_},m_]:=(
v1=bb[{i},m];
v2=bb[{j},m];
v3=bb[{k},m];
v12=bb[{i,j},m];
v13=bb[{i,k},m];
v23=bb[{j,k},m];
v123=bb[{i,j,k},m];
v={v123,v12 v3,v13 v2,v23 v1,v1 v2 v3};
w={a,b,b,b,c};
v.w)

(*This computes the double sum, which in the paper is advertised as less than 315.*)

Psi[i_,j_]:=(
Print[i," ",j];
a1=Table[Table[PhiHat[{i,j,k},m],{k,1,7}],{m,1,10}];
a2=Flatten[a1];
a3=Sum[a2[[jj]],{jj,1,Length[a2]}])

(*Here we compute the L2 norm of the resulting matrix.
This is the main routine.  The output is the expression FINAL,
a polynomial in a,b,c.  This expression should match FINAL2,
which is the polynomial for "Upsilon^2" that is advertised
in our Variation Bound II in the paper.  The routine FinalPoly
is fairly slow to run, but it outputs steadily to show its progress.*)

FinalPoly:=(
Clear[a,b,c];
matx=Table[Table[Psi[i,j],{i,1,7}],{j,1,7}];
m2=Flatten[matx];
m3=Expand[m2.m2];

(*m3 is the main expression, the square of the L2 norm.  Here we get the
   coefficients and round up to the nearest integer.*)

Caa=Coefficient[m3,a,2];
Cab=Coefficient[Coefficient[m3,a],b];
Cbb=Coefficient[m3,b,2];
Cac=Coefficient[Coefficient[m3,a],c];
Cbc=Coefficient[Coefficient[m3,b],c];
Ccc=Coefficient[m3,c,2];

Caa=Floor[Caa]+1;
Cab=Floor[Cab]+1;
Cbb=Floor[Cbb]+1;
Cac=Floor[Cac]+1;
Cbc=Floor[Cbc]+1;
Ccc=Floor[Ccc]+1;

(*Here is the final expression*)
FINAL=Caa a a  + Cab a b + Cbb b b+ Cac a c + Cbc b c + Ccc c c
)

(*The final expression should match the following expression,
which is in the paper.*)

FINAL2=19336 a a + 19036 a b + 4922 b b + 1474 a c + 772 b c + 31 c c

(**Now we substitute in the maximum derivatives and evaluate*)

(**Our final bound*)
FINAL