[Based on SDPB.m by David Simmons-Duffin and a notebook by Kay Wiese.]

In the numerical conformal bootstrap we solve Polynomial Matrix Programs (PMP). This means we want to find a functional α such that α · V[x] is positive for all vectors V from some list and for all positive x. Such a PMP can be converted to a Semi-Definite-Problem (SDP) which in turn can be solved by SDPB.

Historically SDPB used to take .xml files and therefore as an intermediate step the PMP gets written as a .xml file. Nowadays SDPB also takes JSON so this step can be skipped. 

Either way some preprocessing of the files is required to convert the PMP to the SDP format that SDPB accepts as input.

Some terminology:
	Dual Feasible: a functional exists that is semi-positive definite on all polynomial matrix vectors.
	Primal Feasible: No functional will be found. A point cannot be excluded using the current set of assumptions and constraints.
So points that are allowed are Primal Feasible and points that are excluded are Dual Feasible.
The term feasible (without being preceded by Primal or Dual) is also used and means that a solution to crossing could exist.

Installation instructions

I guess this notebook is not going to work on Windows since you won’t have access to bash? There is probably some equivalent thing to make it work. If worst comes to worst you might have to run Ubuntu in a virtual environment. (I did this once and can give some pointers. If you really want to get working with numerics bash is a must.)

Installing Docker +SDPB +scalar_blocks

Installing Docker:
1.) Install Docker. (For example for Ubuntu 20.04: https://linuxhint.com/install_configure_docker_ubuntu/). If you are using Mac be careful to not use the latest version and instead use say 3.2.2. There is a bug in the newer versions.
2.) Set up permissions:
https://docs.docker.com/engine/install/linux-postinstall/#manage-docker-as-a-non-root-user 
(otherwise you will encounter the error “Got permission denied while trying...”).
3.) Set up SDPB (and scalar_blocks) with docker pull wlandry/sdpb:2.4.0
See also the instructions on https://github.com/davidsd/sdpb/blob/master/docs/Docker.md
The code below assumes that you are using the suggested directory /usr/local/share/sdpb/ as the mount point in Docker for the directory with your notebook and output folder. If you want to put it somewhere else just replace instances of the above with your directory.

4.) Some code below uses a non-standard shell that you might need to install, probably with : sudo apt-get install tcsh
If this is a problem for anyone I can fix a work around without it.

Initialization

User directories

Please set outputDataDirectory to a directory path where you want to store output (blocks and SDPB files). Also set dockerExecutable to something that runs docker. If docker is in your $PATH simply putting “docker” should work. If not put the explicit path to docker there.

(* user defined directories*)
outputDataDirectory=FileNameJoin@{NotebookDirectory[],"OutputDir"}; (*You can pick any directory for this. By default I put it in the directory of this notebook.*)
Quiet[CreateDirectory[outputDataDirectory],CreateDirectory::filex];

sdpb="sdpb:2.4.0";  (* the sdpb version we want to use, 2.4.0 is the most recent as of 4/2021 *)

(*dockerExecutable="/usr/local/bin/docker";*)
(* If you have your paths set up correctly docker should be callable simply by docker.*)
dockerExecutable="docker";

SetDirectory[NotebookDirectory[]];
(* Print more if verbatim is set to True*)
VPrint[X___]:=If[verbatim,Print[X]]

Parameters for the numerics

(* the precision for all our calculations in bits *)
precision=768;
(* the same in digits of 10, plus a little reserve, only used by Mathematica (does not majorly affect any timings) *)
digits=Round[N[precision Log[2]/Log[10]]+5]; 
(* physical variables *)
dim=3;(* or N[30/10,digits]*) (* the dimension *)
delta12=delta34=0; (* for 4 identical scalar fields *)

(* spin range*)
spinrangeMathematica=Range[0,20];
(* the same in scalar_blocks notation*)
spinrange=StringRiffle[ToString/@MinMax@Range[0,20],"-"];
nmax=6; (* include derivatives upto n+m<2nmax-1=Λ*)
(* variables for the generation of the conformal blocks *)
poles=18;
order=55; (* needs to be sufficiently high. Used for generation of the conformal blocks *)

(* internal variables for sdpb *)
(*termination conditions:*)
maxIterations=500; (*terminate if not converged after this many points*)
maxRuntime=864000; (* terminate after this many seconds.*)
findPrimalFeasible=1; (* terminate if primal solution found, i.e. the point is allowed*)
detectPrimalFeasibleJump=1;  (*terminate if the primal error sharply decreases in one step*)
findDualFeasible=1; (*terminate if a functional excluding the point is found*)
detectDualFeasibleJump=1;(*terminate if the dual error sharply decreases in one step*)
maxComplementarity=(10^100); (*terminate if the task gets too far from converging. This indicates the numerics are unstable.*)

(* a problem is considered dual if the dual error is smaller than this:*)
primalErrorThreshold=(10^(-70)); (* should be small, although ideally we always terminate on so called "jumps"*)
(* a problem is considered primal if the primal error is smaller than this:*)
dualErrorThreshold=(10^(-20));

dualityGapThreshold=(10^(-5)); (* terminate if both the primal and dual objective converge to the same number within this threshold*)

checkpointInterval=3600; (* frequency of writing progress to a check point*)
noFinalCheckpoint=1;
procsPerNode=6;      (* the number of kernels to use *)

(* some parameters to improve numeric stability:*)
initialMatrixScalePrimal=(10^30);
initialMatrixScaleDual=(10^30);
feasibleCenteringParameter=0.1;
infeasibleCenteringParameter=0.3;
stepLengthReduction=0.7;

