(* %%%%%%%%%%%%%%%%%%  *)
(* Graphs on 4 vertices *)
(* %%%%%%%%%%%%%%%%%% *)

(* List of all graphs *)
allEdges[n_]:=Flatten[Table[UndirectedEdge[i,j],{i,1,n},{j,i+1,n}]];
allConnected[n_]:=Select[Map[Graph[Range[n],#,VertexLabels->None]&,Subsets[allEdges[n]]],ConnectedGraphQ];
allConnectedUpToIso[n_]:=DeleteDuplicates[allConnected[n],IsomorphicGraphQ];
allConnectedUpToIso[4]

(* Path P_4 *)
(* Constants computed using Maximize *)
Clear[f,g,i,j,k,x1,x2,x3,x4];
av[1]=(({{1, 0, 0, 0}, {1/2, 1/2, 0, 0}, {1/3, 1/3, 1/3, 0}, {1/4, 1/4, 1/4, 1/4}}));
av[2]=(({{0, 1, 0, 0}, {1/3, 1/3, 1/3, 0}, {1/4, 1/4, 1/4, 1/4}}));
av[3]=(({{0, 0, 1, 0}, {0, 1/3, 1/3, 1/3}, {1/4, 1/4, 1/4, 1/4}}));
av[4]=(({{0, 0, 0, 1}, {0, 0, 1/2, 1/2}, {0, 1/3, 1/3, 1/3}, {1/4, 1/4, 1/4, 1/4}}));
var[{x1_,x2_,x3_,x4_}]=Abs[x1-x2]+Abs[x2-x3]+Abs[x3-x4];
g[k_]:=Max[av[k].{x1,x2,x3,x4}];
Maximize[{var[{g[1],g[2],g[3],g[4]}],var[{x1,x2,x3,x4}]==1},{x1,x2,x3,x4}]

(* Rigorous verification for P_4 *)
values={};
con={x1>=x2>=x3>=x4};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x2>=x4>=x3,x4<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x2>=x4>=x3,x4>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x1+x2-2 x3>=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x1+x2-2 x3<=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x4>=x2,(1/3) (x1+x2+x3)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x4>=x2,(1/3) (x1+x2+x3)>=(1/4) (x1+x2+x3+x4),((x3+x4)/2)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x4>=x2,(1/3) (x1+x2+x3)>=(1/4) (x1+x2+x3+x4),((x3+x4)/2)<=(1/4) (x1+x2+x3+x4),x1+x2-2 x3>=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x4>=x2,(1/3) (x1+x2+x3)>=(1/4) (x1+x2+x3+x4),((x3+x4)/2)<=(1/4) (x1+x2+x3+x4),x1+x2-2 x3<=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x4>=x2>=x3,x4<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x4>=x2>=x3,x4>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x4>=x3>=x2,x1+x2+x3>=3 x4};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x4>=x3>=x2,x1+x2+x3<=3 x4};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x3>=x4,(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x3>=x4,(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x4>=x3,(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x4>=x3,(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4),x4<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x4>=x3,(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4),x4>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x1>=x4,((x1+x2)/2)<=(1/3) (x1+x2+x3),(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x1>=x4,((x1+x2)/2)<=(1/3) (x1+x2+x3),(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x1>=x4,((x1+x2)/2)>=(1/3) (x1+x2+x3),(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x1>=x4,((x1+x2)/2)>=(1/3) (x1+x2+x3),(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x4>=x1,((x1+x2)/2)>=(1/3) (x1+x2+x3)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x4>=x1,((x1+x2)/2)<=(1/3) (x1+x2+x3),(1/3) (x1+x2+x3)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x4>=x1,((x1+x2)/2)<=(1/3) (x1+x2+x3),(1/3) (x1+x2+x3)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x1>=x3,(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x1>=x3,(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4),x4<=(1/3) (x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x1>=x3,(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4),x4>=(1/3) (x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x3>=x1,x4<=(1/3) (x2+x3+x4),((x1+x2)/2)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x3>=x1,x4<=(1/3) (x2+x3+x4),((x1+x2)/2)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x3>=x1,x4>=(1/3) (x2+x3+x4),((x1+x2)/2)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x3>=x1,x4>=(1/3) (x2+x3+x4),((x1+x2)/2)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
Print[Length[values]," linear optimization problems"];
Print["All local maxima: ",Union[values]];
Print["Global maximum = ",Max[values]];

(* Cycle C_4 *)
(* Constant computed using Maximize *)
Clear[f,g,i,j,k,x1,x2,x3,x4];
av[1]=(({{1, 0, 0, 0}, {1/3, 1/3, 0, 1/3}, {1/4, 1/4, 1/4, 1/4}}));
av[2]=(({{0, 1, 0, 0}, {1/3, 1/3, 1/3, 0}, {1/4, 1/4, 1/4, 1/4}}));
av[3]=(({{0, 0, 1, 0}, {0, 1/3, 1/3, 1/3}, {1/4, 1/4, 1/4, 1/4}}));
av[4]=(({{0, 0, 0, 1}, {1/3, 0, 1/3, 1/3}, {1/4, 1/4, 1/4, 1/4}}));
var[{x1_,x2_,x3_,x4_}]=Abs[x1-x2]+Abs[x2-x3]+Abs[x3-x4]+Abs[x4-x1];
g[k_]:=Max[av[k].{x1,x2,x3,x4}];
Maximize[{var[{g[1],g[2],g[3],g[4]}],var[{x1,x2,x3,x4}]==1},{x1,x2,x3,x4}]

(* Rigorous verification for C_4 *)
values={};
con={x1>=x2>=x3>=x4,(1/3) (x1+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x2>=x3>=x4,(1/3) (x1+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x2>=x4>=x3};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x3<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x3>=(1/4) (x1+x2+x3+x4),x1-2 x3+x4>=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x3>=(1/4) (x1+x2+x3+x4),x1-2 x3+x4<=0,x1+x2-2 x3<=0,(1/3) (x1+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x3>=(1/4) (x1+x2+x3+x4),x1-2 x3+x4<=0,x1+x2-2 x3<=0,(1/3) (x1+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x3>=(1/4) (x1+x2+x3+x4),x1-2 x3+x4<=0,x1+x2-2 x3>=0,(1/3) (x1+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x3>=(1/4) (x1+x2+x3+x4),x1-2 x3+x4<=0,x1+x2-2 x3>=0,(1/3) (x1+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
Print[Length[values]," linear optimization problems"];
Print["All local maxima: ",Union[values]];
Print["Global maximum = ",Max[values]];

(* Paw Y *)
(* Constant computed using Maximize *)
Clear[f,g,i,j,k,x1,x2,x3,x4];
av[1]=(({{1, 0, 0, 0}, {1/4, 1/4, 1/4, 1/4}}));
av[2]=(({{0, 1, 0, 0}, {1/2, 1/2, 0, 0}, {1/4, 1/4, 1/4, 1/4}}));
av[3]=(({{0, 0, 1, 0}, {1/3, 0, 1/3, 1/3}, {1/4, 1/4, 1/4, 1/4}}));
av[4]=(({{0, 0, 0, 1}, {1/3, 0, 1/3, 1/3}, {1/4, 1/4, 1/4, 1/4}}));
var[{x1_,x2_,x3_,x4_}]=Abs[x1-x2]+Abs[x1-x3]+Abs[x1-x4]+Abs[x3-x4];
g[k_]:=Max[av[k].{x1,x2,x3,x4}];
Maximize[{var[{g[1],g[2],g[3],g[4]}],var[{x1,x2,x3,x4}]==1},{x1,x2,x3,x4}]

(* Rigorous verification for Y *)
values={};
con={x1>=x2>=x3>=x4,(1/3) (x1+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x2>=x3>=x4,(1/3) (x1+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,(1/3) (x1+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,(1/3) (x1+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x4>=x2,((x1+x2)/2)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x4>=x2,((x1+x2)/2)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
con={x2>=x1>=x3>=x4,x1<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x3>=x4,x1>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x1>=x4,x1>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x1>=x4,x1<=(1/4) (x1+x2+x3+x4),x3>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x1>=x4,x1<=(1/4) (x1+x2+x3+x4),x3<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x4>=x1,x3<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x3>=x4>=x1,x3>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x1>=x2>=x4,(1/3) (x1+x3+x4)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x1>=x2>=x4,(1/3) (x1+x3+x4)>=(1/4) (x1+x2+x3+x4),x1<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x1>=x2>=x4,(1/3) (x1+x3+x4)>=(1/4) (x1+x2+x3+x4),x1>=(1/4) (x1+x2+x3+x4),2 x1-x3-x4<=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x1>=x2>=x4,(1/3) (x1+x3+x4)>=(1/4) (x1+x2+x3+x4),x1>=(1/4) (x1+x2+x3+x4),2 x1-x3-x4>=0,((x1+x2)/2)<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x1>=x2>=x4,(1/3) (x1+x3+x4)>=(1/4) (x1+x2+x3+x4),x1>=(1/4) (x1+x2+x3+x4),2 x1-x3-x4>=0,((x1+x2)/2)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x1>=x4>=x2,x1<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x1>=x4>=x2,x1>=(1/4) (x1+x2+x3+x4),2 x1-x3-x4<=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x1>=x4>=x2,x1>=(1/4) (x1+x2+x3+x4),2 x1-x3-x4>=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x2>=x1>=x4,x1>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x2>=x1>=x4,x1<=(1/4) (x1+x2+x3+x4),x2>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x2>=x1>=x4,x1<=(1/4) (x1+x2+x3+x4),x2<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x2>=x4>=x1,x2<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x2>=x4>=x1,x2>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x4>=x1>=x2,x1<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x4>=x1>=x2,x1>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x4>=x2>=x1,x2<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x3>=x4>=x2>=x1,x2>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
Print[Length[values]," linear optimization problems"];
Print["All local maxima: ",Union[values]];
Print["Global maximum = ",Max[values]];

(* Diamond D *)
(* Constant computed using Maximize *)
Clear[f,g,i,j,k,x1,x2,x3,x4];
av[1]=(({{1, 0, 0, 0}, {1/3, 1/3, 0, 1/3}, {1/4, 1/4, 1/4, 1/4}}));
av[2]=(({{0, 1, 0, 0}, {1/4, 1/4, 1/4, 1/4}}));
av[3]=(({{0, 0, 1, 0}, {0, 1/3, 1/3, 1/3}, {1/4, 1/4, 1/4, 1/4}}));
av[4]=(({{0, 0, 0, 1}, {1/4, 1/4, 1/4, 1/4}}));
var[{x1_,x2_,x3_,x4_}]=Abs[x1-x2]+Abs[x1-x4]+Abs[x2-x3]+Abs[x2-x4]+Abs[x3-x4];
g[k_]:=Max[av[k].{x1,x2,x3,x4}];
Maximize[{var[{g[1],g[2],g[3],g[4]}],var[{x1,x2,x3,x4}]==1},{x1,x2,x3,x4}]

(* Rigorous verification for D *)
values={};
con={x1>=x2>=x3>=x4,x2<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x2>=x3>=x4,x2>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x2>=x4>=x3,x2<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x2>=x4>=x3,x2>=(1/4) (x1+x2+x3+x4),x4<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x2>=x4>=x3,x2>=(1/4) (x1+x2+x3+x4),x4>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x2>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x2<=(1/4) (x1+x2+x3+x4),x3<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x1>=x3>=x2>=x4,x2<=(1/4) (x1+x2+x3+x4),x3>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x3>=x4};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x4>=x3,(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x4>=x3,(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4),x4<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x1>=x4>=x3,(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4),x4>=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x1>=x3,x4<=(1/4) (x1+x2+x3+x4)};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x1>=x3,x4>=(1/4) (x1+x2+x3+x4),(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4),x1+x2-2 x4<=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x1>=x3,x4>=(1/4) (x1+x2+x3+x4),(1/3) (x2+x3+x4)<=(1/4) (x1+x2+x3+x4),x1+x2-2 x4>=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x1>=x3,x4>=(1/4) (x1+x2+x3+x4),(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4),x1+x2-2 x4<=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x1>=x3,x4>=(1/4) (x1+x2+x3+x4),(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4),x1+x2-2 x4>=0,x2+x3-2 x4<=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
con={x2>=x4>=x1>=x3,x4>=(1/4) (x1+x2+x3+x4),(1/3) (x2+x3+x4)>=(1/4) (x1+x2+x3+x4),x1+x2-2 x4>=0,x2+x3-2 x4>=0};
con2=Simplify[var[{x1,x2,x3,x4}]==1,Assumptions->con];
obj=Simplify[var[{g[1],g[2],g[3],g[4]}],Assumptions->con];
v=obj/.LinearOptimization[-obj,Join [con,{con2}],{x1,x2,x3,x4}];
values=Append[values,v];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
Print["Maximize ",obj," subject to ",Simplify[Join [con,{con2}]]," returns ",v,"."];
Print[Length[values]," linear optimization problems"];
Print["All local maxima: ",Union[values]];
Print["Global maximum = ",Max[values]];

(* %%%%%%%%%%%%%%%%%%  *)
(* Graphs on 5 vertices *)
(* %%%%%%%%%%%%%%%%%% *)

(* List of all graphs *)
allConnectedUpToIso[5]

(* Constants computed using Maximize *)
allgraphs=allConnectedUpToIso[5];
For[it=1,it<=Length[allgraphs],it++,
adj=AdjacencyMatrix[allgraphs[[it]]];
dis=GraphDistanceMatrix[allgraphs[[it]]];
vf=Sum[If[adj[[i,j]]==1,Abs[x[i]-x[j]],0],{i,1,5},{j,i+1,5}];
gv={};
For[i=1,i<=5,i++,
tmv={};
For[r=0,r<=4,r++,
bl=0;tms=0;
For[j=1,j<=5,j++,
If[dis[[i,j]]<=r,
bl=bl+1;tms=tms+x[j]
]];
tmv=Append[tmv,tms/bl]
];
gv=Append[gv,tmv]
];
g=Simplify[Max/@gv];
vg=Sum[If[adj[[i,j]]==1,Abs[g[[i]]-g[[j]]],0],{i,1,5},{j,i+1,5}];
Print[Maximize[{vg,x[1]>=0,x[2]>=0,x[3]>=0,x[4]>=0,x[5]>=0,vf==1},{x[1],x[2],x[3],x[4],x[5]}]]
];