Implementation of 3SAT to SUBSETSUM reduction (and also via 3D Matching)

reduce_to_change.py

## Factoring to 3-SAT
## https://www.cs.indiana.edu/cgi-pub/sabry/cnf.html

import math

class ThreeSatInstance(object):
    def __init__( self, numVariables ):
        self.numVariables = numVariables
        self.clauses = []

    def addClauseFromFile( self, line ):
        tokens = line.split( " " )
        clause = tuple( int( x ) for x in tokens )
        assert clause[-1] == 0
        self.clauses.append( clause[:-1] )

    def addClause( self, x, y, z ):
        self.clauses.append( (x, y, z) )

    def dump( self ):
        print "|vars| = ", self.numVariables
        print "|clauses| =", len( self.clauses )
        for c in self.clauses:
            print "  ", " ".join( str(x) for x in c )
        
def readDimacs( filename ):
    tsi = None
    with open( filename, "r" ) as f:
        for line in f:
            if line[0] == "c":
                continue
            elif line[0] == "p":
                tokens = line.split( " " )
                assert tokens[1] == "cnf"
                numVars = int( tokens[2] )
                numClauses = int( tokens[3] )
                tsi = ThreeSatInstance( numVars )
            elif line.strip() == "":
                continue
            else:
                tsi.addClauseFromFile( line )
                
    assert len( tsi.clauses ) == numClauses
    return tsi

def simple3SAT():
    sat = ThreeSatInstance( 5 )
    sat.addClause( 1, -3, -5 )
    sat.addClause( 2, 4, 5 )
    #sat.addClause( 2, -3, -1 )
    return sat
    
# Following the reduction here:
# https://www.cs.cmu.edu/~ckingsf/bioinfo-lectures/3dm.pdf
class ThreeDimensionalMatching(object):
    def __init__( self ):
        self.X = []
        self.Y = []
        self.Z = []
        self.triples = []
        self.variableToTips = {}
        self.nextClauseGadget = 1

    def addVariableGadget( self, i, numClauses ):
        # The slide give "a_ij" twice, the second should be "b_ij"
        # A_i = a_i1, a_i2, ..., a_{i,2k}  where k is the number of clauses
        # a_ij goes in X for even j
        # a_ij goes in Y for odd j
        # B_i = b_i1, b_i2, ..., b_{i,2k}
        # b_i goes in Z
        # Triples are ( a_ij, a_ij+1, b_ij )
        core = [ "a_" + str( i ) + "," + str( j )
                 for j in xrange( 1, 2 * numClauses + 1 ) ]
        tips = [ "b_" + str( i ) + "," + str( j )
                 for j in xrange( 1, 2 * numClauses + 1 ) ]

        for j in xrange( 0, 2 * numClauses ):
            a = core[j]
            b = core[(j+1)%(2 * numClauses)]
            if j % 2 == 0:
                # Odd core elements (even indexes in the list)
                self.Y.append( a )
                self.triples.append( (b,a,tips[j]) )
            else:
                self.X.append( a )
                self.triples.append( (a,b,tips[j]) )
            self.Z.append( tips[j] )
    
        self.variableToTips[i] = tips
            
    def addClauseGadget( self, clause ):
        j = self.nextClauseGadget
        self.nextClauseGadget += 1
        nodes = ( "c_" + str( j ), "c'_" + str( j ) )
        self.X.append( nodes[0] )
        self.Y.append( nodes[1] )

        for v in clause:
            if v < 0:
                # Not x_i
                # Pick tip b_i,1 for clause 1
                #          b_i,3 for clause 2
                vi = -v
                whichTip = 2 * j - 1
            else:
                # x_i
                # Pick tip b_i,2 for clause 1
                #      tip b_i,4 for clause 2
                vi = v
                whichTip = 2 * j
                
            tip = self.variableToTips[vi][whichTip-1]
            self.triples.append( ( nodes[0], nodes[1], tip ) )

    def finishCoveringTips( self, numClauses, numVars ):
        # After core is covered, numClauses * numVars tips remain open
        # The clause gadgets cover numClauses of them.
        cleanupGadgetsNeeded = numClauses * numVars - numClauses
        for i in xrange( 1, cleanupGadgetsNeeded + 1 ):
            x = "q_" + str( i )
            y = "q'_" + str( i )
            self.X.append( x )
            self.Y.append( y )
            for tipList in self.variableToTips.values():
                for b in tipList:
                    self.triples.append( (x, y, b) )
                    
        assert len( self.X ) == len( self.Y )
        assert len( self.Y ) == len( self.Z )
        
    def dump( self ):
        print "|T| =", len( self.triples )
        for t in self.triples:
            print( "{:10} {:10} {:10}".format( t[0], t[1], t[2] ) )
        print "|X| =", len( self.X )
        print "|Y| =", len( self.Y )
        print "|Z| =", len( self.Z )
            
            
