[sdpb] 147/233: Added Bootstrap2dExample.m
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:31 UTC 2017
This is an automated email from the git hooks/post-receive script.
thansen pushed a commit to branch master
in repository sdpb.
commit 1f7e30066cd27a04b8f08a9a4e44c0d11dd8b9ff
Author: David Simmons-Duffin <dsd at ssh1.sns.ias.edu>
Date: Thu Feb 5 18:18:04 2015 -0500
Added Bootstrap2dExample.m
---
mathematica/Bootstrap2dExample.m | 163 +++++++++++++++++++++++++++++++++++++++
mathematica/SDPB.m | 35 +++++++--
2 files changed, 191 insertions(+), 7 deletions(-)
diff --git a/mathematica/Bootstrap2dExample.m b/mathematica/Bootstrap2dExample.m
new file mode 100644
index 0000000..ffd7692
--- /dev/null
+++ b/mathematica/Bootstrap2dExample.m
@@ -0,0 +1,163 @@
+(* ::Package:: *)
+
+<<"SDPB.m";
+
+half = SetPrecision[1/2,prec];
+
+rho[z_] := z/(1+Sqrt[1-z])^2;
+
+(* An SL_2 conformal block in the rho coordinate *)
+chiralBlock[x_, rho_] := rho^(x/2) Hypergeometric2F1[1/2,x/2,(x+1)/2,rho^2];
+
+(* Turn a series a_0 + a_1 x + ... into a list of rules { name[0] ->
+ 0! a_0, name[1] -> 1! a_1, ... } *)
+seriesToRules[name_, series_] := Table[
+ (* series[[3]] is the list of coefficients *)
+ name[n] -> n! series[[3, n+1]],
+ {n, 0, Length[series[[3]]] - 1}];
+
+(* A list of replacement rules zDeriv[n] -> sum of rhoDeriv[k]'s
+giving the derivatives of a function with respect to z around 1/2 in
+terms of derivatives with respect to rho around rho[1/2] *)
+zDerivTable[order_] :=
+ (zDerivTable[order] =
+ Module[{dRhoSeries = Series[rho[half+dz] - rho[half], {dz, 0, order}]},
+ seriesToRules[zDeriv, Sum[rhoDeriv[n] dRhoSeries^n / n!, {n, 0, order}]]]);
+
+(* A table of the form { prefactor, { zDeriv[n] -> polynomial_n(x),
+... } }, where prefactor * polynomial_n(x) approximates the n-th
+derivative of a chiral conformal bock around z=1/2 *)
+chiralBlockTable[derivativeOrder_, keptPoleOrder_] :=
+ (chiralBlockTable[derivativeOrder, keptPoleOrder] =
+ Module[
+ {
+ approx,
+ numerator,
+ prefactor,
+ poles,
+ derivTable
+ },
+ poles = Table[n, {n, 1, keptPoleOrder-1, 2}];
+ numerator = (Series[chiralBlock[x,rho], {rho, 0, keptPoleOrder}] Product[x+n, {n, poles}]) // Normal;
+ prefactor = DampedRational[1,Table[-n,{n,poles}],rhoCrossing^(1/2),x];
+ (* We compute derivatives with respect to rho and use
+ zDerivTable to convert to zDerivatives *)
+ derivTable = zDerivTable[derivativeOrder] /. Table[
+ rhoDeriv[n] -> ((rho^(-x/2) D[numerator, {rho,n}] // Expand) /. rho->rhoCrossing),
+ {n, 0, derivativeOrder}];
+ {prefactor,Expand[derivTable]}]);
+
+(* Replacement rules giving derivatives of (1-z)^deltaPhi f[z] around
+z=1/2 in terms of derivatives of f[z] *)
+withDeltaPhiDerivTable[deltaPhi_,order_] :=
+ seriesToRules[
+ withDeltaPhiDeriv,
+ Series[(half-dz)^deltaPhi f[half+dz], {dz, 0, order}] /. {
+ Derivative[j_][f][_] :> zDeriv[j],
+ f[_] :> zDeriv[0]}];
+
+(* A list of derivative such that m >= n, m+n <= derivativeOrder, and
+m+n is odd *)
+oddDerivs[derivativeOrder_] :=
+ Flatten[Table[zzbDeriv[m,n],
+ {m,0,derivativeOrder},
+ {n, 1 - Mod[m,2], Min[m, derivativeOrder-m], 2}]
+ ,1];
+
+(* Test whether the point (deltaPhiLowPrecision,
+deltaPhi2LowPrecision) is allowed in a Z_2-symmetric CFT.
+derivativeOrder gives the number of derivatives to use, keptPoleOrder
+controls the accuracy of the conformal block approximation, and Lmax
+sets the number of included spins *)
+singletAllowed2d[deltaPhiLowPrecision_, deltaPhi2LowPrecision_, derivativeOrder_, keptPoleOrder_, Lmax_] :=
+ Module[
+ {
+ deltaPhi = SetPrecision[deltaPhiLowPrecision,prec],
+ deltaPhi2 = SetPrecision[deltaPhi2LowPrecision,prec],
+ chiralBlocksPrefactor,
+ chiralBlockPols,
+ chiralBlocksWithDeltaPhi,
+ unitOp,
+ pols,
+ norm,
+ obj
+ },
+ {chiralBlocksPrefactor, chiralBlockPols} = chiralBlockTable[derivativeOrder, keptPoleOrder];
+ chiralBlocksWithDeltaPhi = withDeltaPhiDerivTable[deltaPhi, derivativeOrder] /. chiralBlockPols;
+ pols=Table[
+ PositiveMatrixWithPrefactor[
+ (chiralBlocksPrefactor/.x->x+2L) (chiralBlocksPrefactor),
+ (* These are 1x1 matrices of polynomial vectors *)
+ {{oddDerivs[derivativeOrder]/.{
+ zzbDeriv[m_,n_] :> (
+ (withDeltaPhiDeriv[m]/.chiralBlocksWithDeltaPhi /. x -> x+2L)
+ (withDeltaPhiDeriv[n]/.chiralBlocksWithDeltaPhi) +
+ (withDeltaPhiDeriv[n]/.chiralBlocksWithDeltaPhi /. x -> x+2L)
+ (withDeltaPhiDeriv[m]/.chiralBlocksWithDeltaPhi))
+ }//Expand}}],
+ {L, 0, Lmax, 2}];
+ (* Replace x -> x + deltaPhi2 for scalar operators *)
+ pols = MapAt[# /. x -> x + deltaPhi2 &, pols, 1];
+ unitOp = oddDerivs[derivativeOrder] /. {
+ zzbDeriv[m_,n_] :> 2 withDeltaPhiDeriv[m] withDeltaPhiDeriv[n]
+ } /. withDeltaPhiDerivTable[deltaPhi, derivativeOrder] /. {
+ zDeriv[0] :> 1, zDeriv[_] :> 0};
+ norm = unitOp;
+ obj = 0*norm;
+ SolveBootstrapSDP[SDP[obj,norm,pols]]];
+
+(* This is not a recommended long-term solution for evaluating SDPs
+for the following reasons: 1) different sdpFiles should have unique
+names so that they can be sovled in parallel and so that their
+checkpoints and output files don't overwrite each other 2) The
+Run[...] command forces Mathematica to be running until sdpb
+finishes. It is better to use WriteBootstrapSDP instead of
+SolveBootstrapSDP and run sdpb by hand or with an external script. *)
+SolveBootstrapSDP[sdp_] := Module[
+ {
+ sdpFile = "mySDP.xml",
+ outFile = "mySDP.out"
+ },
+ WriteBootstrapSDP[sdpFile, sdp];
+ (* Most of the defaults are way over the top for this size
+ problem, but we'll use them because it's easy. If you want speed,
+ try fiddling with some of the parameters. *)
+ Run["sdpb", "-s", sdpFile,
+ "--findPrimalFeasible",
+ "--findDualFeasible",
+ "--noFinalCheckpoint"];
+ (* Careful! Simply evaluating the output file will bring a bunch
+ of variables into scope! *)
+ Get[outFile];
+ (* If we forget to clear x, it'll ruin everything *)
+ Clear[x,y];
+ Switch[
+ terminateReason,
+ "found primal feasible solution", True,
+ "found dual feasible solution", False,
+ _, Throw[terminateReason]]];
+
+(* For a binary-valued function f, find the value of x in the interval
+(true, false) where f changes from True to False, to within an error
+of thresh. *)
+binarySearch[f_, true_, false_, thresh_] := Module[
+ {test = (true + false)/2},
+ If[Abs[true - false] < thresh,
+ test,
+ WriteString["stdout", "> trying: ", test, "\n"];
+ If[f[test],
+ binarySearch[f, test, false, thresh],
+ binarySearch[f, true, test, thresh]]]];
+
+(* Compute an approximate upper bound on deltaPhi^2, as a function of
+deltaPhi. This function takes a few minutes to run *)
+bootstrapBound2d[] = Table[
+ WriteString["stdout", "> deltaPhi = ", deltaPhi, "\n"];
+ {deltaPhi,
+ binarySearch[
+ singletAllowed2d[deltaPhi, #, 7, 10, 15] &,
+ 0.1,
+ 2,
+ 0.02]},
+ {deltaPhi, 0.01,0.25,0.01}];
+
diff --git a/mathematica/SDPB.m b/mathematica/SDPB.m
index a25a094..92f0c93 100644
--- a/mathematica/SDPB.m
+++ b/mathematica/SDPB.m
@@ -25,20 +25,42 @@ DampedRational[const_, poles_, base_, a_ /; FreeQ[a, x]] :=
DampedRational/:x DampedRational[const_, poles_ /; MemberQ[poles, 0], base_, x] :=
DampedRational[const, DeleteCases[poles, 0], base, x];
+DampedRational/:DampedRational[c1_,p1_,b1_,x] DampedRational[c2_,p2_,b2_,x] :=
+ DampedRational[c1 c2, Join[p1, p2], b1 b2, 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);
-memoizeGamma[a_,b_]:=memoizeGamma[a,b]=Gamma[a,b];
+(*memoizeGamma[a_,b_]:=memoizeGamma[a,b]=Gamma[a,b];*)
-(* The general DampedRational case *)
-bilinearForm[DampedRational[const_, poles_, base_, x], m_] :=
+(* The case where f[x] has only single poles *)
+(*bilinearForm[DampedRational[const_, poles_, base_, x], m_] :=
const Sum[
((-poles[[i]])^m) ( base^poles[[i]]) Gamma[1 + m] memoizeGamma[-m, poles[[i]] Log[base]]/
Product[poles[[i]] - p, {p, Delete[poles, i]}],
- {i, Length[poles]}];
+ {i, Length[poles]}];*)
+
+(* The case where f[x] can have single or double 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 *)
@@ -55,10 +77,10 @@ orthogonalPolynomials[DampedRational[const_, poles_, base_, x], degree_] :=
{m, 0, 2 degree}]]]];
(* Preparing SDP for Export *)
-rho = SetPrecision[3-2 Sqrt[2], prec];
+rhoCrossing = SetPrecision[3-2 Sqrt[2], prec];
rescaledLaguerreSamplePoints[n_] := Table[
- SetPrecision[\[Pi]^2 (-1+4k)^2/(-64n Log[rho]), prec],
+ SetPrecision[\[Pi]^2 (-1+4k)^2/(-64n Log[rhoCrossing]), prec],
{k,0,n-1}];
maxIndexBy[l_,f_] := SortBy[
@@ -106,7 +128,6 @@ WriteBootstrapSDP[file_, SDP[objective_, normalization_, positiveMatricesWithPre
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[
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/sdpb.git
More information about the debian-science-commits
mailing list