factor graphs, burnside sampling

.#factor.py

@@ -0,0 +1 @@
+sholtzen@wifi-131-179-3-172.host.ucla.edu.13952

evidence.inst

@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<instantiation date="Jun 4, 2005 7:07:21 AM">
+<inst id="0" value="false"/>
+</instantiation>

factor.py

@@ -0,0 +1,106 @@
+from sage.all import *
+import numpy.random
+import random
+from my_graphs import *
+
+def flip():
+    return numpy.random.binomial(1, 0.5)
+
+class FactorGraph:
+    ### graph: a sage graph
+    ### variables: a list of graph vertices which correspond to variables in the factor graph
+    ### factors: a list of list of graph vertices which correspond to colored factors in the factor graph
+    ### potential: state -> real: a function which evaluates the potential on a particular state
+    ###    a state is a dictionary assigning variables to Boolean values
+    def __init__(self, graph, variables, factors, potential):
+        self.graph = graph
+        self.variables = variables
+        self.factors = factors
+        self.potential = potential
+
+    ### converts a variable state into a variable partition
+    def state_to_partition(self, state):
+        var_part_1 = []
+        var_part_2 = []
+        for v,val in state.iteritems():
+            if val == True:
+                var_part_1.append(v)
+            else:
+                var_part_2.append(v)
+        return [var_part_1, var_part_2]
+
+    ### computes a burnside process transition beginning from a state particular state
+    ### n: number of moves
+    def burnside(self, state, n):
+        print("state: %s" % state)
+        for j in range(0,n):
+            var_part = self.state_to_partition(state)
+            partition = var_part + self.factors
+            stab = self.graph.automorphism_group(partition=partition)
+            stab = libgap(stab)
+            # product replacement
+            p = libgap.PseudoRandom(stab).sage()
+            # now convert p into a state
+            state = dict()
+            for cyc in p.to_cycles():
+                # sample a value for cyc and fill in the new current value
+                v = flip()
+                for idx in cyc:
+                    if (idx - 1) in self.variables:
+                        #-1, zero indexed
+                        state[idx-1] = v
+        return state
+
+    ### perform a single step of orbit jumping
+    ### returns: a pair, (the ratio of transition probabilities, new state)
+    ### n: number of burnside steps to take
+    def orbitjump(self, state, n):
+        hatx = self.burnside(state, n)
+        probx = self.potential(state)
+        probhatx = self.potential(hatx)
+        orbx = self.graph.automorphism_group(partition=self.state_to_partition(state)).order()
+        orbhatx = self.graph.automorphism_group(partition=self.state_to_partition(hatx)).order()
+        transitionprob = (probhatx * orbx) / (probx * orbhatx)
+        return (transitionprob, hatx)
+
+    ### draw n samples according to orbit jump MCMC with no orbital MCMC
+    ### burnsidesize: number of burnside steps to take
+    ### returns a list of samples
+    def orbitjumpmcmc(self, burnsidesize, n):
+        samples = []
+        # set up initial random state
+        cur_state = dict()
+        for v in self.variables:
+            cur_state[v] = flip()
+
+        for i in range(0, n):
+            (ratio, new_state) = self.orbitjump(cur_state, burnsidesize)
+            # accept with transition probability
+            acceptprob = numpy.minimum(1, ratio)
+            if numpy.random.binomial(1, acceptprob):
+                samples.append(new_state)
+                cur_state = new_state
+            else:
+                samples.append(cur_state)
+
+        return samples
+
+
+
+def gen_complete_pairwise_factorgraph(n):
+    (g, (v, factors)) = gen_complete_pairwise_factor(n)
+    def potential(state):
+        print("potential got state: %s" % state)
+        p = 0.0
+        for v in state.itervalues():
+            if v:
+                p += 1
+        return p
+    return FactorGraph(g, v, [factors], potential)
+
+def main():
+    graph = gen_complete_pairwise_factorgraph(4)
+    print(graph.orbitjumpmcmc(4, 10))
+
+if __name__ == "__main__":
+    main()

