(* Sample code of Numerical Bootstrap in 2109.02701[hep-th] by Takeshi Morita *)
(* This code is for Mathematica ver.12 or later. *)

(*Define the Hamiltonian*)
(* Here, op[m,n]:=exp(imx)*p^n. Note that op[0,0]=1. *)

θ=0;
a=1/2;

H=op[0,2]/2+θ*op[0,1]/2/Pi+θ^2/8/Pi*op[0,0]-a*((op[1,0]+op[-1,0])/2-op[0,0])

(* 
Define the optimal function:
Write an equation that you want to minimize. 
For example, optimal=a*op[k,l]+b*op[n,m]+...
*)

optimal=op[1,0]

(*Set the parameters for the numerical bootstrap.*)

Kx=4
Kp=3

(* Set the range of the energy you want to test.
The energy is taken as the range from "eneini" to "enemax" divided by "Ne". 
 *)

eneini=0
enemax=2.5

Ne=100

(* Optional parameter *)
(* If you have any trouble, try setting K larger. *)

K=2*Max[{Kx,Kp}]+4;



(* -------------------MAIN CODE---------------------- *)
(* You do not need to change the following code. *)

(* Set some parameters. *)

enelist=Subdivide[eneini,enemax,Ne-1];
ansatz={op[0,0]==1,op[1,0]==xvev,op[0,1]==pvev};

(* DEFINE FUNCTIONS *)

(* Taking the operators op[m,n] from an equation. *)

constants={a,θ,ene,xvev,pvev};
vars[eq_]:=Complement[Variables[eq],constants]
(* Define a product of two operators: op[k,l]*op[m,n] *)

OPprod[op1_,op2_]:=Module[{k,l,m,n},
k=op1[[1]];
l=op1[[2]];
m=op2[[1]];
n=op2[[2]];
If[m==0,
op[k,l+n],
Sum[l!/r!/(l-r)!*m^r*op[k+m,l-r+n],{r,0,l}]
]]


(* 
Define the product of two equations. 
eq1, eq2 should be taken the form: 
eq1=a*op[k,l]+b*op[m,n]+....
*)

EQprod[eq1_,eq2_]:=Module[{ops1,ops2,coeff1,coeff2},
ops1=vars[eq1];
coeff1=Coefficient[eq1,ops1];
ops2=vars[eq2];
coeff2=Coefficient[eq2,ops2];
Expand[Flatten[Outer[Times,coeff1,coeff2]].Flatten[Outer[OPprod,ops1,ops2]]]
]

(*Define the commutator: [op[k,l],op[m,n]] *)

com[eq1_,eq2_]:=Expand[EQprod[eq1,eq2]-EQprod[eq2,eq1]]

(*CONSTRUCTING THE BOOTSTRAP MATRIX*)
(*Define the seed operators for the bootstrap matrix *)

ops=Flatten[Table[op[i,j],{i,0,Kx},{j,0,Kp}]];

(*Define the Hermite conjugate of "ops". *)

ops2=Flatten[Table[OPprod[op[0,j],op[-i,0]],{i,0,Kx},{j,0,Kp}]];

(*Define the bootstrap matrix*)
bootstrapmatrix=Table[EQprod[ops2[[i]],ops[[j]]],{i,Length[ops]},{j,Length[ops]}];

(*THE CONSTRAINT EQUATIONS*)

(* Deriving the constraint equation [H,op[m,n]]=0 and [H,op[-m,n]]=0. *)

Ex=K;
Ep=K;


Eq1=Table[{com[H,op[m,n]],com[H,op[-m,n]]},{m,0,Ex},{n,0,Ep}];

(* Deriving the constraint equation H*op[m,n]= ene*op[m,n] and H*op[-m,n]= ene*op[-m,n]. *)


Eq2=Table[{EQprod[H,op[m,n]]-ene*op[m,n],EQprod[H,op[-m,n]]-ene*op[-m,n]},{m,0,Ex},{n,0,Ep}];

(* Prepare the constraint equations *)

Eqs={};
For[n=1,n<K+1,n++,For[m=1,m<K+1-n,m++,AppendTo[Eqs,Eq1[[m,n]]]]]
For[n=1,n<K,n++,For[m=1,m<K-n,m++,AppendTo[Eqs,Eq2[[m,n]]]]]
Eqs=DeleteDuplicates[Flatten[Eqs]];
eqvars=vars[Eqs];
Eqs2=Thread[Eqs==0];
PrependTo[Eqs2,ansatz];


(* Solving the constraint equations. *)

Print["Solving the constraint equations: HO=EO and [H,O]=0."]
sol=Solve[Flatten[Eqs2],eqvars][[1]];
(* Testing the solutions of the constraint equations. *)

bootstrapmatrix2=Expand[bootstrapmatrix/.sol];
Print["Is the bootstrap matrix Hermitian?"]
If[HermitianMatrixQ[bootstrapmatrix2/.Thread[Variables[bootstrapmatrix2]->RandomReal[]]],Print["Yes."],
Print["No. Try setting K larger."]]

Print["Is the number of the constraint equations sufficient?"]
If[Length[Complement[Variables[bootstrapmatrix2],Variables[constants]]]==0
,Print["Yes."]
,Print["No. Try setting K larger."]]

(* Calcutations for the optimization. *)

Print["Computing the optimization."]
boot=Table[
Module[{rule,BootstrapMatrix,Bootstrapvars,constraint,result},
rule={ene->enelist[[i]]};

BootstrapMatrix=bootstrapmatrix2/.rule;
Bootstrapvars=Variables[BootstrapMatrix];
constraint= VectorGreaterEqual[{BootstrapMatrix,0}, {"SemidefiniteCone",Length[BootstrapMatrix[[1]]]}];
result=Quiet[SemidefiniteOptimization[optimal/.sol/.rule,{constraint},Bootstrapvars]];
{enelist[[i]],optimal/.sol/.rule/.result}],{i,1,Length[enelist]}];
boot2={};
Do[If[NumberQ[ib[[2]]],AppendTo[boot2,ib]],{ib,boot}]

(* Plot the result. *)

plboot=ListPlot[boot2,AxesLabel->{"E","<"<>ToString[optimal]<>">"},PlotRange->All]