forked from tdupu/root-unitary
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathnonsimples_from_simples.sage
254 lines (243 loc) · 11 KB
/
nonsimples_from_simples.sage
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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
#####################################################
#### Process for computing and uploading new data ###
#####################################################
## In the following, I assume that the root-unitary
## and lmfdb folders are contained in a common root.
# `cd /PATH/TO/lmfdb`
# `./warwick.sh &`
# `sage`
## From the sage prompt:
# run `import lmfdb`
# `cd ../root-unitary`
# run `load('nonsimples_from_simples.sage')`
# run `make_simples(g, q)` for new (g, q=p^m)
## this will create new weil-simple-g%-q%.txt
## Make sure that you still have the files
## weil-simple-g%-q%.txt for smaller g (same q) and
## weil-all-g%-q%.txt for divisors of q (same g)
# run `nonsimples_from_simples(g, q)`
## this will create new weil-all-g%-q%.txt
## Alternatively, you can run
## all_nonsimples_from_simples with a range of (g, q)
# create dictionary qs with
# `qs[g] = [p^r for r in range(1,m+1)]`
# create dictionary models with
# `models = _make_models_dict(qs)`
# run `_fill_primitive_models(models, qs)`
## this will fill in the primitive models field
## Add more g and q to qs to fill more models
# run `_fill_angle_ranks(g, q)`
## this will fill in the angle_ranks field
# run `load("../lmfdb/lmfdb/import_scripts/abvar/fq_isog.py")`
# run `do_import_all('weil-all-g%s-q%s.txt'%(g, q))`
load("table_functions.sage")
load("simple_maker.sage")
from itertools import combinations_with_replacement, imap, izip_longest
import shutil
import datetime
def multiplicity_dict(L):
L = list(L)
return {a: L.count(a) for a in set(L)}
class CombosWithReplacement(object):
def __init__(self, L, m):
self.L = L
self.m = m
def __len__(self):
return binomial(len(self.L)+self.m-1,self.m)
def __iter__(self):
return combinations_with_replacement(self.L, self.m)
def nonsimples_from_simples(g,q,D=None,rootdir=None):
if rootdir is None:
rootdir = os.path.abspath(os.curdir)
print "Generating from simples g=%s, q=%s"%(g, q)
simplefilename = os.path.join(rootdir, "weil-simple-g%s-q%s.txt"%(g, q))
filename = os.path.join(rootdir, "weil-all-g%s-q%s.txt"%(g, q))
shutil.copy(simplefilename, filename)
p,r = q.is_prime_power(get_data=True)
polyRingF.<t> = PolynomialRing(Qp(p))
if D is None:
D = {}
for a in range(1,g):
D.update(load_previous_polys(q, a, rootdir))
with open(filename, 'a') as target:
for split in Partitions(g):
if len(split) == 1:
continue
split_mD = multiplicity_dict(split)
#for a,c in sorted(split_mD.iteritems()):
# print a, c, len(list(combinations_with_replacement(D[a,q],c)))
it = cartesian_product_iterator([CombosWithReplacement(D[a,q],c) for a,c in sorted(split_mD.items())])
for factors in it:
if len(factors) == 1:
factors = factors[0]
else:
factors = sum(factors,())
Lpoly = prod(factor.poly for factor in factors)
if Lpoly.degree() != 2*g:
print factors
raise RuntimeError
Lpoly_mD = multiplicity_dict([factor.label for factor in factors])
seen = set()
labels_uniq = [factor.label for factor in factors if not (factor.label in seen or seen.add(factor.label))] # keeps order
decomposition = [[lbl, Lpoly_mD[lbl]] for lbl in labels_uniq]
label = make_label(g, q, Lpoly)
angle_numbers = sorted(sum([factor.angle_numbers for factor in factors], []))
p_rank = sum(factor.p_rank for factor in factors)
slopes = [str(a) for a in sorted(sum([[QQ(str(s)) for s in factor.slopes] for factor in factors], []))]
invs = sum([factor.invs for factor in factors], [])
places = sum([factor.places for factor in factors], [])
jac = know_nonsimple_jac(g, q, p, r, [factor.poly for factor in factors], p_rank)
polydata = NonsimplePolyData(g, q, p, r, label, angle_numbers, p_rank, slopes, jac, decomposition, invs, places, polyRingF)
target.write(create_line(Lpoly, polydata, False))
def all_nonsimples_from_simples(qstart=None, gstart=None, rootdir=None):
if rootdir is None:
rootdir = os.path.abspath(os.curdir)
gs = defaultdict(list)
for filename in os.listdir(rootdir):
match = simplematcher.match(filename)
if match:
gf, qf = map(Integer, match.groups())
gs[qf].append(gf)
for q in gs:
if qstart is not None and q < qstart:
continue
gs[q].sort()
if gs[q][-1] > 1:
D = {}
for a in range(1,gs[q][-1]):
D.update(load_previous_polys(q, a, rootdir))
for g in gs[q]:
if gstart is not None:
if g < gstart:
continue
gstart = None
if g == 1:
simplefilename = "weil-simple-g%s-q%s.txt"%(g, q)
filename = "weil-all-g%s-q%s.txt"%(g, q)
shutil.copy(simplefilename, filename)
else:
nonsimples_from_simples(g, q, D, rootdir)
def _fill_angle_rank(g, q, rootdir=None):
rootdir = rootdir or os.path.abspath(os.curdir)
R = PolynomialRing(ZZ,'x')
all_filename = os.path.join(rootdir, 'weil-all-g%s-q%s.txt'%(g, q))
tmp_filename = os.path.join(rootdir, 'weil-tmp-g%s-q%s.txt'%(g, q))
print "Starting Angle rank (g=%s, q=%s) computation"%(g, q)
with open(all_filename) as Fall:
for num_lines, line in enumerate(Fall,1r):
pass
# Now num_lines is the number of lines in Fall
with open(all_filename) as Fall:
# a+ creates file if it doesn't exist.
with open(tmp_filename, 'a+') as Fwrite:
with open(tmp_filename) as Ftmp:
# Angle ranks take a long time to compute,
# so we keep sum_of_times and sum_of_squares
# to give good progress estimates
print_next_time = walltime()
sum_of_times = float(0r)
sum_of_squares = float(0r)
start_line = None
for cur_line, (line_all, line_tmp) in enumerate(izip_longest(Fall, Ftmp, fillvalue=None),1r):
if line_tmp is not None:
if line_all[:line_all.find(',')] != line_tmp[:line_tmp.find(',')]:
raise RuntimeError("Label mismatch")
if cur_line % 1000 == 0:
print "Skipping existing computations (%s/%s)"%(cur_line, num_lines)
continue
if start_line is None:
start_line = cur_line
data = json.loads(line_all.strip())
if data[5] != "":
Fwrite.write(line_all.strip() + "\n")
continue
t = walltime()
Ppoly = R(list(reversed(data[3])))
data[5] = int(num_angle_rank(Ppoly))
t = walltime(t)
sum_of_times += t
sum_of_squares += t^2
if walltime() > print_next_time:
print_next_time = walltime() + 15
to_print = "Angle rank (g=%s, q=%s) %s/%s."%(g, q, cur_line, num_lines)
if cur_line - start_line > 10:
elapsed_lines = cur_line - start_line + 1
scaling = float(num_lines - elapsed_lines) / elapsed_lines
# sigma^2 = 1/n sum(xj^2) - 2mu * 1/n sum(xj) + 1/n sum(mu^2)
# = 1/n sum(xj^2) - 2mu * mu + mu^2
# = 1/n sum(xj^2) - 1/n^2 sum(xj)^2
# sigma = sqrt( (1/n sum(xj^2) - 1/n^2 sum(xj)^2) )
# actual mean is sum(xj)/n +- 2sigma/sqrt(n)
sigma = sqrt(sum_of_squares - sum_of_times^2 / elapsed_lines)
lower_bound = max(0r, sum_of_times*scaling - 2*sigma*sqrt(scaling))
upper_bound = sum_of_times*scaling + 2*sigma*sqrt(scaling)
lower_bound = datetime.timedelta(seconds=int(lower_bound))
upper_bound = datetime.timedelta(seconds=int(upper_bound))
to_print += " %s to %s remaining."%(lower_bound, upper_bound)
print to_print
Fwrite.write(json.dumps(data) + "\n")
shutil.move(tmp_filename, all_filename)
def fill_angle_ranks(rootdir=None):
rootdir = rootdir or os.path.abspath(os.curdir)
gs = defaultdict(list)
for filename in os.listdir(rootdir):
match = allmatcher.match(filename)
if match:
gf, qf = map(int, match.groups())
gs[qf].append(gf)
for q in sorted(gs.keys()):
gs[q].sort()
for g in gs[q]:
_fill_angle_rank(g, q, rootdir)
def _fill_primitive_models(models, qs, rootdir=None):
rootdir = rootdir or os.path.abspath(os.curdir)
R = PolynomialRing(QQ,'x')
for g in qs:
for q in qs[g]:
filename = os.path.join(rootdir, "weil-all-g%s-q%s.txt"%(g, q))
s = ""
print filename
with open(filename) as F:
for line in F.readlines():
data = json.loads(line.strip())
Lpoly = R(data[3])
if Lpoly in models:
data[15] = models[Lpoly]
else:
data[15] = [data[0]]
s += json.dumps(data) + "\n"
try:
with open(filename, 'w') as F:
F.write(s)
except KeyboardInterrupt:
with open(filename, 'w') as F:
F.write(s)
raise
def _make_models_dict(qs, rootdir=None):
rootdir = rootdir or os.path.abspath(os.curdir)
models = defaultdict(list)
D = {}
for g in qs:
qs[g].sort()
for q in qs[g]:
if q^2 in qs[g]:
D.update(load_previous_polys(q,g,rootdir,True))
r = 2
while q^r in qs[g]:
for factor in D[g, q]:
if factor.poly in models: # already a base change
continue
models[base_change(factor.poly, r)].append(factor.label)
r += 1
return models
def fill_primitive_models(rootdir=None):
# Assumes that, for fixed g, if q^r exists then all smaller powers do as well.
rootdir = rootdir or os.path.abspath(os.curdir)
qs = defaultdict(list)
for filename in os.listdir(rootdir):
match = allmatcher.match(filename)
if match:
gf, qf = map(int, match.groups())
qs[gf].append(qf)
models = _make_models_dict(qs, rootdir)
_fill_primitive_models(models, qs, rootdir)