COPYING file for distribution
[anf2cnf:anf2cnf.git] / anf2cnf.py
1 # -*- coding: utf-8 -*-
2 r"""
3 Boolean Polynomial SAT-Solver
4
5 Given an ideal or polynomial system this module performs conversion to
6 the DIMACS CNF format, calls MiniSat2 on that input and parses the
7 output.
8
9 AUHTORS:
10
11 - Martin Albrecht - (2008-09) initial version
12 - Mate SOOS - (2010) updated version
13 """
14
15 import commands
16 import tempfile
17 import re
18
19 from sage.rings.polynomial.pbori import BooleanPolynomial, BooleanMonomial
20 from sage.misc.prandom import shuffle as do_shuffle
21
22 from sage.all import *
23
24 class ANFSatSolver(SageObject):
25     """
26     Solve a boolean polynomial system using MiniSat2.
27     """
28     def __init__(self, ring=None, cutting_num=None):
29         """
30         Setup the SAT-Solver and reset internal data.
31
32         This function is also called from :meth:`__call__()` to pass
33         in parameters.
34
35         INPUT:
36
37         - ``ring`` - a boolean polynomial ring
38
39         - ``cuting_num`` - the cutting number ``>= 2`` (default: ``4``)
40
41         EXAMPLE::
42
43             sage: B = BooleanPolynomialRing(10,'x')
44             sage: ANFSatSolver(B)
45             ANFSatSolver(4) over Boolean PolynomialRing in x0, x1, x2, x3, x4, x5, x6, x7, x8, x9
46
47         """
48         self._i = 1 # the maximal index for literals
49         self.minus = -1
50         self._mon_map = []
51
52         if hasattr(self,"_ring") and self._ring is not None:
53             if ring is not None:
54                 assert(ring is self._ring)
55         else:
56             assert(ring is not None)
57             self._ring = ring
58
59         if hasattr(self,"cutting_num") and cutting_num is None:
60             pass
61         elif cutting_num is None:
62             self.cutting_num = 4
63         else:
64             if cutting_num < 2:
65                 raise TypeError("cutting_num must be >= 2 but is %d."%(cutting_num,))
66             self.cutting_num = cutting_num
67
68         self.cnf_literal_map.clear_cache()
69
70     def _repr_(self):
71         return "ANFSatSolver(%d) over %s"%(self.cutting_num,self._ring)
72
73     def __getattr__(self, name):
74         if name == 'ring':
75             return self._ring
76         else:
77             raise AttributeError("ANFSatSolver does not have an attribute names '%s'"%name)
78
79     def _cnf_literal(self):
80         """
81         Return a new variable for the CNF formulas.
82
83         EXAMPLE::
84
85             sage: B.<x,y,z> = BooleanPolynomialRing()
86             sage: anf2cnf = ANFSatSolver(B)
87             sage: anf2cnf._cnf_literal()
88             2
89             sage: anf2cnf._cnf_literal()
90             3
91
92         .. note::
93
94           Increases the internal literal counter.
95         """
96         self._i += 1
97         return self._i-1
98     
99     @cached_method
100     def cnf_literal_map(self, m):
101         """
102         Given a monomial ``m`` in a boolean polynomial ring return a
103         tuple ``((i,),(c0,c1,...))`` where ``i`` is the index of the
104         monomial ``m`` and ``c0,c1,...`` are clauses which encode the
105         relation between the monomial ``m`` and the variables
106         contained in ``m``.
107
108         EXAMPLE::
109
110             sage: B.<x,y,z> = BooleanPolynomialRing()
111             sage: anf2cnf = ANFSatSolver(B)
112             sage: anf2cnf.cnf_literal_map(B(1))
113             1
114
115             sage: anf2cnf.cnf_literal_map(x)
116             2
117
118             sage: anf2cnf.cnf_literal_map(x*y)
119             (4, [(2, -4), (3, -4), (4, -2, -3)])
120
121             sage: anf2cnf.cnf_literal_map(y)
122             3
123
124         .. note:: 
125         
126           May call :meth:`_cnf_literal()`
127         """
128         minus = self.minus
129         cnf_literal_map = self.cnf_literal_map
130
131         if isinstance(m, BooleanPolynomial):
132             if len(m) == 1:
133                 m = m.lm()
134             else:
135                 raise TypeError("Input must be monomial.")
136         elif not isinstance(m, BooleanMonomial):
137             raise TypeError("Input must be of type BooleanPolynomial.")
138
139         if m.deg() == 0:
140             # adding the clause that 1 has to be True
141             monomial = self._cnf_literal()
142             return monomial
143
144         elif m.deg() == 1:
145             # a variable
146             monomial = self._cnf_literal()
147             return monomial
148
149         else:
150             # we need to encode the relationship between the monomial
151             # and its variables
152             variables = [cnf_literal_map(v) for v in m.variables()]
153             monomial = self._cnf_literal()
154
155             # (a | -w) & (b | -w) & (w | -a | -b) <=> w == a*b
156             this_mon_map = []
157             for v in variables:
158                 this_mon_map.append((v, minus * monomial))
159             this_mon_map.append((monomial,) + sum([(minus*v,) for v in variables],tuple()))
160             self._mon_map.append([this_mon_map, ("monomial %s" % m)])
161
162             return monomial
163
164     @cached_function
165     def cached_permutations(mylen, equal_zero):
166         """
167         Generates permutations for XOR clauses, such as:
168         -1  1  1
169          1 -1  1
170          1  1 -1
171         -1 -1 -1
172
173         or, for the inverse (inverted 'equal_zero'):
174         -1 -1  1
175         -1  1 -1
176          1 -1 -1
177          1  1  1
178         """
179         E = []
180         for numNegated in range(0, mylen+1) :
181             if (((numNegated % 2) ^ ((mylen+1) % 2)) == equal_zero) : continue
182             start = []
183             for i in range(numNegated) :
184                 start.append(1)
185             for i in range(mylen-numNegated) :
186                 start.append(-1)
187             E.extend(Permutations(start))
188         return E
189
190     def _get_monomial_map(self) :
191         """
192         Returns the variable-to-monomial map. These essentially
193         are the clauses that define the monomials
194         """
195         return self._mon_map
196
197     def cnf(self, F,  shuffle=False, format='dimacs', **kwds):
198         """
199         Return CNF for ``F``.
200
201         INPUT:
202
203         - ``shuffle`` - shuffle output list (default: ``False``)
204
205         - ``format`` - either 'dimacs' for DIMACS or ``None`` for
206           tuple list (default: ``dimacs``)
207
208         EXAMPLE:
209             sage: B.<a,b,c,d> = BooleanPolynomialRing()
210             sage: solver = ANFSatSolver(B)
211             sage: print solver.cnf([a*b + c + 1, d + c, a + c])
212             p cnf 5 9
213             c ------------------------------
214             c Next definition: a*b + c + 1
215             -3 -4 0
216             3 4 0
217             c ------------------------------
218             c Next definition: c + d
219             4 -5 0
220             -4 5 0
221             c ------------------------------
222             c Next definition: a + c
223             1 -4 0
224             -1 4 0
225             c ------------------------------
226             c Next definition: a*b
227             1 -3 0
228             2 -3 0
229             3 -1 -2 0
230         """
231         self.__init__( **kwds)
232
233         if get_verbose() >= 1:
234             print "Parameters: c: %d, shuffle: %s"%(self.cutting_num,shuffle)
235
236         C = []
237         for f in F:
238             this_c = []
239             for c in self._process_polynomial(f):
240                 this_c.append(c)
241             C.append([this_c, ("%s" % f)])
242
243         for mono_def,mono_desc in self._mon_map :
244             C.append([mono_def, mono_desc])
245
246         if shuffle:
247             do_shuffle(C)
248
249         if format == 'dimacs':
250             return self.to_dimacs(C)
251         else:
252             return C
253
254     def delimited(self, F):
255         out = ""
256         for f in F:
257             for c in f:
258                 for m in c:
259                     if (m != self._ring(1).lm()) :
260                         for v in m.variables() :
261                             out += "%d," %(v.index()+1)
262                 out += "|"
263             out += "\n"
264
265         return out
266
267     def solveDelimited(self, F):
268         tmpf, tmpfname = tempfile.mkstemp(text = True)
269         os.write(tmpf, self.delimited(F))
270         os.close(tmpf)
271
272         tmpf2, tmpfname2 = tempfile.mkstemp(text = True)
273         os.close(tmpf2)
274         toexec = "SATsolve %s > %s" %(tmpfname, tmpfname2)
275         os.system(toexec)
276         os.unlink(tmpfname)
277
278         f = open(tmpfname2, "r")
279         text = f.read()
280         mylines = text.splitlines()
281         f.close()
282
283         sol = {}
284         time = 0.0
285         gen = self._ring.gens()
286         foundTime = False
287         foundSol = False
288         for line in mylines:
289             line = line.rstrip().lstrip()
290             if (re.search('CPU time', line) != None):
291                 time = float(line.lstrip('c CPU time').lstrip().lstrip(':').rstrip('s').strip())
292                 foundTime = True
293                 continue
294             if (line[0] == "[" and len(line) != 1):
295                 num = 0
296                 for v in line.split() :
297                     if (v == "[" or v == "]"): continue
298                     v = v.strip(',')
299                     sol[gen[num]] = int(v)
300                     num += 1
301                 foundSol = True
302                 continue
303
304         assert foundTime, "Didn't find time that took to solve the instance! File: %s" %(tmpfname2)
305         assert foundSol, "Didn't find solution!"
306         os.unlink(tmpfname2)
307
308         return sol, time
309
310
311     def to_dimacs(self, C):
312         """
313         Prints header + clause descriptions + clauses
314
315         Calculates number of clauses and number of variables to print the
316         P VARS CLAUSES
317         header. Then it just prints the clauses one after the other
318         """
319         numVars = self._i-1
320         nClauses = 0
321         for clause_set, desc in C :
322             nClauses += len(clause_set)
323
324         out = ["p cnf %d %d\n"%(numVars,nClauses)]
325         for clause_set, desc in C :
326             out.append("c ------------------------------\n");
327             out.append("c Next definition: %s\n" %  desc);
328             out.extend(" ".join(map(str,clause)) + " 0\n"  for clause in clause_set)
329         return "".join(out)
330
331     def _process_polynomial(self, f):
332         """
333         EXAMPLE:
334             sage: B.<a,b,c,d> = BooleanPolynomialRing()
335             sage: solver = ANFSatSolver(B)
336             sage: solver._process_polynomial(a*b + c + 1)
337             [(2, -4), (3, -4), (4, -2, -3), (1,), (4, 5), (-4, -5)]
338         """
339         M, E = [], []
340         #M is the XOR
341         #E is the helper clauses (to define the monomials)
342         equal_zero = True
343         cnf_literal_map = self.cnf_literal_map
344         for m in f:
345             if (m == self._ring(1).lm()) :
346                 equal_zero = not equal_zero
347             else :
348                 mbar = cnf_literal_map( m )
349                 M.append( mbar )
350
351         if len(M) > self.cutting_num:
352             for Mbar, this_equal_zero in self._split_xor_list(M, equal_zero):
353                 E.extend(self._cnf_for_xor_list(Mbar, this_equal_zero))
354         else:
355             E.extend( self._cnf_for_xor_list(M, equal_zero) )
356
357         return E
358
359     def _split_xor_list(self, monomial_list, equal_zero):
360         """
361         Splits a list of monomials into sublists and introduces
362         connection variables. 
363
364         INPUT:
365
366         - ``monomial_list`` - a list of indices already registered
367           with ``self``.
368
369         EXAMPLE::
370
371             sage: B = BooleanPolynomialRing(3,'x')
372             sage: asolver = ANFSatSolver(B)
373             sage: l = [asolver._cnf_literal() for _ in range(5)]; l
374             [2, 3, 4, 5, 6]
375             sage: asolver._split_xor_list(l)
376             [[2, 3, 7], [7, 4, 5, 8], [8, 6]]
377         """
378         c = self.cutting_num
379
380         nm = len(monomial_list)
381         part_length =  ceil((c-2)/ZZ(nm) * nm)
382         M = []
383
384         new_variables = []
385         for j in range(0, nm, part_length):
386             m =  new_variables + monomial_list[j:j+part_length]
387             if (j + part_length) < nm:
388                 new_variables = [self._cnf_literal()]
389                 m += new_variables
390             M.append([m, equal_zero])
391             equal_zero = True
392         return M
393
394     def _cnf_for_xor_list(self, M, equal_zero):
395         minus = self.minus
396
397         E = []
398         ll = len( M )
399         for p in self.cached_permutations(ll, equal_zero):
400             E.append( sum([(p[i] * M[i],) for i in range(ll)],tuple()) )
401         return E
402             
403     def __call__(self, F, **kwds):
404         """
405         EXAMPLE:
406             sage: sr = mq.SR(1,1,1,4,gf2=True)
407             sage: F,s = sr.polynomial_system()
408             sage: B = BooleanPolynomialRing(F.ring().ngens(), F.ring().variable_names())
409             sage: F = [B(f) for f in F if B(f)]
410             sage: solver = ANFSatSolver(B)
411             sage: solution, t = solver(F)
412             sage: solution
413             {k001: 0, k002: 0, s003: 1, k000: 1, k003: 0, 
414              x103: 0, w100: 0, w101: 0, w102: 1, s000: 1, 
415              w103: 1, s002: 1, s001: 1, x102: 1, x101: 1, 
416              x100: 1, k103: 0, k101: 0, k102: 0, k100: 1}
417
418             sage: B.ideal(F).groebner_basis()
419             [k100 + k003 + 1, 
420              k101 + k003, 
421              k102, 
422              k103 + k003, 
423              x100 + 1, 
424              ...
425              k000 + k003 + 1, 
426              k001, 
427              k002 + k003]
428         """
429         fn = tmp_filename()
430         have_poly = False
431
432         if isinstance(F, str):
433             fh = open(fn,"w")
434             fh.write(self.cnf(F, format='dimacs', **kwds))
435             fh.close()
436
437         else:
438             p = iter(F).next()
439             if isinstance(p, BooleanPolynomial):
440                 have_poly = True
441                 self._ring = p.parent()
442
443                 fh = open(fn,"w")
444                 fh.write(self.cnf(F, format='dimacs', **kwds))
445                 fh.close()
446             elif isinstance(p, tuple):
447                 fh = open(fn,"w")
448                 fh.write(self.to_dimacs(F))
449                 fh.close()
450             else:
451                 raise TypeError("Type '%s' not supported."%(type(p),))
452
453         # call MiniSat
454         on = tmp_filename()
455         s =  commands.getoutput("minisat2 %s %s"%(fn,on))
456
457         if get_verbose() >= 2:
458             print 
459             print 
460             print s
461
462         s = s.splitlines()
463         for l in s:
464             if "CPU time" in l:
465                 t = float(l[l.index(":")+1:l.rindex("s")])
466                 break
467
468         res =  open(on).read()
469         if res.startswith("UNSAT"):
470             return False, t
471         res = res[4:]
472         res = map(int, res.split(" "))
473
474         # parse result
475         if not have_poly:
476             return res, t
477         else:
478             return self.map_solution_to_variables(res), t
479
480     def map_solution_to_variables(self, res, gd=None):
481         """
482         """
483         if gd is None:
484             gens = self._ring.gens()
485             gd = {}
486             cnf_literal_map = self.cnf_literal_map
487             for gen in gens:
488                 try:
489                     im = cnf_literal_map(gen.lm())
490                     gd[im] = gen
491                 except KeyError:
492                     pass
493
494         solution = {}
495         for r in res:
496             if abs(r) in gd:
497                 if r>0:
498                     solution[gd[abs(r)]] = 1
499                 else:
500                     solution[gd[abs(r)]] = 0
501         return solution
502
503 def beta(F):
504     B = F.ring()
505     
506     n = F.nvariables()
507     m = F.ngens()
508     k = sum([len(f) for f in F])
509     d = max([f.total_degree() for f in F])
510     return k / float(m * sum([binomial(n,i) for i in range(0,d+1)]))
511
512