Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/interfaces/phc.py
4045 views
1
r"""
2
Interface to PHC.
3
4
PHC computes numerical information about systems of polynomials over
5
the complex numbers.
6
7
PHC implements polynomial homotopy methods to exploit structure in
8
order to better approximate all isolated solutions. The package also
9
includes extra tools to handle positive dimensional solution
10
components.
11
12
AUTHOR:
13
-- PHC was written by J. Verschelde, R. Cools, and many others (?)
14
-- William Stein and Kelly ?? -- first version of interface to PHC
15
-- Marshall Hampton -- second version of interface to PHC
16
-- Marshall Hampton and Alex Jokela -- third version, path tracking
17
18
"""
19
20
########################################################################
21
# Copyright (C) 2006 William Stein <[email protected]>
22
# Copyright (C) 2008 Marshall Hampton <[email protected]>
23
#
24
# Distributed under the terms of the GNU General Public License (GPL)
25
#
26
# The full text of the GPL is available at:
27
#
28
# http://www.gnu.org/licenses/
29
########################################################################
30
31
import os
32
33
import sage.misc.misc
34
from sage.rings.real_mpfr import RR
35
from sage.rings.all import CC
36
from sage.rings.integer import Integer
37
from sage.plot.plot import *
38
import re
39
40
import pexpect
41
42
def get_solution_dicts(output_file_contents, input_ring, get_failures = True):
43
"""
44
Returns a list of dictionaries of variable:value (key:value)
45
pairs. Only used internally; see the solution_dict function in
46
the PHC_Object class definition for details.
47
48
INPUT:
49
output_file_contents -- phc solution output as a string
50
input_ring -- a PolynomialRing that variable names can be coerced into
51
52
OUTPUT:
53
a list of dictionaries of solutions
54
55
EXAMPLES:
56
sage: from sage.interfaces.phc import * #optional
57
sage: R2.<x1,x2> = PolynomialRing(QQ,2) #optional
58
sage: test_sys = [(x1-1)^5-x2, (x2-1)^5-1] #optional
59
sage: sol = phc.blackbox(test_sys, R2) #optional
60
sage: test = get_solution_dicts(sol.output_file_contents,R2) #optional
61
sage: str(sum([q[x1].real() for q in test]))[0:4] #optional
62
'25.0'
63
"""
64
output_list = output_file_contents.splitlines()
65
test = 'False'
66
solution_dicts = []
67
for solution_line in range(len(output_list)-1,-1,-1):
68
if output_list[solution_line].find('THE SOLUTIONS') == 0:
69
break
70
try:
71
var_number = int(output_list[solution_line+2].split(' ')[1])
72
sol_number = int(output_list[solution_line+2].split(' ')[0])
73
except IndexError:
74
var_number = int(output_list[solution_line+1].split(' ')[1])
75
sol_number = int(output_list[solution_line+1].split(' ')[0])
76
for i in range(solution_line + 1,len(output_list)):
77
if output_list[i].count('the solution for t') == 1:
78
if output_list[i-3].count('success') > 0 or get_failures == True:
79
temp_dict = {}
80
for j in range(1,var_number+1):
81
rawsplit = output_list[i+j].split(': ')[1].split(' ')
82
for extras in range(rawsplit.count('')):
83
rawsplit.remove('')
84
temp_var = output_list[i+j].split(': ')[0].replace(' ','')
85
temp_dict[input_ring(temp_var)] = CC(rawsplit[0],rawsplit[1])
86
solution_dicts.append(temp_dict)
87
return solution_dicts
88
89
def get_classified_solution_dicts(output_file_contents, input_ring, get_failures = True):
90
"""
91
Returns a dictionary of lists of dictionaries of variable:value (key:value)
92
pairs. Only used internally; see the classified_solution_dict function in
93
the PHC_Object class definition for details.
94
95
INPUT:
96
output_file_contents -- phc solution output as a string
97
input_ring -- a PolynomialRing that variable names can be coerced into
98
99
OUTPUT:
100
a dictionary of lists if dictionaries of solutions, classifies by type
101
102
EXAMPLES:
103
sage: from sage.interfaces.phc import * #optional
104
sage: R2.<x1,x2> = PolynomialRing(QQ,2) #optional
105
sage: test_sys = [(x1-2)^5-x2, (x2-1)^5-1] #optional
106
sage: sol = phc.blackbox(test_sys, R2) #optional
107
sage: sol_classes = get_classified_solution_dicts(sol.output_file_contents,R2) #optional
108
sage: len(sol_classes['real']) #optional
109
1
110
"""
111
output_list = output_file_contents.splitlines()
112
test = 'False'
113
solution_dicts = {}
114
solution_types = ['complex', 'real','failure']
115
for sol_type in solution_types:
116
solution_dicts[sol_type] = []
117
for solution_line in range(len(output_list)-1,-1,-1):
118
if output_list[solution_line].find('THE SOLUTIONS') == 0:
119
break
120
var_number = int(output_list[solution_line+2].split(' ')[1])
121
sol_number = int(output_list[solution_line+2].split(' ')[0])
122
for i in range(solution_line + 1,len(output_list)):
123
if output_list[i].count('the solution for t') == 1:
124
phc_type = output_list[i+var_number+1].split(' = ')[-1]
125
if phc_type.find('complex') != -1:
126
phc_type = 'complex'
127
elif phc_type.find('real') != -1:
128
phc_type = 'real'
129
else:
130
phc_type = 'failure'
131
temp_dict = {}
132
for j in range(1,var_number+1):
133
rawsplit = output_list[i+j].split(': ')[1].split(' ')
134
for extras in range(rawsplit.count('')):
135
rawsplit.remove('')
136
temp_var = output_list[i+j].split(': ')[0].replace(' ','')
137
if phc_type == 'real':
138
temp_dict[input_ring(temp_var)] = RR(rawsplit[0])
139
else:
140
temp_dict[input_ring(temp_var)] = CC(rawsplit[0],rawsplit[1])
141
solution_dicts[phc_type].append(temp_dict)
142
return solution_dicts
143
144
def get_variable_list(output_file_contents):
145
"""
146
Returns the variables, as strings, in the order in which PHCpack has processed them.
147
148
EXAMPLES:
149
sage: from sage.interfaces.phc import * #optional
150
sage: R2.<x1,x2> = PolynomialRing(QQ,2) #optional
151
sage: test_sys = [(x1-2)^5-x2, (x2-1)^5-1] #optional
152
sage: sol = phc.blackbox(test_sys, R2) #optional
153
sage: get_variable_list(sol.output_file_contents) #optional
154
['x1', 'x2']
155
"""
156
output_list = output_file_contents.splitlines()
157
for solution_line in range(len(output_list)-1,-1,-1):
158
if output_list[solution_line].find('THE SOLUTIONS') == 0:
159
break
160
var_number = int(output_list[solution_line+2].split(' ')[1])
161
varlist = []
162
for var_ind in range(var_number):
163
var = output_list[solution_line + 8 + var_ind].split(' ')[1]
164
varlist.append(var)
165
return varlist
166
167
class PHC_Object:
168
169
def __init__(self, output_file_contents, input_ring):
170
"""
171
A container for data from the PHCpack program - lists of float
172
solutions, etc. Currently the file contents are kept as a string;
173
for really large outputs this would be bad.
174
175
INPUT:
176
output_file_contents: the string output of PHCpack
177
input_ring: for coercion of the variables into the desired ring.
178
179
EXAMPLES:
180
sage: from sage.interfaces.phc import phc #optional
181
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
182
sage: start_sys = [(x-1)^2+(y-1)-1, x^2+y^2-1] #optional
183
sage: sol = phc.blackbox(start_sys, R2) #optional
184
sage: str(sum([x[0] for x in sol.solutions()]).real())[0:3] #optional
185
'2.0'
186
"""
187
self.output_file_contents = output_file_contents
188
self.input_ring = input_ring
189
190
191
def save_as_start(self, start_filename = None, sol_filter = ''):
192
"""
193
Saves a solution as a phcpack start file. The usual output is
194
just as a string, but it can be saved to a file as well. Even
195
if saved to a file, it still returns the output string.
196
197
EXAMPLES:
198
sage: from sage.interfaces.phc import phc #optional
199
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
200
sage: start_sys = [x^3-y^2,y^5-1] #optional
201
sage: sol = phc.blackbox(start_sys, R2) #optional
202
sage: start_save = sol.save_as_start() #optional
203
sage: end_sys = [x^7-2,y^5-x^2] #optional
204
sage: sol = phc.start_from(start_save, end_sys, R2) #optional
205
sage: len(sol.solutions()) #optional
206
15
207
"""
208
start_data = ''
209
output_list = self.output_file_contents.splitlines()
210
for a_line in output_list:
211
if a_line.find('ROOT COUNTS') != -1 or a_line.find('START SOLUTIONS') != -1:
212
break
213
else:
214
start_data += a_line + '\n'
215
for index in range(len(output_list)-1,0,-1):
216
a_line = output_list[index]
217
if a_line.find('THE SOLUTIONS') != -1:
218
found_solutions = index
219
break
220
start_data += output_list[found_solutions] + '\n\n'
221
try:
222
var_number = int(output_list[found_solutions+1].split(' ')[1])
223
except:
224
# bad error handling
225
var_number = int(output_list[found_solutions+2].split(' ')[1])
226
sol_count = 0
227
sol_data = ''
228
for i in range(found_solutions + 2, len(output_list)):
229
if output_list[i].count('the solution for t') == 1 and output_list[i+1+var_number].find(sol_filter) != -1:
230
phc_type = output_list[i+var_number+1].split(' = ')[-1]
231
if phc_type.find('no solution') == -1:
232
sol_count += 1
233
for ind2 in range(i-3,i+var_number+2):
234
sol_data += output_list[ind2] + '\n'
235
jan_bar = '===========================================================================\n'
236
sol_data += jan_bar
237
start_data += str(sol_count) + ' ' + str(var_number) + '\n'
238
start_data += jan_bar + sol_data
239
if start_filename != None:
240
start_file = file(start_filename,'w')
241
start_file.write(start_data)
242
start_file.close()
243
return start_data
244
245
def classified_solution_dicts(self):
246
"""
247
Returns a dictionary of lists of dictionaries of solutions.
248
Its not as crazy as it sounds; the keys are the types of solutions as
249
classified by phcpack: regular vs. singular, complex vs. real
250
251
INPUT:
252
None
253
254
OUTPUT:
255
A dictionary of lists of dictionaries of solutions
256
257
EXAMPLES:
258
sage: from sage.interfaces.phc import phc #optional
259
sage: R.<x,y> = PolynomialRing(CC,2) #optional
260
sage: p_sys = [x^10-y,y^2-1] #optional
261
sage: sol = phc.blackbox(p_sys,R) #optional
262
sage: classifieds = sol.classified_solution_dicts() #optional
263
sage: str(sum([q[y] for q in classifieds['real']]))[0:3] #optional
264
'2.0'
265
"""
266
try:
267
return self.__classified_sols
268
except AttributeError:
269
pass
270
classified_sols = get_classified_solution_dicts(self.output_file_contents, self.input_ring)
271
self.__classified_sols = classified_sols
272
return classified_sols
273
274
def solution_dicts(self, get_failures = False):
275
"""
276
Returns a list of solutions in dictionary form: variable:value.
277
278
INPUT:
279
self -- for access to self_out_file_contents, the string
280
of raw PHCpack output.
281
282
get_failures (optional) -- a boolean. The default (False)
283
is to not process failed homotopies. These either lie on
284
positive-dimensional components or at infinity.
285
286
OUTPUT:
287
solution_dicts: a list of dictionaries. Each dictionary
288
element is of the form variable:value, where the variable
289
is an element of the input_ring, and the value is in
290
ComplexField.
291
292
EXAMPLES:
293
sage: from sage.interfaces.phc import * #optional
294
sage: R.<x,y,z> = PolynomialRing(QQ,3) #optional
295
sage: fs = [x^2-1,y^2-x,z^2-y] #optional
296
sage: sol = phc.blackbox(fs,R) #optional
297
sage: s_list = sol.solution_dicts() #optional
298
sage: s_list.sort() #optional
299
sage: s_list[0] #optional
300
{y: 1.00000000000000, z: -1.00000000000000, x: 1.00000000000000}
301
"""
302
try:
303
return self.__solution_dicts
304
except AttributeError:
305
pass
306
solution_dicts = get_solution_dicts(self.output_file_contents, self.input_ring, get_failures = get_failures)
307
self.__solution_dicts = solution_dicts
308
return solution_dicts
309
310
def solutions(self, get_failures = False):
311
"""
312
Returns a list of solutions in the ComplexField. Use the variable_list function to get the order of variables used by PHCpack, which is usually different than the term order of the input_ring.
313
314
INPUT:
315
self -- for access to self_out_file_contents, the string
316
of raw PHCpack output.
317
get_failures (optional) -- a boolean. The default (False)
318
is to not process failed homotopies. These either lie on
319
positive-dimensional components or at infinity.
320
321
OUTPUT:
322
solutions: a list of lists of ComplexField-valued solutions.
323
324
EXAMPLES:
325
sage: from sage.interfaces.phc import * #optional
326
sage: R2.<x1,x2> = PolynomialRing(QQ,2) #optional
327
sage: test_sys = [x1^5-x1*x2^2-1, x2^5-x1*x2-1] #optional
328
sage: sol = phc.blackbox(test_sys, R2) #optional
329
sage: len(sol.solutions()) #optional
330
25
331
"""
332
try:
333
return self.__solutions
334
except AttributeError:
335
pass
336
solution_dicts = get_solution_dicts(self.output_file_contents, self.input_ring, get_failures = get_failures)
337
self.__solution_dicts = solution_dicts
338
solutions = [sol_dict.values() for sol_dict in solution_dicts]
339
self.__solutions = solutions
340
return solutions
341
342
def variable_list(self):
343
"""
344
Returns the variables, as strings, in the order in which
345
PHCpack has processed them.
346
347
EXAMPLES:
348
sage: from sage.interfaces.phc import * #optional
349
sage: R2.<x1,x2> = PolynomialRing(QQ,2) #optional
350
sage: test_sys = [x1^5-x1*x2^2-1, x2^5-x1*x2-1] #optional
351
sage: sol = phc.blackbox(test_sys, R2) #optional
352
sage: sol.variable_list() #optional
353
['x1', 'x2']
354
"""
355
try:
356
return self.__var_list
357
except AttributeError:
358
pass
359
var_list = get_variable_list(self.output_file_contents)
360
self.__var_list = var_list
361
return var_list
362
363
class PHC:
364
"""
365
A class to interface with PHCpack, for computing numerical
366
homotopies and root counts.
367
368
EXAMPLES:
369
sage: from sage.interfaces.phc import phc #optional
370
sage: R.<x,y> = PolynomialRing(CDF,2) #optional
371
sage: testsys = [x^2 + 1, x*y - 1] #optional
372
sage: phc.mixed_volume(testsys) # optional -- you must have phc install
373
2
374
sage: v = phc.blackbox(testsys, R) # optional
375
sage: sols = v.solutions() # optional
376
sage: sols.sort() # optional
377
sage: sols # optional
378
[[-1.00000000000000*I, 1.00000000000000*I], [1.00000000000000*I, -1.00000000000000*I]]
379
sage: sol_dict = v.solution_dicts() # optional
380
sage: x_sols_from_dict = [d[x] for d in sol_dict] # optional
381
sage: x_sols_from_dict.sort(); x_sols_from_dict # optional
382
[-1.00000000000000*I, 1.00000000000000*I]
383
sage: residuals = [[test_equation.change_ring(CDF).subs(sol) for test_equation in testsys] for sol in v.solution_dicts()] # optional
384
sage: residuals # optional
385
[[0, 0], [0, 0]]
386
"""
387
388
def _output_from_command_list(self, command_list, polys, verbose = False):
389
"""
390
A pexpect interface to phcpack, given a command list for
391
interactive dialogs. The input file is supplied from the
392
polynomial list, output file is also supplied. This is
393
only used as a building block for the interface.
394
395
INPUT:
396
command_list -- a list of commands to phc
397
polys -- a polynomial system as a list of polynomials
398
399
OUTPUT:
400
an output string from phc
401
402
EXAMPLES:
403
sage: from sage.interfaces.phc import * #optional
404
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
405
sage: start_sys = [(x-1)^2+(y-1)-1, x^2+y^2-1] #optional
406
sage: a = phc._output_from_command_list(['phc -m','4','n','n','n'], start_sys)#optional
407
"""
408
# Get temporary file names (these will be in SAGE_HOME/.sage/tmp/pid)
409
input_filename = sage.misc.misc.tmp_filename()
410
output_filename = sage.misc.misc.tmp_filename()
411
412
# Get the input polynomial text
413
input = self._input_file(polys)
414
if verbose:
415
print "Writing the input file to %s"%input_filename
416
open(input_filename,'w').write(input)
417
418
if verbose:
419
print "The following file will be the input polynomial file to phc."
420
print input
421
422
# Create a phc process
423
child_phc = pexpect.spawn(command_list[0])
424
# feed it the commands
425
child_phc.sendline('y')
426
child_phc.sendline(input_filename)
427
child_phc.sendline(output_filename)
428
for command_string in command_list[1:]:
429
if verbose:
430
print command_string
431
child_phc.sendline(command_string)
432
child_phc.expect('results')
433
read_stuff = child_phc.read()
434
if verbose:
435
print read_stuff
436
child_phc.close()
437
if not os.path.exists(output_filename):
438
raise RuntimeError, "The output file does not exist; something went wrong running phc."
439
440
# Delete the input file
441
os.unlink(input_filename)
442
443
# Return the output filename
444
return output_filename
445
446
def _input_file(self, polys):
447
"""
448
This is used internally to implement the PHC interface.
449
450
INPUT:
451
polys -- a list of polynomials in a Sage polynomial ring
452
over a field that embeds into the complex
453
numbers.
454
455
OUTPUT:
456
-- a PHC input file (as a text string) that
457
describes these polynomials.
458
459
EXAMPLES:
460
sage: from sage.interfaces.phc import * #optional
461
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
462
sage: start_sys = [(x-1)^2+(y-1)-1, x^2+y^2-1] #optional
463
sage: phc._input_file(start_sys) #optional
464
'2\nx^2 - 2*x + y - 1;\nx^2 + y^2 - 1;\n'
465
"""
466
if not isinstance(polys, (list, tuple)):
467
raise TypeError, 'polys must be a list or tuple'
468
s = '%s\n'%len(polys)
469
for f in polys:
470
s += f._repr_() + ';\n' # note the semicolon *terminators*
471
472
return s
473
474
def _parse_path_file(self, input_filename, verbose = False):
475
"""
476
Takes a phpack output file containing path tracking information
477
and parses it into a list of lists of dictionaries - i.e. a
478
list of solutions paths, where each solution path is a list of
479
dictionaries of variable and homotopy parameter values.
480
481
INPUT:
482
input_filename -- file must have path-tracking information
483
484
OUTPUT:
485
a list of lists of dictionaries, described above
486
487
EXAMPLES:
488
sage: from sage.interfaces.phc import * #optional
489
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
490
sage: start_sys = [x^5-y^2,y^5-1] #optional
491
sage: sol = phc.blackbox(start_sys, R2) #optional
492
sage: start_save = sol.save_as_start() #optional
493
sage: end_sys = [x^5-2,y^5-x^2] #optional
494
sage: path_track_filename = phc._path_track_file(start_save, end_sys, R2, c_skew = .001) #optional
495
sage: sol_paths = phc._parse_path_file(path_track_filename) #optional
496
sage: len(sol_paths) #optional
497
25
498
"""
499
500
if not os.path.exists(input_filename):
501
raise RuntimeError, "The file containing output from phc (" + input_filename + ") cannot be found"
502
503
fh = open(input_filename)
504
line_idx = 0
505
begin = 0
506
count = 0
507
solutions_dicts = []
508
steps_dicts = []
509
510
# regular expressions for matching certain output types
511
var_cnt_regex = re.compile('^ +([0-9]+)')
512
output_regex = re.compile('^OUTPUT INFORMATION DURING')
513
t_regex = re.compile('(^t +: +(-{0,1}[0-9]+\.[0-9]+E[-+][0-9]+) +(-{0,1}[0-9]+\.[0-9]+E[-+][0-9]+)$)', re.IGNORECASE)
514
sols_regex = re.compile('(^ *(([a-z]|[0-9])+) +: +(-?[0-9]+\.[0-9]+E[-+][0-9]+) +(-?[0-9]+\.[0-9]+E[-+][0-9]+)$)', re.IGNORECASE)
515
complete_regex= re.compile('^TIMING INFORMATION')
516
517
breakfast = False
518
a_line = fh.readline()
519
end_test = ''
520
while a_line:
521
# processing....
522
a_line = a_line.replace("\n", '')
523
if line_idx == 0:
524
m = var_cnt_regex.match(a_line)
525
if m:
526
count = Integer(m.group(1))
527
if count > 0:
528
m = output_regex.match(a_line)
529
if m:
530
begin = 1
531
if begin:
532
m = t_regex.match(a_line)
533
if m:
534
# put the t-values into a dict
535
# m.group(2) contains the real val
536
# m.group(3) contains the imaginary val
537
# fh_w.write( "T=> G1(" + m.group(2) + '),G2(' + m.group(3) + ")\n")
538
# read off two lines - this should be 'm' and 'the solution for t :'
539
a_line = fh.readline()
540
end_test = a_line # store this to check for end of solution
541
a_line = fh.readline()
542
t_val = CC(m.group(2), m.group(3))
543
temp_dict = {}
544
temp_dict["t"] = t_val
545
for i in range(0, count):
546
a_line = fh.readline()
547
m = sols_regex.match(a_line)
548
if m:
549
# m.group(2) contains our var name
550
# m.group(4) contains our real val
551
# m.group(5) contains our imaginary val
552
temp_dict[m.group(2)] = CC(m.group(4),m.group(5))
553
steps_dicts.append(temp_dict)
554
# check if its the end of a solution
555
if end_test.find('Length of path') != -1:
556
if verbose: print "recording sol"
557
if steps_dicts != []:
558
solutions_dicts.append(steps_dicts)
559
steps_dicts = []
560
m = complete_regex.match(a_line)
561
if m:
562
breakfast = True
563
if breakfast:
564
break
565
line_idx += 1
566
a_line = fh.readline()
567
fh.close()
568
return solutions_dicts
569
570
def _path_track_file(self, start_filename_or_string, polys, input_ring, c_skew = 0.001, verbose = False):
571
"""
572
Returns the filename which contains path tracking output.
573
574
EXAMPLES:
575
sage: from sage.interfaces.phc import * #optional
576
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
577
sage: start_sys = [x^6-y^2,y^5-1] #optional
578
sage: sol = phc.blackbox(start_sys, R2) #optional
579
sage: start_save = sol.save_as_start() #optional
580
sage: end_sys = [x^7-2,y^5-x^2] #optional
581
sage: path_track_filename = phc._path_track_file(start_save, end_sys, R2, c_skew = .001) #optional
582
sage: sol_paths = phc._parse_path_file(path_track_filename) #optional
583
sage: len(sol_paths) #optional
584
30
585
"""
586
# Probably unnecessarily redundant from the start_from function
587
if start_filename_or_string.find('THE SOLUTIONS') != -1:
588
start_filename = sage.misc.misc.tmp_filename()
589
start_file = file(start_filename,'w')
590
start_file.write(start_filename_or_string)
591
start_file.close()
592
elif os.path.exists(start_filename_or_string):
593
start_filename = start_filename_or_string
594
else:
595
raise RuntimeError, "There is something wrong with your start string or filename"
596
597
return self._output_from_command_list(['phc','0','0','A',start_filename, 'y','1','0','n','k','2','a','1',str(c_skew),'0','0','2'], polys, verbose = verbose)
598
599
def path_track(self, start_sys, end_sys, input_ring, c_skew = .001, saved_start = None):
600
"""
601
This function computes homotopy paths between the solutions of start_sys and end_sys.
602
603
INPUT:
604
start_sys -- a square polynomial system, given as a list of polynomials
605
end_sys -- same type as start_sys
606
input_ring -- for coercion of the variables into the desired ring.
607
c_skew -- optional. the imaginary part of homotopy multiplier; nonzero values
608
are often necessary to avoid intermediate path collisions
609
saved_start -- optional. A phc output file. If not given, start system solutions
610
are computed via the phc.blackbox function.
611
612
OUTPUT:
613
a list of paths as dictionaries, with the keys variables and t-values on the path.
614
615
EXAMPLES:
616
sage: from sage.interfaces.phc import * #optional
617
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
618
sage: start_sys = [x^6-y^2,y^5-1] #optional
619
sage: sol = phc.blackbox(start_sys, R2) #optional
620
sage: start_save = sol.save_as_start() #optional
621
sage: end_sys = [x^7-2,y^5-x^2] #optional
622
sage: sol_paths = phc.path_track(start_sys, end_sys, R2, saved_start = start_save) #optional
623
sage: len(sol_paths) #optional
624
30
625
"""
626
627
if not saved_start:
628
sol = phc.blackbox(start_sys, input_ring)
629
saved_start = sol.save_as_start()
630
path_track_filename = phc._path_track_file(saved_start, end_sys, input_ring = input_ring, c_skew = c_skew)
631
sol_paths = phc._parse_path_file(path_track_filename)
632
os.unlink(path_track_filename)
633
return sol_paths
634
635
def plot_paths_2d(self, start_sys, end_sys, input_ring, c_skew = .001, endpoints = True, saved_start = None, rand_colors = False):
636
"""
637
This returns a graphics object of solution paths in the complex plane.
638
639
INPUT:
640
start_sys -- a square polynomial system, given as a list of polynomials
641
end_sys -- same type as start_sys
642
input_ring -- for coercion of the variables into the desired ring.
643
c_skew -- optional. the imaginary part of homotopy multiplier; nonzero values
644
are often necessary to avoid intermediate path collisions
645
endpoints -- optional. Whether to draw in the ends of paths as points.
646
saved_start -- optional. A phc output file. If not given, start system solutions
647
are computed via the phc.blackbox function.
648
649
OUTPUT:
650
lines and points of solution paths
651
652
EXAMPLES:
653
sage: from sage.interfaces.phc import * #optional
654
sage: from sage.structure.sage_object import SageObject
655
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
656
sage: start_sys = [x^5-y^2,y^5-1] #optional
657
sage: sol = phc.blackbox(start_sys, R2) #optional
658
sage: start_save = sol.save_as_start() #optional
659
sage: end_sys = [x^5-25,y^5-x^2] #optional
660
sage: testing = phc.plot_paths_2d(start_sys, end_sys, R2) #optional
661
sage: type(testing) #optional (normally use plot here)
662
<class 'sage.plot.graphics.Graphics'>
663
"""
664
paths = phc.path_track(start_sys, end_sys, input_ring, c_skew = c_skew, saved_start = saved_start)
665
path_lines = []
666
sol_pts = []
667
if rand_colors:
668
r_color = {}
669
for a_var in input_ring.gens():
670
var_name = str(a_var)
671
r_color[var_name] = (random(),random(),random())
672
for a_sol in paths:
673
for a_var in input_ring.gens():
674
var_name = str(a_var)
675
temp_line = []
676
for data in a_sol:
677
temp_line.append([data[var_name].real(), data[var_name].imag()])
678
if rand_colors:
679
path_lines.append(line(temp_line, rgbcolor = r_color[var_name]))
680
else:
681
path_lines.append(line(temp_line))
682
if endpoints:
683
sol_pts = []
684
for a_sol in paths:
685
for a_var in input_ring.gens():
686
var_name = str(a_var)
687
sol_pts.append(point([a_sol[0][var_name].real(), a_sol[0][var_name].imag()]))
688
sol_pts.append(point([a_sol[-1][var_name].real(), a_sol[-1][var_name].imag()]))
689
return sum(sol_pts) + sum(path_lines)
690
else:
691
return sum(path_lines)
692
693
def mixed_volume(self, polys, verbose=False):
694
"""
695
Computes the mixed volume of the polynomial system given by the input polys.
696
697
INPUT:
698
polys -- a list of multivariate polynomials (elements of a multivariate
699
polynomial ring).
700
verbose -- print lots of verbose information about what this function does.
701
702
OUTPUT:
703
The mixed volume.
704
705
EXAMPLES:
706
sage: from sage.interfaces.phc import * #optional
707
sage: R2.<x,y,z> = PolynomialRing(QQ,3) #optional
708
sage: test_sys = [(x+y+z)^2-1,x^2-x,y^2-1] #optional
709
sage: phc.mixed_volume(test_sys) #optional
710
4
711
"""
712
output_filename = self._output_from_command_list(['phc -m','4','n','n','n'], polys, verbose = verbose)
713
714
out = open(output_filename).read()
715
# All done
716
out_lines = out.split('\n')
717
for a_line in out_lines:
718
# the two conditions below are necessary because of changes in output format
719
if a_line.find('The mixed volume equals :') == 0 or a_line.find('common mixed volume :') == 0:
720
if verbose: print 'found line: ' + a_line
721
mixed_vol = Integer(a_line.split(':')[1])
722
break
723
724
try:
725
return mixed_vol
726
except NameError:
727
raise RuntimeError, "Mixed volume not found in output; something went wrong running phc."
728
729
730
def start_from(self, start_filename_or_string, polys, input_ring, path_track_file = None, verbose = False):
731
"""
732
This computes solutions starting from a phcpack solution file.
733
734
INPUT:
735
start_filename_or_string -- the filename for a phcpack start system, or the contents of such a file as a string. Variable names must match the inputring variables. The value of the homotopy variable t should be 1, not 0.
736
polys -- a list of multivariate polynomials (elements of a multivariate
737
polynomial ring).
738
input_ring: for coercion of the variables into the desired ring.
739
path_track_file: whether to save path-tracking information
740
verbose -- print lots of verbose information about what this function does.
741
742
OUTPUT:
743
A solution in the form of a PHCObject.
744
745
EXAMPLES:
746
sage: from sage.interfaces.phc import * #optional
747
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
748
sage: start_sys = [x^6-y^2,y^5-1] #optional
749
sage: sol = phc.blackbox(start_sys, R2) #optional
750
sage: start_save = sol.save_as_start() #optional
751
sage: end_sys = [x^7-2,y^5-x^2] #optional
752
sage: sol = phc.start_from(start_save, end_sys, R2) #optional
753
sage: len(sol.solutions()) #optional
754
30
755
"""
756
input_filename = sage.misc.misc.tmp_filename()
757
output_filename = sage.misc.misc.tmp_filename()
758
759
if start_filename_or_string.find('THE SOLUTIONS') != -1:
760
start_filename = sage.misc.misc.tmp_filename()
761
start_file = file(start_filename,'w')
762
start_file.write(start_filename_or_string)
763
start_file.close()
764
elif os.path.exists(start_filename_or_string):
765
start_filename = start_filename_or_string
766
else:
767
raise RuntimeError, "There is something wrong with your start string or filename"
768
769
# Get the input polynomial text
770
input = self._input_file(polys)
771
if verbose:
772
print "Writing the input file to %s"%input_filename
773
open(input_filename,'w').write(input)
774
775
if verbose:
776
print "The following file will be the input polynomial file to phc."
777
print input
778
779
# Create a phc process
780
child_phc = pexpect.spawn('phc')
781
child_phc.sendline('y')
782
child_phc.sendline(input_filename)
783
child_phc.sendline(output_filename)
784
child_phc.sendline('0')
785
child_phc.sendline('0')
786
child_phc.expect('Nonlinear Reduction')
787
child_phc.sendline('A')
788
child_phc.sendline(start_filename)
789
child_phc.sendline('y')
790
child_phc.sendline('1')
791
child_phc.sendline('0')
792
if verbose:
793
phc_dialog = child_phc.read(size = 40)
794
print phc_dialog
795
child_phc.sendline('n')
796
child_phc.sendline('0')
797
if verbose:
798
child_phc.expect('CURRENT CONTINUATION')
799
phc_dialog = child_phc.read(size = 40)
800
print phc_dialog
801
child_phc.sendline('0')
802
if path_track_file == None:
803
child_phc.sendline('0')
804
else:
805
child_phc.sendline('2')
806
child_phc.expect('results')
807
dots = child_phc.read()
808
if verbose:
809
print "should be . : " + dots
810
811
#close down the process:
812
child_phc.close()
813
if not os.path.exists(output_filename):
814
raise RuntimeError, "The output file does not exist; something went wrong running phc."
815
816
# Read the output produced by PHC
817
out = open(output_filename).read()
818
819
# Delete the temporary files
820
os.unlink(output_filename)
821
os.unlink(input_filename)
822
823
# All done
824
return PHC_Object(out, input_ring)
825
826
def blackbox(self, polys, input_ring, verbose = False):
827
"""
828
Returns as a string the result of running PHC with the given polynomials
829
under blackbox mode (the '-b' option).
830
831
INPUT:
832
polys -- a list of multivariate polynomials (elements of a multivariate
833
polynomial ring).
834
input_ring: for coercion of the variables into the desired ring.
835
verbose -- print lots of verbose information about what this function does.
836
837
OUTPUT:
838
a PHC_Object object containing the phcpack output string.
839
840
EXAMPLES:
841
sage: from sage.interfaces.phc import * #optional
842
sage: R2.<x,y> = PolynomialRing(QQ,2) #optional
843
sage: start_sys = [x^6-y^2,y^5-1] #optional
844
sage: sol = phc.blackbox(start_sys, R2) #optional
845
sage: len(sol.solutions()) #optional
846
30
847
"""
848
849
# Get three temporary file names (these will be in SAGE_HOME/.sage/tmp/pid)
850
input_filename = sage.misc.misc.tmp_filename()
851
output_filename = sage.misc.misc.tmp_filename()
852
log_filename = sage.misc.misc.tmp_filename()
853
854
# Get the input polynomial text
855
input = self._input_file(polys)
856
if verbose:
857
print "Writing the input file to %s"%input_filename
858
open(input_filename,'w').write(input)
859
860
if verbose:
861
print "The following file will be the input polynomial file to phc."
862
print input
863
864
# Create the phc command line>
865
cmd = 'phc -b %s %s'%(input_filename, output_filename)
866
867
if verbose:
868
print "The phc command line is:"
869
print cmd
870
871
# Do it -- make the system call.
872
e = os.system(cmd)
873
874
# Was there an error?
875
if e:
876
from sage.misc.sage_ostools import have_program
877
if not have_program('phc'):
878
print os.system('which phc') + ' PHC needs to be installed and in your path'
879
raise RuntimeError
880
# todo -- why? etc.
881
raise RuntimeError, open(log_filename).read() + "\nError running phc."
882
883
if not os.path.exists(output_filename):
884
raise RuntimeError, "The output file does not exist; something went wrong running phc."
885
886
# Read the output produced by PHC
887
out = open(output_filename).read()
888
889
# Delete the temporary files
890
os.unlink(input_filename)
891
os.unlink(output_filename)
892
893
# All done
894
return PHC_Object(out, input_ring)
895
896
897
################################
898
899
# The unique phc interface instance.
900
phc = PHC()
901
902
903