Skip to content

Commit b85906d

Browse files
author
tomoki.ohtsuki
committed
init for bitbucket
1 parent fbcabe2 commit b85906d

4 files changed

Lines changed: 186 additions & 88 deletions

File tree

‎examples/3dIsing.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
from subprocess import Popen, PIPE # To execute sdpb
55
from sage.cboot.scalar import context_for_scalar
66

7-
sdpb = None
7+
sdpb = "sdpb"#None
88

99
if not sdpb:
1010
raise RuntimeError("The path for sdpb not specified!")
@@ -64,9 +64,10 @@ def binary_seach_ising(delta,lower=1.0,upper=3.0):
6464
print("trying for Delta = {0}".format(try_delta))
6565
prob=ising_singlet_bound_SDP(delta,{0:try_delta})
6666
prob.write("3dIsing.xml")
67-
sdpbres=Popen([sdpb,"-s","test.xml","--findPrimalFeasible",\
67+
sdpbres=Popen([sdpb,"-s","3dIsing.xml","--findPrimalFeasible",\
6868
"--findDualFeasible","--noFinalCheckpoint"],\
6969
stderr=PIPE,stdout=PIPE).communicate()
70+
print sdpbres
7071
sol=re.compile(r"found ([^ ]) feasible").search(sdpbres[0])
7172
if not sol:
7273
raise RuntimeError("unexpected sdpb output")

‎examples/3dOn.py‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
import sage.cboot as cb
2+
import numpy as np
3+
from subprocess import Popen, PIPE
4+
import re
5+
6+
context=cb.context_for_scalar(epsilon=0.5,Lambda=13)
7+
lmax=25
8+
nu_max=12
9+
cbs={}
10+
for spin in range(0,lmax):
11+
g=context.approx_cb(nu_max,spin)
12+
cbs.update({spin:g})
13+
14+
def make_F(delta,sector,spin,gap_dict,NSO):
15+
delta=context(delta)
16+
try:
17+
gap=context(gap_dict[(sector,spin)])
18+
except KeyError:
19+
if spin==0:
20+
gap=context.epsilon
21+
else:
22+
gap=2*context.epsilon+spin
23+
g_shift=cbs[spin].shift(gap)
24+
25+
g_num=g_shift.matrix
26+
g_pref=g_shift.prefactor
27+
F=context.make_F_minus_matrix(delta).dot(g_num)
28+
H=context.make_F_plus_matrix(delta).dot(g_num)
29+
30+
if sector=="S":
31+
num=np.concatenate((context.null_ftype,F,H))
32+
33+
elif sector=="T":
34+
num=np.concatenate((F,(1-2/context(NSO))*F,-(1+2/context(NSO))*H))
35+
36+
elif sector=="A":
37+
num=np.concatenate((-F,F,-H))
38+
39+
return context.prefactor_numerator(g_pref,num)
40+
41+
def make_SDP(delta,gap_dict,NSO=2):
42+
delta=context(delta)
43+
pvms=[]
44+
for sector in ("S","T","A"):
45+
if sector is not "A":
46+
spins=[spin for spin in cbs.keys() if not spin%2]
47+
else:
48+
spins=[spin for spin in cbs.keys() if spin%2]
49+
for spin in spins:
50+
pvms.append(make_F(delta,sector,spin,gap_dict,NSO))
51+
52+
norm_F=context.make_F_minus_matrix(delta).dot(context.gBlock(0,0,0,0))
53+
norm_H=context.make_F_plus_matrix(delta).dot(context.gBlock(0,0,0,0))
54+
norm=np.concatenate((context.null_ftype,norm_F,norm_H))
55+
56+
obj=norm*0
57+
return context.SDP(norm,obj,pvms)
58+
59+
sdpb="sdpb"
60+
sdpbparams=["--findPrimalFeasible","--findDualFeasible","--noFinalCheckpoint"]
61+
62+
def bs(delta,upper=3,lower=1,sector="S",sdp_method=make_SDP,NSO=2):
63+
upper=context(upper)
64+
lower=context(lower)
65+
while upper - lower > 0.001:
66+
D_try=(upper+lower)/2
67+
prob=sdp_method(delta,{(sector,0):D_try},NSO=NSO)
68+
prob.write("3d_Ising_binary.xml")
69+
sdpbargs=[sdpb,"-s","3d_Ising_binary.xml"]+sdpbparams
70+
out, err=Popen(sdpbargs,stdout=PIPE,stderr=PIPE).communicate()
71+
sol=re.compile(r'found ([^ ]+) feasible').search(out).groups()[0]
72+
if sol=="dual":
73+
print("(Delta_phi, Delta_{1})={0} is excluded."\
74+
.format((float(delta),float(D_try)),sector))
75+
upper=D_try
76+
elif sol=="primal":
77+
print("(Delta_phi, Delta_{1})={0} is permitted."\
78+
.format((float(delta),float(D_try)),sector))
79+
lower=D_try
80+
else:
81+
raise RuntimeError("Unexpected return from sdpb")
82+
return upper
83+
84+
85+
if __name__=="__main__":
86+
# The default example
87+
print bs(0.52)
88+
89+
# ======================================
90+
# if you want to derive the bound on Delta_T
91+
#print bs(0.52,sector="T")
92+
93+

‎examples/3dOn2.py‎

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
import sage.cboot as cb
2+
import numpy as np
3+
from subprocess import Popen, PIPE
4+
import re
5+
6+
context=cb.context_for_scalar(epsilon=0.5,Lambda=13)
7+
lmax=25
8+
nu_max=12
9+
cbs={}
10+
for spin in range(0,lmax):
11+
g=context.approx_cb(nu_max,spin)
12+
cbs.update({spin:g})
13+
14+
def make_F(delta,sector,spin,gap_dict,NSO,Delta=None):
15+
delta=context(delta)
16+
if Delta==None:
17+
try:
18+
gap=context(gap_dict[(sector,spin)])
19+
except KeyError:
20+
if spin==0:
21+
gap=context.epsilon
22+
else:
23+
gap=2*context.epsilon+spin
24+
g=cbs[spin].shift(gap)
25+
else:
26+
Delta=context(Delta)
27+
g=context.gBlock(spin,Delta,0,0)
28+
F=context.dot(context.make_F_minus_matrix(delta),g)
29+
H=context.dot(context.make_F_plus_matrix(delta),g)
30+
31+
if sector=="S":
32+
return [0,F,H]
33+
34+
elif sector=="T":
35+
return [F,(1-2/context(NSO))*F,-(1+2/context(NSO))*H]
36+
37+
elif sector=="A":
38+
return [-F,F,-H]
39+
40+
def make_SDP(delta,gap_dict,NSO=2):
41+
delta=context(delta)
42+
pvms=[]
43+
for sector in ("S","T","A"):
44+
if sector is not "A":
45+
spins=[spin for spin in cbs.keys() if not spin%2]
46+
else:
47+
spins=[spin for spin in cbs.keys() if spin%2]
48+
for spin in spins:
49+
pvms.append(make_F(delta,sector,spin,gap_dict,NSO))
50+
51+
norm=make_F(delta,"S",0,None,NSO,Delta=0)
52+
obj=0
53+
return context.sumrule_to_SDP(norm,obj,pvms)
54+
55+
56+
sdpb="sdpb"
57+
sdpbparams=["--findPrimalFeasible","--findDualFeasible","--noFinalCheckpoint","--maxThreads","12"]
58+
59+
def bs(delta,upper=3,lower=1,sector="S",sdp_method=make_SDP,NSO=2):
60+
upper=context(upper)
61+
lower=context(lower)
62+
while upper - lower > 0.001:
63+
D_try=(upper+lower)/2
64+
prob=sdp_method(delta,{(sector,0):D_try},NSO=NSO)
65+
prob.write("3d_On_binary.xml")
66+
sdpbargs=[sdpb,"-s","3d_On_binary.xml"]+sdpbparams
67+
out, err=Popen(sdpbargs,stdout=PIPE,stderr=PIPE).communicate()
68+
sol=re.compile(r'found ([^ ]+) feasible').search(out).groups()[0]
69+
if sol=="dual":
70+
print("(Delta_phi, Delta_{1})={0} is excluded."\
71+
.format((float(delta),float(D_try)),sector))
72+
upper=D_try
73+
elif sol=="primal":
74+
print("(Delta_phi, Delta_{1})={0} is permitted."\
75+
.format((float(delta),float(D_try)),sector))
76+
lower=D_try
77+
else:
78+
raise RuntimeError("Unexpected return from sdpb")
79+
return upper
80+
81+
82+
if __name__=="__main__":
83+
# The default example
84+
print bs(0.52)
85+
86+
# ======================================
87+
# if you want to derive the bound on Delta_T
88+
print bs(0.52,sector="T")
89+
90+

‎examples/XY.py‎

Lines changed: 0 additions & 86 deletions
This file was deleted.

0 commit comments

Comments
 (0)