-
Notifications
You must be signed in to change notification settings - Fork 2
/
topy.m
145 lines (124 loc) · 6.18 KB
/
topy.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
BeginPackage["ToPython`"]
pyeval::usage = "pyeval[x]"
createPython::usage = "createPython[secs,scalarnum,vals,rmats,smats,umats,filename,opinfo]"
Begin["`Private`"]
pyeval[n_Integer] := TemplateApply["context(\"``\")", n]
pyeval[a_Times] := StringRiffle[pyeval /@ List @@ a, {"(", " * ", ")"}]
pyeval[a_Plus] := StringRiffle[pyeval /@ List @@ a, {"(", " + ", ")"}]
pyeval[Power[a_, b_]] := TemplateApply["pow(``, ``)", {pyeval[a], pyeval[b]}]
pyeval[Surd[x_, n_]] /; x >= 0 := TemplateApply["pow(``, ``)", {pyeval[x], pyeval[1/n]}]
pyeval[Surd[x_, n_]] := pyeval[-Surd[-x, n]]
pyeval[Sqrt[a_]] := TemplateApply["sqrt(``)", pyeval[a]]
pyeval[Rational[a_, b_]] := TemplateApply["(`` / ``)", {pyeval[a], pyeval[b]}]
pyeval[Pi] = "context(pi)";
pyeval[E] = "context(e)";
pyeval[EulerGamma] = "context(euler_gamma)";
pyeval[Catalan] = "context(catalan)";
pyeval[Khinchin] = "context(khinchin)";
pyeval[Glaisher] = "context(glaisher)";
pyeval[Sin[a_]] := TemplateApply["sin(``)", pyeval[a]]
pyeval[Cos[a_]] := TemplateApply["cos(``)", pyeval[a]]
pyeval[Tan[a_]] := TemplateApply["tan(``)", pyeval[a]]
pyeval[Sec[a_]] := TemplateApply["sec(``)", pyeval[a]]
pyeval[Csc[a_]] := TemplateApply["csc(``)", pyeval[a]]
pyeval[Cot[a_]] := TemplateApply["cot(``)", pyeval[a]]
pyeval[Sinh[a_]] := TemplateApply["sinh(``)", pyeval[a]]
pyeval[Cosh[a_]] := TemplateApply["cosh(``)", pyeval[a]]
pyeval[Tanh[a_]] := TemplateApply["tanh(``)", pyeval[a]]
pyeval[Sech[a_]] := TemplateApply["sech(``)", pyeval[a]]
pyeval[Csch[a_]] := TemplateApply["csch(``)", pyeval[a]]
pyeval[Coth[a_]] := TemplateApply["coth(``)", pyeval[a]]
pyeval[ArcSin[a_]] := TemplateApply["asin(``)", pyeval[a]]
pyeval[ArcCos[a_]] := TemplateApply["acos(``)", pyeval[a]]
pyeval[ArcTan[a_]] := TemplateApply["atan(``)", pyeval[a]]
pyeval[ArcSec[a_]] := TemplateApply["asec(``)", pyeval[a]]
pyeval[ArcCsc[a_]] := TemplateApply["acsc(``)", pyeval[a]]
pyeval[ArcCot[a_]] := TemplateApply["acot(``)", pyeval[a]]
pyeval[ArcSinh[a_]] := TemplateApply["asinh(``)", pyeval[a]]
pyeval[ArcCosh[a_]] := TemplateApply["acosh(``)", pyeval[a]]
pyeval[ArcTanh[a_]] := TemplateApply["atanh(``)", pyeval[a]]
pyeval[ArcSech[a_]] := TemplateApply["asech(``)", pyeval[a]]
pyeval[ArcCsch[a_]] := TemplateApply["acsch(``)", pyeval[a]]
pyeval[ArcCoth[a_]] := TemplateApply["acoth(``)", pyeval[a]]
pyeval[Exp[a_]] := TemplateApply["exp(``)", pyeval[a]]
pyeval[Log[a_]] := TemplateApply["log(``)", pyeval[a]]
pyeval[Log[b_, a_]] := TemplateApply["log(``, ``)", {pyeval[a], pyeval[b]}]
pyeval[I] = "I";
pyeval[x_?ExactNumberQ] := TemplateApply["context(\"``\")", ToString @ FortranForm @ N[x, 300]]
pyeval[x_?NumberQ] := Module[{r = Rationalize[x, 0]}, If[Abs[Denominator[r]] + Abs[Numerator[r]] < 10000, pyeval[r], TemplateApply["context(\"``\")", ToString @ FortranForm[x]]]]
createPython[secs_, scalarnum_, vals_, rmats_, smats_, umats_, filename_, opinfo_] :=
TemplateApply[template, <|"secs"->secs, "scalarnum"->scalarnum, "vals"->vals, "rmats"->rmats, "smats"->smats, "umats"->umats, "filename"->filename, "opinfo"->opinfo|>]
template = "# -*- coding: utf-8 -*-
from __future__ import print_function
import sage.cboot as cb
from sage.misc.cachefunc import cached_function
import os.path
import sys
from sage.rings.rational import Rational
from sage.all import pi, e, euler_gamma, catalan, khinchin, glaisher, sin, cos, tan, sec, csc, cot, sinh, cosh, tanh, sech, csch, coth, asin, acos, atan, asec, acsc, acot, asinh, acosh, atanh, asech, acsch, acoth, sqrt, log, exp, I
# <-- edit from here
context = cb.context_for_scalar(epsilon=0.5, Lambda=15, Prec=800)
spins = list(range(27)) + [49, 50]
nu_max = 14
`opinfo`
# {(r, s): d, ...} means operators in irrep r and spin s must have Delta >= d.
mygap = {}
def gaps(deltas):
return mygap
save_dir = \".\"
def name(deltas):
return \"sdp-`filename`\".format(deltas).replace(\"/\", \"#\")
# edit to here -->
# do not edit from here
secs = {`secs`}
scalarnum = `scalarnum`
@cached_function
def prepare_g(spin, Delta_12, Delta_34):
return context.approx_cb(nu_max, spin, Delta_1_2=Delta_12, Delta_3_4=Delta_34)
def get_shift(gaps, sector, spin):
if (sector, spin) in gaps: return context(gaps[(sector, spin)])
elif spin == 0: return context.epsilon
else: return 2 * context.epsilon + spin
val = `vals`
one = context(1)
def make_F_range(delta, sector, num, spin):
shift = get_shift(gaps(delta), sector, spin)
sign = one if spin % 2 == 0 else -one
F = lambda d1, d2, d3, d4: sign * context.dot(context.F_minus_matrix((d2 + d3) / 2), prepare_g(spin, d1 - d2, d3 - d4).shift(shift))
H = lambda d1, d2, d3, d4: sign * context.dot(context.F_plus_matrix((d2 + d3) / 2), prepare_g(spin, d1 - d2, d3 - d4).shift(shift))
get = lambda fh, d1, d2, d3, d4: fh(delta[d1], delta[d2], delta[d3], delta[d4])
`rmats`
raise RuntimeError(\"unknown sector name\")
def make_F_scalar(delta, num):
F = lambda d1, d2, d3, d4, d: context.dot(context.F_minus_matrix((d2 + d3) / 2), context.gBlock(0, d, d1 - d2, d3 - d4))
H = lambda d1, d2, d3, d4, d: context.dot(context.F_plus_matrix((d2 + d3) / 2), context.gBlock(0, d, d1 - d2, d3 - d4))
get = lambda fh, d1, d2, d3, d4, d: fh(delta[d1], delta[d2], delta[d3], delta[d4], delta[d])
`smats`
raise RuntimeError(\"unknown sector number\")
def make_F_unit(delta):
F = lambda d1, d2, d3, d4: context.dot(context.F_minus_matrix((d2 + d3) / 2), context.gBlock(0, 0, d1 - d2, d3 - d4))
H = lambda d1, d2, d3, d4: context.dot(context.F_plus_matrix((d2 + d3) / 2), context.gBlock(0, 0, d1 - d2, d3 - d4))
get = lambda fh, d1, d2, d3, d4: fh(delta[d1], delta[d2], delta[d3], delta[d4])
`umats`
def make_SDP(delta):
cdel = dict()
for k in delta: cdel[k] = context(delta[k])
pvms = []
for sec in secs:
for spin in spins:
if spin % 2 == sec[1]:
for num in range(secs[sec]): pvms.append(make_F_range(cdel, sec[0], num, spin))
for num in range(scalarnum): pvms.append(make_F_scalar(cdel, num))
norm = make_F_unit(cdel)
obj = 0
return context.sumrule_to_SDP(norm, obj, pvms)
def write_SDP(deltas, dir=None, file=None, overwrite=False):
if dir is None: dir = save_dir
if file is None: file = name(deltas) + \".xml\"
path = os.path.join(dir, file)
if overwrite or not os.path.exists(path): make_SDP(deltas).write(path)
return path
"
Protect[pyeval, createPython]
End[ ]
EndPackage[ ]