(* routine to write parameter file param.sdpb used by sdpb *)
sdpbParamWrite[]:=Module[{aux},
aux={"maxIterations = "<>ToString[maxIterations],
"maxRuntime = "<>ToString[maxRuntime],"checkpointInterval = "<>ToString[checkpointInterval],
"noFinalCheckpoint = "<>ToString[noFinalCheckpoint],
"findPrimalFeasible = "<>ToString[findPrimalFeasible],
"detectPrimalFeasibleJump = "<>ToString[detectPrimalFeasibleJump],
"findDualFeasible = "<>ToString[findDualFeasible],
"detectDualFeasibleJump = "<>ToString[detectDualFeasibleJump],
"precision = "<>ToString[precision],
"procsPerNode = "<>ToString[procsPerNode],
"dualityGapThreshold = "<>ToString[dualityGapThreshold//N//CForm],
"primalErrorThreshold = "<>ToString[primalErrorThreshold//N//CForm],
"dualErrorThreshold = "<>ToString[dualErrorThreshold//N//CForm],
"initialMatrixScalePrimal = "<>ToString[initialMatrixScalePrimal//N//CForm],"initialMatrixScaleDual = "<>ToString[initialMatrixScaleDual//N//CForm],"feasibleCenteringParameter = "<>ToString[feasibleCenteringParameter],"infeasibleCenteringParameter = "<>ToString[infeasibleCenteringParameter],"stepLengthReduction = "<>ToString[stepLengthReduction],"maxComplementarity = "<>ToString[maxComplementarity//N//CForm]};
Export[FileNameJoin@{outputDataDirectory,"param.sdpb"},StringRiffle[aux,"\n"],"Text"];
]

(* For an optimization task with a meaningful objective we use different settings, most notably we want to continue until the primal and dual objectives converge *)
sdpbParamWriteOptimization[]:=Module[{aux},
Block[{findPrimalFeasible=0,detectPrimalFeasibleJump=0,findDualFeasible=0,detectDualFeasibleJump=0},
aux={"maxIterations = "<>ToString[maxIterations],
"maxRuntime = "<>ToString[maxRuntime],"checkpointInterval = "<>ToString[checkpointInterval],
"noFinalCheckpoint = "<>ToString[noFinalCheckpoint],
"findPrimalFeasible = "<>ToString[findPrimalFeasible],
"detectPrimalFeasibleJump = "<>ToString[detectPrimalFeasibleJump],
"findDualFeasible = "<>ToString[findDualFeasible],
"detectDualFeasibleJump = "<>ToString[detectDualFeasibleJump],
"precision = "<>ToString[precision],
"procsPerNode = "<>ToString[procsPerNode],
"dualityGapThreshold = "<>ToString[dualityGapThreshold//N//CForm],
"primalErrorThreshold = "<>ToString[(10^(-20))//N//CForm],
"dualErrorThreshold = "<>ToString[dualErrorThreshold//N//CForm],
"initialMatrixScalePrimal = "<>ToString[initialMatrixScalePrimal//N//CForm],"initialMatrixScaleDual = "<>ToString[initialMatrixScaleDual//N//CForm],"feasibleCenteringParameter = "<>ToString[feasibleCenteringParameter],"infeasibleCenteringParameter = "<>ToString[infeasibleCenteringParameter],"stepLengthReduction = "<>ToString[stepLengthReduction],"maxComplementarity = "<>ToString[maxComplementarity//N//CForm]};
Export[FileNameJoin@{outputDataDirectory,"param_optimization.sdpb"},StringRiffle[aux,"\n"],"Text"];
]
]

sdpbParamWriteOptimization[]
sdpbParamWrite[]

Loading derivatives

cf[x_]:=CForm[x];
error[msg_] := (Print[msg]; Exit[];);
floatToString[x_]:=If[Abs[x]<10^(-10),"0", ToString[CForm[SetAccuracy[x,10]]]];
(* Fail with an error message if a file doesn't exist *)
safeGet[file_] := If[FileExistsQ[file], Get[file], error[file <> " does not exist. Exiting."];];
safeImport[file_] := If[FileExistsQ[file], Import[file], error[file <> " does not exist. Exiting."];];  
derivPolTable[d_, delta12_, delta34_, L_, nmax_, keptPoleOrder_, order_] :=
    safeImport[
        StringJoin[
                   FileNameJoin@{outputDataDirectory,"output"},
                   "/zzbDerivTable-d", ToString[d],
                   "-delta12-",        floatToString[delta12],
                   "-delta34-",        floatToString[delta34],
                   "-L",               ToString[L],
                   "-nmax",            ToString[nmax],
                   "-keptPoleOrder",   ToString[keptPoleOrder],
                   "-order",           ToString[order],
                   ".m"]];

Procedure for creating vectors with derivatives of F

(* Computation Setup *)

seriesDataToRules[ruleName_,series_] := Module[
    {toRule},
    toRule[expr_,{m_,n_}] := ruleName[m-1,n-1]->(m-1)!(n-1)!expr;
    toRule[expr_,_]       := expr;
    Flatten[MapIndexed[toRule,Normal[series],2]]
];


withDeltaPhiToWithoutDeltaPhi[DeltaPhi_][mMax_,nMax_] := Module[
    {series, half},
    half = SetPrecision[1/2,precision];
    series = Series[
        ((1-z) (1-zb))^DeltaPhi f[z,zb]/.{z->half+dz,zb->half+dzb},
    {dz,0,mMax},
    {dzb,0,nMax}
    ] // Normal;

	seriesDataToRules[withDeltaPhiDeriv,
		CoefficientList[series,{dz,dzb}]/.{
			Derivative[j_,k_][f][_,_]:>zzbDeriv[Max[j,k],Min[j,k]],
			f[_,_]:>zzbDeriv[0,0]
		}
    ]
];

oddDerivs[nMax_]  := Flatten[Table[withDeltaPhiDeriv[n+i-1,n-i]/((n+i-1)!(n-i)!),{n,1,nMax},{i,1,n}]];
evenDerivs[nMax_] := Flatten[Table[withDeltaPhiDeriv[n+i-2,n-i]/((n+i-2)!(n-i)!),{n,1,nMax},{i,1,n}]];

Procedure to go from PMP to .xml file

(* from SDPB.m*)

(* A matrix with constant anti-diagonals given by the list bs *)
antiBandMatrix[bs_] := Module[
    {n = Ceiling[Length[bs]/2]},
    Reverse[Normal[
        SparseArray[
            Join[
                Table[Band[{i, 1}] -> bs[[n - i + 1]], {i, n}],
                Table[Band[{1, i}] -> bs[[n + i - 1]], {i, 2, n}]],
            {n, n}]]]];

(* DampedRationals are used to write the positive prefactor which is of the form c b^x / ((x-p1)(x-p2)...)*)
(* DampedRational[c, {p1, p2, ...}, b, x] stands for c b^x / ((x-p1)(x-p2)...) *)
(* It satisfies the following identities *)
DampedRational[const_, poles_, base_, x + a_] := 
    DampedRational[base^a const, # - a & /@ poles, base, x];
DampedRational[const_, poles_, base_, a_ /; FreeQ[a, x]] := 
    const base^a/Product[a - p, {p, poles}];
DampedRational/:x DampedRational[const_, poles_ /; MemberQ[poles, 0], base_, x] :=
    DampedRational[const, DeleteCases[poles, 0], base, x];

(* bilinearForm[f, m] = Integral[x^m f[x], {x, 0, Infinity}] *)
(* The special case when f[x] has no poles *)
bilinearForm[DampedRational[const_, {}, base_, x], m_] :=
    const Gamma[1+m] (-Log[base])^(-1-m);

(* The general DampedRational case for double and single poles*)
bilinearForm[DampedRational[c_, poles_, b_, x_], m_] := Module[
    {
        gatheredPoles = Gather[poles],
        quotientCoeffs = CoefficientList[PolynomialQuotient[x^m, Product[x-p, {p, poles}], x], x],
        integral, p, rest
    },
    integral[a_,1] := b^a Gamma[0, a Log[b]];
    integral[a_,2] := -1/a + b^a Gamma[0, a Log[b]] Log[b];
    c (Sum[
        p = gatheredPoles[[n,1]];
        rest = x^m / Product[x-q, {q, Join@@Delete[gatheredPoles, n]}];
        Switch[Length[gatheredPoles[[n]]],
               1, integral[p,1] rest /. x->p,
               2, integral[p,2] rest + integral[p,1] D[rest, x] /. x->p],
        {n, Length[gatheredPoles]}] + 
       Sum[
           quotientCoeffs[[n+1]] Gamma[1+n] (-Log[b])^(-1-n),
           {n, 0, Length[quotientCoeffs]-1}])];

(* orthogonalPolynomials[f, n] is a set of polynomials with degree 0
through n which are orthogonal with respect to the measure f[x] dx *)
orthogonalPolynomials[const_ /; FreeQ[const, x], 0] := {1/Sqrt[const]};

orthogonalPolynomials[const_ /; FreeQ[const, x], degree_] := 
    error["can't get orthogonal polynomials of nonzero degree for constant measure"];

orthogonalPolynomials[DampedRational[const_, poles_, base_, x], degree_] := 
    Table[x^m, {m, 0, degree}] . Inverse[
        CholeskyDecomposition[
            antiBandMatrix[
                Table[bilinearForm[DampedRational[const, Select[poles, # < 0&], base, x], m],
                      {m, 0, 2 degree}]]]];

(* Preparing SDP for Export *)
rho = SetPrecision[3-2 Sqrt[2], precision];

rescaledLaguerreSamplePoints[n_] := Table[
    SetPrecision[π^2 (-1+4k)^2/(-64n Log[rho]), precision],
    {k,0,n-1}];

maxIndexBy[l_,f_] := SortBy[
    Transpose[{l,Range[Length[l]]}],
    -f[First[#]]&][[1,2]];

(* finds v' such that a . v = First[v'] + a' . Rest[v'] when normalization . a == 1, where a' is a vector of length one less than a *)
reshuffleWithNormalization[normalization_, v_] := Module[
    {j = maxIndexBy[normalization, Abs], const},
    const = v[[j]]/normalization[[j]];
    Prepend[Delete[v - normalization*const, j], const]];

(* XML Exporting *)
nf[x_Integer] := x;
nf[x_] := NumberForm[SetPrecision[x,precision],precision,ExponentFunction->(Null&)];

safeCoefficientList[p_, x_] := Module[
    {coeffs = CoefficientList[p, x]},
    If[Length[coeffs] > 0, coeffs, {0}]];

WriteBootstrapSDP[file_, SDP[objective_, normalization_, positiveMatricesWithPrefactors_]] := Module[
    {
        stream = OpenWrite[file],
        node, real, int, vector, polynomial,
        polynomialVector, polynomialVectorMatrix,
        affineObjective, polynomialVectorMatrices
    },
	

    (* write a single XML node to file.  children is a routine that writes child nodes when run. *)
    node[name_, children_] := (
        WriteString[stream, "<", name, ">"];
        children[];
        WriteString[stream, "</", name, ">\n"];
    );

    real[r_][] := WriteString[stream, nf[r]];
    int[i_][] := WriteString[stream, i];
    vector[v_][] := Do[node["elt", real[c]], {c, v}];
    polynomial[p_][] := Do[node["coeff", real[c]], {c, safeCoefficientList[p,x]}];
    polynomialVector[v_][] := Do[node["polynomial", polynomial[p]], {p, v}];

    polynomialVectorMatrix[PositiveMatrixWithPrefactor[prefactor_, m_]][] := Module[
        {degree = Max[Exponent[m, x]], samplePoints, sampleScalings, bilinearBasis},
        samplePoints   = rescaledLaguerreSamplePoints[degree + 1];
		sampleScalings = Table[prefactor /. x -> a, {a, samplePoints}];
        bilinearBasis  = orthogonalPolynomials[prefactor, Floor[degree/2]];
        node["rows", int[Length[m]]];
        node["cols", int[Length[First[m]]]];
        node["elements", Function[
            {},
            Do[node[
                "polynomialVector",
                polynomialVector[reshuffleWithNormalization[normalization,pv]]],
               {row, m}, {pv, row}]]];
        node["samplePoints", vector[samplePoints]];
        node["sampleScalings", vector[sampleScalings]];
        node["bilinearBasis", polynomialVector[bilinearBasis]];
    ];                                                                                       
    node["sdp", Function[
        {},
        node["objective", vector[reshuffleWithNormalization[normalization, objective]]];
        node["polynomialVectorMatrices", Function[
            {},
            Do[node["polynomialVectorMatrix", polynomialVectorMatrix[pvm]];
	        , {pvm, positiveMatricesWithPrefactors}];
        ]];
    ]];                                          
    Close[stream];
];

Ising single correlator

Gap Δϵ

(* Function used to make an XML file for the single correlator with a gap deltaEpsLowPrec. *)
(* Has hardcoded dimensions d.*)
IsingGapEpsilon[datfile_][
    deltaPhiLowPrec_,
    deltaEpsLowPrec_,
    nmax_,
    keptPoleOrder_,
    Ls_,
    order_] := Module[
    {
		unitVec, pols,norm,obj,
		deltaPhi=SetPrecision[deltaPhiLowPrec,precision],
		deltaEps = SetPrecision[deltaEpsLowPrec,precision],
        addDeltaPhi,
        V,SpinZeroEps,SpinZeroPhi,SpinZeroT,
        F
    },
	(* adds DeltaPhi dependent prefactor.*)
	addDeltaPhi[x_] := (addDeltaPhi[x] = withDeltaPhiToWithoutDeltaPhi[x][2 nmax,2 nmax]);
  
	(*This function loads tables of derivatives of blocks as polynomials in x=Δ-j-1, and turns them into derivatives of F using addDeltaPhi, which involves derivatives of the v^Subscript[Δ, ϕ] terms.*)
	F[derivs_, deltaI_, deltaJ_, deltaK_, deltaL_, L_] := Module[                                                                                                                                                                    
        {delta12, delta34, derivReplacementRules},
        {delta12, delta34}= {deltaI-deltaJ,deltaK-deltaL};                                                                                                                                                                 
        derivReplacementRules = derivPolTable[3, delta12, delta34, L, nmax, keptPoleOrder, order];                                                                                                           
        (derivs[nmax] /. addDeltaPhi[(deltaJ + deltaK)/2] // Expand) /. derivReplacementRules // Expand                                                                                                                              
    ]; 
    
    
    (* The location of the poles in Δ in a conformal block, where order specifies how many poles *)
	cbPoles[nu_, L_, kkeptPoleOrder_] := Join[
	    Table[-(n + L - 1), {n, kkeptPoleOrder}],
	    Table[-(-nu + k - 1), {k, kkeptPoleOrder/2}],
	    Table[-(-L - 2 nu + n - 1), {n, Min[kkeptPoleOrder,L]}]];
    
    (*This is the positive prefactor χ that we multiply by the positive matrix so the blocks will have a natural normalization.*)
	cbPrefactor[delta_, L_, kkkeptPoleOrder_] := DampedRational[1, cbPoles[(3 - 2)/2, L, kkkeptPoleOrder], rho, delta];
    

	(*This function is where you put the crossing equations, written as tables of derivatives of F, 
		for every even spin, the positive prefactor χ is used just for stability of the numerics. 
		If there are more crossing equations, they enter as more lines in the Join function.*)

	V = Table[
			PositiveMatrixWithPrefactor[
				cbPrefactor[x+L+1,L,keptPoleOrder],
				Transpose[
					Join[
						Map[{{#}}&,    F[oddDerivs,  deltaPhi, deltaPhi, deltaPhi, deltaPhi, L]]
                    (* add additional crossing equations here.*)
                ],
                {3,1,2}
            ] 
        ],
        {L, Select[Ls,EvenQ]}
    ];

	(*x=0 defines the unitarity bound for j>0. Here we shift our assumption of the j=0 lower bound to a value deltaEps.*)
    pols = Join[
         MapAt[#/.x -> x+deltaEps-(3-2)  &, V, 1]
    ];

	(*The unit operator has unit OPE coefficient and unit blocks, so the derivatives of F in this case only involve derivatives of (v^Subscript[Δ, ϕ]).*)
	 unitVec = Join[
        oddDerivs[nmax]/.addDeltaPhi[deltaPhi]
    ]/.{zzbDeriv[0,0]->1,zzbDeriv[_,_]->0};


	(*Here we set the normalization of the functional. For scalar bounds we normalize the unit operator.*)
	norm = unitVec;
    
	(*Here we set the value of the objective function that we maximize. For scalar bounds the objective function is trivially a vector of zeros.*)
	obj = 0*norm;
    
    (*This calls the code that will write the XML file*)
	WriteBootstrapSDP[datfile, SDP[obj, norm, pols]];
];

Navigator for gap Δϵ

(* Function used to make XML file for a navigator computation. *)

IsingNavigator[datfile_][
	deltaPhiLowPrec_,
    deltaEpsLowPrec_,
    nmax_,
    keptPoleOrder_,
    Ls_,
    order_] := Module[
    {
		unitVec, pols,norm,obj,
		deltaPhi=SetPrecision[deltaPhiLowPrec,precision],
        deltaEps = SetPrecision[deltaEpsLowPrec,precision],
        deltaGFF,gffOP,
        addDeltaPhi,
        V,SpinZeroEps,SpinZeroPhi,SpinZeroT,
        F
    },
    deltaGFF= 2*deltaPhi;
    
	addDeltaPhi[x_] := (addDeltaPhi[x] = withDeltaPhiToWithoutDeltaPhi[x][2 nmax,2 nmax]);
  
	(*This function loads tables of derivatives of blocks as polynomials in x=Δ-j-1, and turns them into derivatives of F using addDeltaPhi, 
		which involves derivatives of the v^Subscript[Δ, ϕ] terms.*)

	F[derivs_, deltaI_, deltaJ_, deltaK_, deltaL_, L_] := Module[                                                                                                                                                                    
        {delta12, delta34, derivReplacementRules},
        {delta12, delta34}= {deltaI-deltaJ,deltaK-deltaL};                                                                                                                                                                 
        derivReplacementRules = derivPolTable[3, delta12, delta34, L, nmax, keptPoleOrder, order];                                                                                                           
        (derivs[nmax] /. addDeltaPhi[(deltaJ + deltaK)/2] // Expand) /. derivReplacementRules // Expand                                                                                                                              
    ]; 
    
    
    (* The location of the poles in Δ in a conformal block, where order specifies how many poles *)
    (* This has to match the rational approximation that was used to compute the blocks.*)
	cbPoles[nu_, L_, kkeptPoleOrder_] := Join[
		Table[-(n + L - 1), {n, kkeptPoleOrder}],
		Table[-(-nu + k - 1), {k, kkeptPoleOrder/2}],
		Table[-(-L - 2 nu + n - 1), {n, Min[kkeptPoleOrder,L]}]
	];
	
    (*This is the positive prefactor χ that we multiply by the positive matrix so the blocks will have a natural normalization.*)
    cbPrefactor[delta_, L_, kkkeptPoleOrder_] := DampedRational[1, cbPoles[(3 - 2)/2, L, kkkeptPoleOrder], rho, delta];
    

	(*This function is where you put the crossing equations, written as tables of derivatives of F, 
		for every even spin, the positive prefactor χ is used just for stability of the numerics. 
		If there are more crossing equations, they enter as more lines in the Join function.*)

	V = Table[
			PositiveMatrixWithPrefactor[
				cbPrefactor[x+L+1,L,keptPoleOrder],
				Transpose[
					Join[
						Map[{{#}}&,    F[oddDerivs,  deltaPhi, deltaPhi, deltaPhi, deltaPhi, L]]
						],
                {3,1,2}] 
			],
		{L, Select[Ls,EvenQ]}];
    
	(* In order to guarantee boundedness of the correlator we normalize on the GFF-operator that is excluded by our gap assumption:
	Note that this input is just a number because we fixed x, so no need to remove poles and r^Δ to make it a polynomial.*)
    
	gffOP={PositiveMatrixWithPrefactor[DampedRational[1,{},1,1],
			Transpose[
				((rho)^(x+3-2)/((x+3-2-cbPoles[(3-2)/2,0,keptPoleOrder])/.List->Times))   Join[
					Map[{{#}}&,    F[oddDerivs,  deltaPhi, deltaPhi, deltaPhi, deltaPhi, 0]]
                 ]/.x-> deltaGFF-(3-2)
			,{3,1,2}]]
	};

	(*x=0 corresponds to the unitarity bound for j>0. For the scalar we shift so that x=0 corresponds to the j=0 lower bound value deltaEps.*)
    pols = Join[
         MapAt[#/.x -> x+deltaEps-(3-2)  &, V, 1]
    ];

	(*The unit operator has unit OPE coefficient and unit blocks, so the derivatives of F in this case only involve derivatives of (v^Subscript[Δ, ϕ]).*)
	unitVec = Join[
			oddDerivs[nmax]/.addDeltaPhi[deltaPhi]
			]/.{zzbDeriv[0,0]->1,zzbDeriv[_,_]->0};
	
	(* In this case we normalize α.VGFF=-1*)
	norm=-Flatten[gffOP[[1,2]]];
	
	(*Here we set the value of the objective function that we maximize. For scalar bounds the objective function is trivially a vector of zeros.*)
    
    obj = unitVec;
    
    (*This calls the code that will write the XML file using *)
    
	WriteBootstrapSDP[datfile, SDP[obj, norm, pols]];
];

One relevant scalar

(* Function used to make an XML file for the single correlator with a gap deltaEpsLowPrec. *)
(* Has hardcoded dimensions d.*)
IsingOneRelevantScalar[datfile_][
    deltaPhiLowPrec_,
    deltaEpsLowPrec_,
    nmax_,
    keptPoleOrder_,
    Ls_,
    order_] := Module[
    {
		unitVec, pols,norm,obj,
		deltaPhi=SetPrecision[deltaPhiLowPrec,precision],
		deltaEps = SetPrecision[deltaEpsLowPrec,precision],
        addDeltaPhi,fixedEpsVec,
        V,SpinZeroEps,SpinZeroPhi,SpinZeroT,
        F,
        gapEpsilonPrime=3
    },
	(* adds DeltaPhi dependent prefactor.*)
	addDeltaPhi[x_] := (addDeltaPhi[x] = withDeltaPhiToWithoutDeltaPhi[x][2 nmax,2 nmax]);
  
	(*This function loads tables of derivatives of blocks as polynomials in x=Δ-j-1, and turns them into derivatives of F using addDeltaPhi, which involves derivatives of the v^Subscript[Δ, ϕ] terms.*)
	F[derivs_, deltaI_, deltaJ_, deltaK_, deltaL_, L_] := Module[                                                                                                                                                                    
        {delta12, delta34, derivReplacementRules},
        {delta12, delta34}= {deltaI-deltaJ,deltaK-deltaL};                                                                                                                                                                 
        derivReplacementRules = derivPolTable[3, delta12, delta34, L, nmax, keptPoleOrder, order];                                                                                                           
        (derivs[nmax] /. addDeltaPhi[(deltaJ + deltaK)/2] // Expand) /. derivReplacementRules // Expand                                                                                                                              
    ]; 
    
    
    (* The location of the poles in Δ in a conformal block, where order specifies how many poles *)
	cbPoles[nu_, L_, kkeptPoleOrder_] := Join[
	    Table[-(n + L - 1), {n, kkeptPoleOrder}],
	    Table[-(-nu + k - 1), {k, kkeptPoleOrder/2}],
	    Table[-(-L - 2 nu + n - 1), {n, Min[kkeptPoleOrder,L]}]];
    
    (*This is the positive prefactor χ that we multiply by the positive matrix so the blocks will have a natural normalization.*)
	cbPrefactor[delta_, L_, kkkeptPoleOrder_] := DampedRational[1, cbPoles[(3 - 2)/2, L, kkkeptPoleOrder], rho, delta];
    

	(*This function is where you put the crossing equations, written as tables of derivatives of F, 
		for every even spin, the positive prefactor χ is used just for stability of the numerics. 
		If there are more crossing equations, they enter as more lines in the Join function.*)

	V = Table[
			PositiveMatrixWithPrefactor[
				cbPrefactor[x+L+1,L,keptPoleOrder],
				Transpose[
					Join[
						Map[{{#}}&,    F[oddDerivs,  deltaPhi, deltaPhi, deltaPhi, deltaPhi, L]]
                    (* add additional crossing equations here.*)
                ],
                {3,1,2}
            ] 
        ],
        {L, Select[Ls,EvenQ]}
    ];

	(* fixed x=Δϵ operator.*)
	fixedEpsVec =  {PositiveMatrixWithPrefactor[DampedRational[1,{},1,1],
		Transpose[
			((rho)^(x+3-2)/((x+3-2-cbPoles[(3-2)/2,0,keptPoleOrder])/.List->Times))   Join[
                    Map[{{#}}&,    F[oddDerivs,  deltaPhi, deltaPhi, deltaPhi, deltaPhi, 0]]
                 ]/.x-> deltaEps-(3-2)
		,{3,1,2}]]
	};
	
	(*x=0 corresponds to the unitarity bound for j>0. For teh scalar we shift so that x=0 corresponds to the j=0 lower bound value deltaEps.*)
    pols = Join[
         MapAt[#/.x -> x+gapEpsilonPrime-(3-2)  &, V, 1],fixedEpsVec
    ];

	(*The unit operator has unit OPE coefficient and unit blocks, so the derivatives of F in this case only involve derivatives of (v^Subscript[Δ, ϕ]).*)
	 unitVec = Join[
        oddDerivs[nmax]/.addDeltaPhi[deltaPhi]
    ]/.{zzbDeriv[0,0]->1,zzbDeriv[_,_]->0};


	(*Here we set the normalization of the functional. For scalar bounds we normalize the unit operator.*)
	norm = unitVec;
    
	(*Here we set the value of the objective function that we maximize. For scalar bounds the objective function is trivially a vector of zeros.*)
	obj = 0*norm;
    
    (*This calls the code that will write the XML file*)
	WriteBootstrapSDP[datfile, SDP[obj, norm, pols]];
];

Navigator for one relevant scalar

(* Function used to make an XML file for the single correlator with a gap deltaEpsLowPrec. *)
(* Has hardcoded dimensions d.*)
IsingNavigatorOneRelevantScalar[datfile_][
    deltaPhiLowPrec_,
    deltaEpsLowPrec_,
    nmax_,
    keptPoleOrder_,
    Ls_,
    order_] := Module[
    {
		unitVec, pols,norm,obj,
		deltaPhi=SetPrecision[deltaPhiLowPrec,precision],
		deltaEps = SetPrecision[deltaEpsLowPrec,precision],
        addDeltaPhi,fixedEpsVec,
        V,SpinZeroEps,SpinZeroPhi,SpinZeroT,
        F,deltaGFF,gffOP,
        gapEpsilonPrime=3
    },
    
    deltaGFF= 2*deltaPhi;
	(* adds DeltaPhi dependent prefactor.*)
	addDeltaPhi[x_] := (addDeltaPhi[x] = withDeltaPhiToWithoutDeltaPhi[x][2 nmax,2 nmax]);
  
	(*This function loads tables of derivatives of blocks as polynomials in x=Δ-j-1, and turns them into derivatives of F using addDeltaPhi, which involves derivatives of the v^Subscript[Δ, ϕ] terms.*)
	F[derivs_, deltaI_, deltaJ_, deltaK_, deltaL_, L_] := Module[                                                                                                                                                                    
        {delta12, delta34, derivReplacementRules},
        {delta12, delta34}= {deltaI-deltaJ,deltaK-deltaL};                                                                                                                                                                 
        derivReplacementRules = derivPolTable[3, delta12, delta34, L, nmax, keptPoleOrder, order];                                                                                                           
        (derivs[nmax] /. addDeltaPhi[(deltaJ + deltaK)/2] // Expand) /. derivReplacementRules // Expand                                                                                                                              
    ]; 
    
    
    (* The location of the poles in Δ in a conformal block, where order specifies how many poles *)
    (* nu=((dim-2)/2), and is often used in equations like this because it simplifies them a bit.*)
	cbPoles[nu_, L_, kkeptPoleOrder_] := Join[
	    Table[-(n + L - 1), {n, kkeptPoleOrder}],
	    Table[-(-nu + k - 1), {k, kkeptPoleOrder/2}],
	    Table[-(-L - 2 nu + n - 1), {n, Min[kkeptPoleOrder,L]}]];
    
    (*This is the positive prefactor χ that we multiply by the positive matrix so the blocks will have a natural normalization.*)
	cbPrefactor[delta_, L_, kkkeptPoleOrder_] := DampedRational[1, cbPoles[(3 - 2)/2, L, kkkeptPoleOrder], rho, delta];
    

	(*This function is where you put the crossing equations, written as tables of derivatives of F, 
		for every even spin, the positive prefactor χ is used just for stability of the numerics. 
		If there are more crossing equations, they enter as more lines in the Join function.*)

	V = Table[
			PositiveMatrixWithPrefactor[
				cbPrefactor[x+L+1,L,keptPoleOrder],
				Transpose[
					Join[
						Map[{{#}}&,    F[oddDerivs,  deltaPhi, deltaPhi, deltaPhi, deltaPhi, L]]
                    (* add additional crossing equations here.*)
                ],
                {3,1,2}
            ] 
        ],
        {L, Select[Ls,EvenQ]}
    ];
	
	gffOP={PositiveMatrixWithPrefactor[DampedRational[1,{},1,1],
			Transpose[
				((rho)^(x+3-2)/((x+3-2-cbPoles[(3-2)/2,0,keptPoleOrder])/.List->Times))   Join[
					Map[{{#}}&,    F[oddDerivs,  deltaPhi, deltaPhi, deltaPhi, deltaPhi, 0]]
                 ]/.x-> deltaGFF-(3-2)
			,{3,1,2}]]
	};


	(* fixed x=Δϵ operator.*)
	fixedEpsVec =  {PositiveMatrixWithPrefactor[DampedRational[1,{},1,1],
		Transpose[
			((rho)^(x+3-2)/((x+3-2-cbPoles[(3-2)/2,0,keptPoleOrder])/.List->Times))   Join[
                    Map[{{#}}&,    F[oddDerivs,  deltaPhi, deltaPhi, deltaPhi, deltaPhi, 0]]
                 ]/.x-> deltaEps-(3-2)
		,{3,1,2}]]
	};
	
	(*x=0 corresponds to the unitarity bound for j>0. For teh scalar we shift so that x=0 corresponds to the j=0 lower bound value deltaEps.*)
    pols = Join[
         MapAt[#/.x -> x+gapEpsilonPrime-(3-2)  &, V, 1],fixedEpsVec
    ];

	(*The unit operator has unit OPE coefficient and unit blocks, so the derivatives of F in this case only involve derivatives of (v^Subscript[Δ, ϕ]).*)
	 unitVec = Join[
        oddDerivs[nmax]/.addDeltaPhi[deltaPhi]
    ]/.{zzbDeriv[0,0]->1,zzbDeriv[_,_]->0};


	(* In this case we normalize α.VGFF=-1*)
	norm=-Flatten[gffOP[[1,2]]];
    
	(*Here we set the value of the objective function that we maximize. For scalar bounds the objective function is trivially a vector of zeros.*)
	obj = unitVec;
    
    (*This calls the code that will write the XML file*)
	WriteBootstrapSDP[datfile, SDP[obj, norm, pols]];
];

Generating the blocks, writing to input files

(* MAC Version*)
(* generate the blocks via docker *)
GenerateBlocks[]:=Module[{options,command,out},
(* options="--output-poles";*)
options="";
command=dockerExecutable<>" run -v $PWD:/usr/local/share/sdpb wlandry/"<>sdpb<>" scalar_blocks --dim "<>ToString[dim] <>" --order "<>ToString[order]<> " --max-derivs "<>ToString[nmax]<>" --spin-ranges "<>spinrange <>" --poles "<>ToString[poles]<>" --delta-12 "<>ToString[delta12]<>"  --delta-34 "<>ToString[delta34]<>" --num-threads="<>ToString[procsPerNode]<>" -o /usr/local/share/sdpb/output --precision="<>ToString[precision]<>" "<>options<>"\n";
Export[FileNameJoin@{outputDataDirectory,"temp-script"},command,"Text"];
VPrint["####    running scalar-blocks"];
out=RunProcess["/bin/tcsh","StandardOutput","cd "<>outputDataDirectory<>"; rm -r output; set path=(/usr/local/share/ $path); bash temp-script; ls output/*L0*"] ;
If[StringFreeQ[out,"L0"],Print["failure."],Print["sucess, we generated (L0 only):"];Print[out]]
]

(* write xml file *)
WriteXMLfile[deltaphi_,deltaeps_,spectrumTest_:"GapEpsilon"]:=Module[{},
	VPrint["####    writing xml-file"];
	Switch[spectrumTest,
		"GapEpsilon",
			IsingGapEpsilon[FileNameJoin@{outputDataDirectory,"torun4sdpb.xml"}][deltaphi,deltaeps,nmax,poles,spinrangeMathematica,order],
		"NavigatorGapEpsilon",
			IsingNavigator[FileNameJoin@{outputDataDirectory,"torun4sdpb_experimental.xml"}][deltaphi,deltaeps,nmax,poles,spinrangeMathematica,order],
		"OneRelevantScalar", 
			IsingOneRelevantScalar[FileNameJoin@{outputDataDirectory,"torun4sdpb.xml"}][deltaphi,deltaeps,nmax,poles,spinrangeMathematica,order],
		"NavigatorOneRelevantScalar",
			IsingNavigatorOneRelevantScalar[FileNameJoin@{outputDataDirectory,"torun4sdpb.xml"}][deltaphi,deltaeps,nmax,poles,spinrangeMathematica,order],
		_, Print["Bug in WriteXMLfile."]
	];
    VPrint["    done."]
]

Here we assume a really simple workflow. You always want to perform the next step on the last input generated, i.e. convert the latest xml file the preprocessed SDP files, or running SDPB on the latest SDP files.

(* preprocess xml to input sdpb via docker *)
Clear[ConvertXml2SDPB]
ConvertXml2SDPB[xmlFileName_:"torun4sdpb.xml"]:=Module[{command,out},
command=dockerExecutable<>" run -v "<>outputDataDirectory<>":/usr/local/share/sdpb wlandry/"<>sdpb<>" mpirun --allow-run-as-root -n "<>ToString[procsPerNode]<>" pvm2sdp "<>ToString[precision] <>" /usr/local/share/sdpb/"<>xmlFileName<>" /usr/local/share/sdpb/input\n";
Export[FileNameJoin@{outputDataDirectory,"temp-script"},command,"Text"];
VPrint["####    running pvm2sdp"];
out=RunProcess["/bin/tcsh","ExitCode","cd "<>outputDataDirectory<>"; set path=(/usr/local/share/ $path); bash temp-script;"] ;
If[out>0,Print["    error."],VPrint["    done."]];
out
]

Running sdpb

(* Mac Version*)
(* We remove the timings each time. This is not that efficient but avoids errors when you change the setup.*)
RunSDPB[paramFile_:"param.sdpb"]:=Module[{command,out,outFileName,bashLauncher},
(* writes progress to sdpb_progress.txt, however I don't know whether tee works on MAC.*)
If[$SystemID==="Linux-x86-64",
command="cd "<>outputDataDirectory<>"; "<>dockerExecutable<>" run -v $PWD:/usr/local/share/sdpb wlandry/"<>sdpb<>" mpirun --allow-run-as-root -n "<>ToString[procsPerNode]<>" sdpb -p /usr/local/share/sdpb/"<>paramFile<>" -s /usr/local/share/sdpb/input 2>&1 | tee sdpb_progress.txt;\n"
,
(* without tee for MAC*)
command="cd "<>outputDataDirectory<>"; "<>dockerExecutable<>" run -v $PWD:/usr/local/share/sdpb wlandry/"<>sdpb<>" mpirun --allow-run-as-root -n "<>ToString[procsPerNode]<>" sdpb -p /usr/local/share/sdpb/"<>paramFile<>" -s /usr/local/share/sdpb/input ;\n"
];
Export[FileNameJoin@{outputDataDirectory,"temp-script.command"},command,"Text"];
outFileName=FileNameJoin@{outputDataDirectory,"/input_out/out.txt"};
Quiet[DeleteFile[outFileName]];
VPrint["####    running sdpb"];
bashLauncher=If[$SystemID==="Linux-x86-64","bash","open"];
out=RunProcess["/bin/sh","StandardOutput","cd "<>outputDataDirectory<>"; chmod u+x temp-script.command; "<>bashLauncher<>" temp-script.command;"] ;
If[out>0,Print["    error."],Print["    done."]];
Print["... waiting for sdpb"];
While[!FileExistsQ[outFileName ]|| StringLength[(out=Import[outFileName])]<1,Pause[1]];
Print["... sdpb is done."];
out=Import[outFileName]; 
out
]

(*run[command_String]:=RunProcess[{"bash", "-c", command}]*)

Clear[feasibilityResults,CheckFeasible];

feasibilityResults=<||>;

CheckFeasible[DeltaPhi_,DeltaEpsilon_]:=Module[{feasibility,result,primalFeasible,dualFeasible,test},
test=feasibilityResults[{DeltaPhi,DeltaEpsilon}];
If[test==True||test==False,
Print["... known as ",test];Return[test]];
WriteXMLfile[DeltaPhi,DeltaEpsilon];
ConvertXml2SDPB[];
result=RunSDPB[];
VPrint[result];
feasibility=Which[
	StringContainsQ[result,"found primal feasible solution"], "Primal",
	StringContainsQ[result,"found dual feasible solution"], "Dual",
	StringContainsQ[result,"found primal-dual optimal solution"],"PrimalDualOptimal",
	True,Print["Error: SDPB did not converge for Δϕ=",DeltaPhi,", Δϵ=",DeltaEpsilon,". Aborting..."];"Error"];
feasibilityResults[{DeltaPhi,DeltaEpsilon}]=feasibility;
feasibility
];

Clear[Bisect];
Bisect[Δphi_,startlowerbound_,startupperbound_,err_:1/10000]:=Module[{X,i,lowerbound,upperbound,Δguess,res,feasibility},
lowerbound=startlowerbound;
upperbound=startupperbound;
(* Check that the lowerbound is indeed primal and that upperbound is indeed Dual:*)
(*If[
	And[CheckFeasible[Δphi,lowerbound]=="Primal",CheckFeasible[Δphi,upperbound]=="Dual"],
		Return["error"]
];*)
While[upperbound-lowerbound>err,
		Δguess=((lowerbound+upperbound)/2)//N//N;
		feasibility=CheckFeasible[Δphi,Δguess];
		VPrint[feasibility," at ",Δguess];
		Switch[feasibility,
			"Primal",
				lowerbound=Δguess,
			"Dual",
				upperbound=Δguess,
			"PrimalDualOptimal",
				Print["Warning bound ", {Δphi,Δguess}, " is primal-dual optimal. Returning the guess."];Return[Δguess],
			_,
			Print["Error: SDPB did not converge for ",{Δphi,Δguess}, " Aborting..."];Abort[]
		];
	];
Δguess//N
]

Clear[navigatorResults,ComputeNavigator];

navigatorResults=<||>;

ComputeNavigator[DeltaPhi_,DeltaEpsilon_]:=Module[{result,primalFeasible,dualFeasible,test,navigator},
	test=navigatorResults[{DeltaPhi,DeltaEpsilon}];
	If[test==True||test==False,
		Print["... known as ",test];Return[test]];
	IsingNavigator[FileNameJoin@{outputDataDirectory,"torun4sdpb_experimental.xml"}][DeltaPhi,DeltaEpsilon,nmax,poles,spinrangeMathematica,order];
	ConvertXml2SDPB["torun4sdpb_experimental.xml"];
	result=RunSDPB["param_optimization.sdpb"];
	VPrint[result];
	navigator=StringCases[result,Shortest["primalObjective ="~~obj___~~";"]:>ToExpression[sdpbNumStringToMMANumber[obj]]][[1]];
	navigatorResults[{DeltaPhi,DeltaEpsilon}]=navigator;
	navigator
];

Clear[feasibilityResultsOneRelevantScalar,CheckFeasibleOneRelevantScalar];

feasibilityResultsOneRelevantScalar=<||>;

CheckFeasibleOneRelevantScalar[DeltaPhi_,DeltaEpsilon_]:=Module[{feasibility,result,primalFeasible,dualFeasible,test},
test=feasibilityResultsOneRelevantScalar[{DeltaPhi,DeltaEpsilon}];
If[test==True||test==False,
Print["... known as ",test];Return[test]];
WriteXMLfile[DeltaPhi,DeltaEpsilon,"OneRelevantScalar"];
ConvertXml2SDPB[];
result=RunSDPB[];
VPrint[result];
feasibility=Which[
	StringContainsQ[result,"found primal feasible solution"], "Primal",
	StringContainsQ[result,"found dual feasible solution"], "Dual",
	StringContainsQ[result,"found primal-dual optimal solution"],"PrimalDualOptimal",
	True,Print["Error: SDPB did not converge for Δϕ=",DeltaPhi,", Δϵ=",DeltaEpsilon,". Aborting..."];"Error"];
feasibilityResultsOneRelevantScalar[{DeltaPhi,DeltaEpsilon}]=feasibility;
feasibility
];

Clear[navigatorResultsOneRelevantScalar,ComputeNavigatorOneRelevantScalar];

navigatorResultsOneRelevantScalar=<||>;

ComputeNavigatorOneRelevantScalar[DeltaPhi_,DeltaEpsilon_]:=Module[{result,primalFeasible,dualFeasible,test,navigator},
test=navigatorResultsOneRelevantScalar[{DeltaPhi,DeltaEpsilon}];
If[test==True||test==False,
Print["... known as ",test];Return[test]];
WriteXMLfile[DeltaPhi,DeltaEpsilon,"NavigatorOneRelevantScalar"];
ConvertXml2SDPB[];
result=RunSDPB["param_optimization.sdpb"];
VPrint[result];
navigator=StringCases[result,Shortest["primalObjective ="~~obj___~~";"]:>ToExpression[sdpbNumStringToMMANumber[obj]]][[1]];
	navigatorResultsOneRelevantScalar[{DeltaPhi,DeltaEpsilon}]=navigator;
	navigator
];

sdpbNumStringToMMANumber[string_]:=StringReplace[string,{"e+":>"*^","e-":>"*^-"}]

Of course one of the major advantages of the navigator function is that you can compute its gradient for free. Unfortunately the tool to do that is not yet ready for release although it will be soon. I do not feel it is the best idea to share an experimental version just before the release of the stable one. So the computation of the gradient is not covered in the tutorial.

Main

Comment: The blocks in scalar_block are not ideal (I believe the prefactor contains a pole not present in the original which causes this [as found by Ning Su]. I am not getting any dual jumps. Everything still works but now convergence is much slower. You will notice that you never get dual-jumps. Dual and primal jumps are moments in the SDPB optimization where either the primal error or dual error sharply reduces in one iteration step. When this happens the task should eventually become primal or dual respectively and you can immediately terminate. (Also hot-starting will probably not work properly.)

First you will need to generate the blocks. It should compute them and then print the available ones.

GenerateBlocks[]

Examples Δϵ gap

verbatim=True;
Δϕ=.52;
Δϵ=1.42;
CheckFeasible[Δϕ,Δϵ]

verbatim=True;
Bisect[0.52,1.4,1.5,0.01]

feasibilityResults

Examples Navigator

verbatim=True;
ComputeNavigator[0.52,1.45]

verbatim=True;
ComputeNavigator[0.52,1.4]

navigatorResults//N

Example One Relevant Scalar

Please delete the directory input.ck when changing setup. Unfortunately I could not get docker to play nicely with file permissions. It wants to do everything as root and create files that can only be modified by root. So I cannot delete it from Mathematica.

CheckFeasibleOneRelevantScalar[0.52,1.4]

CheckFeasibleOneRelevantScalar[0.52,0.7]

CheckFeasibleOneRelevantScalar[0.52,1.3]

CheckFeasibleOneRelevantScalar[0.52,1.55]

feasibilityResultsOneRelevantScalar

Example One Relevant Scalar Navigator

ComputeNavigatorOneRelevantScalar[0.52,1.5]

ComputeNavigatorOneRelevantScalar[0.52,0.7]

ComputeNavigatorOneRelevantScalar[0.52,1.3]

navigatorResults

nav052=ComputeNavigatorOneRelevantScalar[0.52,#]&/@Range[0.65,1.85,0.1]

{#1[[2]],#2}&@@@Normal[navigatorResultsOneRelevantScalar]//ListPlot

Negative Navigator means that a point is allowed. Positive means it is excluded. We see the results match the results above.