# -*- coding: utf-8 -*- # Runs the Karp-Luby estimator import bisect; import random; # Variables and probabilities p = { 'X1': 0.09, 'X2': 0.8, 'Y1': 0.85, 'Y2': 0.75, 'Y3': 0.065 } # The DNF formula # Literal = array [positive?, variable] # Clause = array of literals (conjunctive) # Formula = array of clauses (disjunctive) phi = [ [ [True, 'X1'], [True, 'Y1'] ], [ [True, 'X1'], [True, 'Y2'] ], [ [True, 'X2'], [True, 'Y3'] ] ] m = len(phi) # compute clause probabilities q = [] for Ci in phi: qi = 1 for positive,var in Ci: if positive: qi *= p[var] else: qi *= 1-p[var] q.append(qi) # sum them up (to facilitate clause sampling) qcum = q for i in range(m-1): qcum[i+1] += qcum[i] Q = qcum[-1] # checks whether conjunct i is satisfied def issat(theta, i): Ci = phi[i] for positive,var in Ci: if positive != theta[var]: return False; return True # checks whether the conjunct i is the first satisfied def isfirstsat(theta, i): for j in range(i): if issat(theta, j): return False return issat(theta, i) # Number of KL steps to run n = 10000 # Array of results Z = [] # Current valuation (initialized empty) theta = { } # Run MC for k in range(n): # generate theta for var,pvar in p.iteritems(): theta[var] = random.random() < pvar # pick a clause u = random.random()*Q i = bisect.bisect(qcum, u) # make the clause satisfied Ci = phi[i] for positive,var in Ci: theta[var] = positive # compute Z_k Z.append(isfirstsat(theta, i)) # print information about current step print 'k : ', k+1 print 'i, Ci : ', i+1, phi[i] print 'theta_k: ', theta print 'Z_k : ', Z[k] print 'phat : ', Q/(k+1.)*sum(Z) print '---'