def SAT_to_TDM( sat ):
    tdm = ThreeDimensionalMatching()
    for i in xrange( 1, sat.numVariables + 1 ):
        tdm.addVariableGadget( i, len( sat.clauses ) )
    for clause in sat.clauses:
        tdm.addClauseGadget( clause )
    tdm.finishCoveringTips( len( sat.clauses ), sat.numVariables )
    return tdm

# SAT -> SubsetSum reduction given here
# https://www.cs.mcgill.ca/~lyepre/pdf/assignment2-solutions/subsetSumNPCompleteness.pdf
class SubsetSum(object):
    def __init__( self ):
        self.w = []
        self.targetW = None
                  
    def zeros( self ):
        return [0] * (self.m + self.n)
                  
    def initFromSAT( self, numClauses, numVariables ):
        # Maximum column sum is 5 (three variables in a clause, plus the
        # two auxiliary columns.)
        # But, 5 doesn't produce a valid result base 4, it's 1 + a carry,
        # we only care if it's 3 + a carry
        self.d = 4
        self.n = numVariables
        self.m = numClauses

        self.tDigits = []
        self.fDigits = []
        for i in xrange( 1, numVariables + 1 ):
            # i'th digit =1
            ti = self.zeros()
            ti[i-1] = 1
            self.tDigits.append( ti )
                  
            fi = self.zeros()
            fi[i-1] = 1
            self.fDigits.append( fi )

        self.xDigits = []
        self.yDigits = []
        for j in xrange( 1, numClauses + 1 ):
            # n + j'th digit  = 1
            xj = self.zeros()
            xj[ self.n + j - 1 ] = 1
            # these won't be modified so they can be shared
            self.xDigits.append( xj )
            self.yDigits.append( xj )
        self.j = 1
                
    def addClause( self, c ):
        # Handling clause j, marking the appropriate variables
        j = self.j
        for v in c:
            # Variable v is at offset v-1
            # it's negated if we want it false
            if v < 0:
                self.fDigits[-v - 1][self.n + j-1] = 1
            else:
                self.tDigits[v-1][self.n + j-1] = 1
                  
        self.j += 1

    def digitsToNumber( self, digits ):
        base = self.d
        ret = 0
        for (p, x) in enumerate( reversed( digits ) ):
            ret += x * ( base ** p )
        return ret        
                  
    def finishSAT( self ):
        for i in xrange( self.n ):
            self.w.append( self.digitsToNumber( self.tDigits[i] ) )
            self.w.append( self.digitsToNumber( self.fDigits[i] ) )

        for j in xrange( self.m ):
            self.w.append( self.digitsToNumber( self.xDigits[j] ) )
            self.w.append( self.digitsToNumber( self.yDigits[j] ) )

        targetDigits = [1] * self.n + [3] * self.m
        self.targetW = self.digitsToNumber( targetDigits )
                  
    def initFromTdm( self, d, n ):
        self.d = d
        self.n = n
        # W = sum_{i=0}^{3n-1} d^i
        # which is just a geometric series
        self.targetW = ( d ** (3*n) - 1 ) / ( d - 1 )
        
    def assignOffset( self, X, Y, Z ):
        self.X = dict( (x,i+2*self.n) for (i,x) in enumerate( X ) )
        self.Y = dict( (y,i+self.n) for (i,y) in enumerate( Y ) )
        self.Z = dict( (z,i) for (i,z) in enumerate( Z ) )

    def addTriple( self, x, y, z ):
        self.w.append( self.d**(self.X[x]) +
                       self.d**(self.Y[y]) +
                       self.d**(self.Z[z]) )

    def dump( self ):
        numDigits = int( math.log( max( self.w ), 10 ) ) + 1
        fmtString = " {:>" + str( numDigits ) + "}"
        
        print "W =", self.targetW
        print "|w| =", len( self.w )
        print "w = {"
        for x in self.w:
            print fmtString.format( x )
        print "}"

def TDM_to_SUBSETSUM( tdm ):
    d = len( tdm.triples) + 1
    n = len( tdm.X )
    ssum = SubsetSum()
    ssum.initFromTDM( d, n )
    ssum.assignOffset( tdm.X, tdm.Y, tdm.Z )
    for (x,y,z) in tdm.triples:
        ssum.addTriple( x, y, z )
    return ssum

def SAT_to_SUBSETSUM( sat ):
    ssum = SubsetSum()
    ssum.initFromSAT( len( sat.clauses ), sat.numVariables )
    for c in sat.clauses:
        ssum.addClause( c )
    ssum.finishSAT()
    return ssum
                                
def test():
    #sat = simple3SAT()
    sat = readDimacs( "factor-91-3sat.dimacs" )
    print "### satisfiability ###"
    sat.dump()

    if False:    
        tdm = SAT_to_TDM( sat )
        print "### 3-dimensional matching ###"
        tdm.dump()

        ssum = TDM_to_SUBSETSUM( tdm )
        print "### subset sum ###"
        ssum.dump()

    if True:
        ssum = SAT_to_SUBSETSUM( sat )
        print "### subset sum ###"
        ssum.dump()