gen_pairwise_fg.py

@@ -0,0 +1,24 @@
+import sys
+import itertools
+
+def findsubsets(S,m):
+    return set(itertools.combinations(S, m))
+
+
+n = 24
+print("MARKOV")
+print(n)
+for i in range(0,n):
+    sys.stdout.write("2 " )
+print("")
+# print factor description
+subsets = findsubsets(range(0,n), 2)
+print(len(subsets))
+for (a,b) in subsets:
+    print("2 %s %s " % (a,b))
+print("")
+# print tables
+for i in range(0, len(subsets)):
+    print("4")
+    print("  0.2 0.3")
+    print("  0.3 0.2")

orbitgen.py

@@ -53,6 +53,7 @@ def bfs_foldrep(graph, colors, fixcolors, acc, f, orders=False):
     # set of representative graph colorings
     reps = set()
     while len(queue) != 0:
+        # print("queue: %s" % queue)
         c = queue.popleft()
         # TODO: turn this into a single GI call by having it return both the
         # automorphism group and the canonical fora
@@ -73,7 +74,7 @@ def bfs_foldrep(graph, colors, fixcolors, acc, f, orders=False):
         if (gcanon, (tuple(c_canon[0]), tuple(c_canon[1]))) in reps:
             continue
 
-        A, orbits = graph.automorphism_group(partition=c + fixcolors, orbits=True,
+        A, orbits = graph.automorphism_group(partition=part, orbits=True,
                                              algorithm="bliss")
         if orders:
             acc = f(acc, graph, c, g_order, A.order())
@@ -157,4 +158,5 @@ def find_representatives():
 
 
 if __name__ == "__main__":
+    set_gap_memory_pool_size(900000000)
     find_representatives()

query.py

@@ -59,7 +59,7 @@ def mpe(G, variables, factors, potential):
     print("G: %s, vars: %s, factors: %s" % (G, variables, factors))
     def max_prob(acc, g, color):
         (cur_max_state, cur_max_value) = acc
-        if len(color[0]) < len(g.vertices()) / 2.0:
+        if len(color[0]) < len(variables) / 2.0:
             # check both possibilities
             pot = potential(color[::-1])
             if cur_max_value < pot:
@@ -77,48 +77,51 @@ def max_prob(acc, g, color):
 # normalizing constant on a particular orbit
 def prob(G, variables, factors, potential):
     # folding function
-    print("G: %s, vars: %s, factors: %s" % (G, variables, factors))
     def sum_prob(acc, g, color, g_order, cur_order):
         cur_sz = g_order / cur_order # yay orbit stabilizer theorem!
         (cur_prob, cur_z) = acc
-        if len(color[0]) < len(g.vertices()) / 2.0:
+        if len(color[0]) < len(variables) / 2.0:
             # check both possibilities
             (prob, z) = potential(color[::-1])
             cur_prob += prob * cur_sz
             cur_z += z * cur_sz
-        (prob, z) = potential(color[::-1])
+        (prob, z) = potential(color)
         cur_prob += prob * cur_sz
         cur_z += z * cur_sz
         return (cur_prob, cur_z)
     return bfs_foldrep(G, variables, factors, (0.0, 0.0), sum_prob, orders=True)
 
 
-def compute_mpe_complete_2factor():
-    G = gen_complete_pairwise_factor(100)
+def compute_prob_complete_2factor():
+    G = gen_complete_pairwise_factor(55)
     # add evidence: the first variable is false
     G[0].add_vertex(name="e")
     G[0].add_edge(("e", G[1][0][0]))
     def potential(assgn):
+        # print("potential eval: %s" % assgn)
         p = 0.0
         # check for evidence
-        if G[1][0][0] in assgn[0]:
-            return 0.0
         for factor in G[1][1]:
             connected = G[0].neighbors(factor)
-            if connected[0] != connected[1]:
+            val1 = connected[0] in assgn[0]
+            val2 = connected[1] in assgn[1]
+            if val1 != val2:
                 p += 0.3
             else:
                 p += 0.2
             for v in connected:
                 if v in assgn[0]:
                     p += 1
-        return p
+        if G[1][0][0] in assgn[0]:
+            return (0.0, p)
+        else:
+            return (p, p)
 
     pr = cProfile.Profile()
     pr.enable()
 
-    res = mpe(G[0], G[1][0], [G[1][1] + ["e"]], potential)
-    print("max state: %s, max value: %s" % (res[0], res[1]))
+    res = prob(G[0], G[1][0], [G[1][1], ["e"]], potential)
+    print("prob: %f" % (res[0] / res[1]))
     # print("brute forced: %d" % (bruteforce_partition(G, partfun)))
 
     pr.disable()
@@ -129,7 +132,7 @@ def potential(assgn):
     print s.getvalue()
 
 def compute_prob_big_factor():
-    G = gen_single_big_factor(10)
+    G = gen_single_big_factor(2)
     # add evidence: the first variable is false
     G[0].add_vertex(name="e")
     G[0].add_edge(("e", G[1][0][0]))
@@ -146,7 +149,7 @@ def potential(assgn):
     pr = cProfile.Profile()
     pr.enable()
 
-    res = prob(G[0], G[1][0], [G[1][1] + ["e"]], potential)
+    res = prob(G[0], G[1][0], [G[1][1], ["e"]], potential)
     print("prob %f" % (res[0] / res[1]))
     # print("brute forced: %d" % (bruteforce_partition(G, partfun)))
 
@@ -157,6 +160,68 @@ def potential(assgn):
     ps.print_stats()
     print s.getvalue()
 
+def compute_mpe_big_factor():
+    G = gen_single_big_factor(300)
+    # add evidence: the first variable is false
+    G[0].add_vertex(name="e")
+    G[0].add_edge(("e", G[1][0][0]))
+    def potential(assgn):
+        p = 0.0
+        for v in assgn[0]:
+            p += 1
+
+        if G[1][0][0] in assgn[0]:
+            return 0.0
+        else:
+            return p
+
+    pr = cProfile.Profile()
+    pr.enable()
+
+    res = mpe(G[0], G[1][0], [G[1][1], ["e"]], potential)
+    print("mpe %s, %s" % (res[0], res[1]))
+    # print("brute forced: %d" % (bruteforce_partition(G, partfun)))
+
+    pr.disable()
+    s = StringIO.StringIO()
+    sortby = 'cumulative'
+    ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
+    ps.print_stats()
+    print s.getvalue()
+
+def compute_mpe_full_pairwise_factor():
+    G = gen_complete_pairwise_factor(50)
+    # add evidence: the first variable is false
+    G[0].add_vertex(name="e")
+    G[0].add_edge(("e", G[1][0][0]))
+    def potential(assgn):
+        p = 0.0
+        for v in assgn[0]:
+            p += 1
+
+        if G[1][0][0] in assgn[0]:
+            return 0.0
+        else:
+            return p
+
+    pr = cProfile.Profile()
+    pr.enable()
+
+    res = mpe(G[0], G[1][0], [G[1][1], ["e"]], potential)
+    print("mpe %s, %s" % (res[0], res[1]))
+    # print("brute forced: %d" % (bruteforce_partition(G, partfun)))
+
+    pr.disable()
+    s = StringIO.StringIO()
+    sortby = 'cumulative'
+    ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
+    ps.print_stats()
+    print s.getvalue()
+
+
 
 if __name__ == "__main__":
-    compute_prob_big_factor()
+    set_gap_memory_pool_size(5000000)
+    compute_prob_complete_2factor()
+    # compute_prob_big_factor()
+    # compute_mpe_big_factor()