mpe experiment 1

orbitgen.py

@@ -4,10 +4,46 @@
 from collections import deque
 import my_bliss
 
+### generate all orbits of `elements` under `gens`
+def generate_orbits(elements, gens):
+    vset = set(elements)
+    seenset = set()
+    orbits = []
+    def step_point(point):
+        res = []
+        for gen in gens:
+            # a generator is a list of cycles, so we go one level deeper
+            for cyc in gen:
+                try:
+                    i = cyc.index(point)
+                    res.append(cyc[(i + 1) % len(cyc)])
+                except:
+                    continue
+        return res
+
+    def dfs_find_orbit(gens, point, cur_orbit):
+        frontier = step_point(point)
+        cur_orbit.update(frontier)
+        # iterate until fixed point
+        for new_point in cur_orbit - set(frontier):
+            cur_orbit.update(dfs_find_orbit(gens, new_point, cur_orbit))
+        return cur_orbit
+
+    while vset != seenset:
+        diff = vset - seenset
+        cur_point = diff.pop()
+        new_orbit = dfs_find_orbit(gens, cur_point, set([cur_point]))
+        seenset.update(new_orbit)
+        orbits.append(new_orbit)
+
+    return orbits
+
 ### bfs_foldrep
 ### visits all possible isomorphic graph colorings in `colors`. The colors in
 ### `fixcolors` are fixed throughout.
-def bfs_foldrep(graph, colors, fixcolors, acc, f):
+### if `orders` is True, provide `f` with the orders of the colored and
+### uncolored graph automorphism group
+def bfs_foldrep(graph, colors, fixcolors, acc, f, orders=False):
     # queue of colorings yet to be considered
     queue = deque()
     queue.append([[], colors])
@@ -18,8 +54,6 @@ def bfs_foldrep(graph, colors, fixcolors, acc, f):
         # TODO: turn this into a single GI call by having it return both the
         # automorphism group and the canonical fora
         part = c + fixcolors
-        print("partition: %s" % part)
-        print("graph: %s" % (graph))
         gcanon, cert = my_bliss.canonical_form(graph, partition=part, certificate=True,
                                                return_graph=False)
         gcanon = tuple(gcanon)
@@ -38,40 +72,14 @@ def bfs_foldrep(graph, colors, fixcolors, acc, f):
 
         A, orbits = graph.automorphism_group(partition=c + fixcolors, orbits=True,
                                              algorithm="bliss")
-        acc = f(acc, graph, c, A)
+        acc = f(acc, graph, c)
         reps.add((gcanon, (tuple(c_canon[0]), tuple(c_canon[1]))))
         print(len(reps))
         if len(c_canon[0]) + 1 <= len(colors) / 2.0:
             # if we can add more colors, try to
 
             # gens = my_bliss.raw_automorphism_generators(graph, partition=c + fixcolors)
-            # # big speed hack: compute orbits from generators. avoid calling GAP for this.
-            # vset = set(colors)
-            # seenset = set()
-            # orbits = []
-            # def step_point(point):
-            #     res = []
-            #     for gen in gens:
-            #         # a generator is a list of cycles, so we go one level deeper
-            #         for cyc in gen:
-            #             try:
-            #                 i = cyc.index(point)
-            #                 res.append(cyc[(i + 1) % len(cyc)])
-            #             except:
-            #                 continue
-            #     return res
-            # def dfs_find_orbit(gens, point, cur_orbit):
-            #     frontier = step_point(point)
-            #     cur_orbit.update(frontier)
-            #     for new_point in cur_orbit - set(frontier):
-            #         cur_orbit.update(dfs_find_orbit(gens, new_point, cur_orbit))
-            #     return cur_orbit
-            # while vset != seenset:
-            #     diff = vset - seenset
-            #     cur_point = diff.pop()
-            #     new_orbit = dfs_find_orbit(gens, cur_point, set([cur_point]))
-            #     seenset.update(new_orbit)
-            #     orbits.append(new_orbit)
+            # orbits = generate_orbits(colors, gens)
             # if we need the full automorphism group, we can compute it here.
             # expand this node and add it to the queue
             for o in orbits:
@@ -87,7 +95,7 @@ def bfs_foldrep(graph, colors, fixcolors, acc, f):
 def genrep_bfs(G):
     colors = G.vertices()
     # folding function
-    def add_vertex(acc, g, color, A):
+    def add_vertex(acc, g, color):
         if len(color[0]) < len(g.vertices()) / 2.0:
             return acc + [color] + [color[::-1]]
 
@@ -97,45 +105,13 @@ def add_vertex(acc, g, color, A):
 
 def genrep_bfs_factor(G, varnodes, colornodes):
     # folding function
-    def add_vertex(acc, g, color, A):
+    def add_vertex(acc, g, color):
         if len(color[0]) < len(varnodes) / 2.0:
             return acc + [color] + [color[::-1]]
 
         return acc + [color]
     return bfs_foldrep(G, varnodes, colornodes, [], add_vertex)
 
-### partition_function
-### computes the partition of a fully symmetric MLN the specified parameter
-###
-### potential: a function which maps variable truth assignments to
-### probabilities. This function *must* be invariant wrt. the automorphism
-### group of the graph.
-def partition(G, variables, factors, potential):
-    # folding function
-    print("G: %s, vars: %s, factors: %s" % (G, variables, factors))
-    def sum_prob(acc, g, color, A):
-        # TODO: avoid recomputing these two automorphism groups
-        g_order = g.automorphism_group().order()
-        aut_order = A.order()
-        orbit_sz = g_order / aut_order # yay orbit stabilizer theorem!
-        if len(color[0]) < len(g.vertices()) / 2.0:
-            return acc + (orbit_sz * potential(color)) + (orbit_sz * potential(color[::-1]))
-
-        return acc + (orbit_sz * potential(color))
-    return bfs_foldrep(G, variables, factors, 0.0, sum_prob)
-
-def bruteforce_partition(G, potential):
-    # evaluate the potential on every assignment to variables
-    def recurse(c1, c2):
-        if len(c1) == 0:
-            return 0
-        tot = potential([c1, c2])
-        for x in c1:
-            new_c1 = c1[:]
-            new_c1.remove(x)
-            tot += recurse(new_c1, c2 + [x])
-        return tot
-    return recurse(G.vertices(), [])
 
 # given a graph g, count the number of distinct 2-colorings using Polya's
 # theorem.
@@ -144,22 +120,6 @@ def recurse(c1, c2):
 def count_num_distinct(g):
     return g.automorphism_group().cycle_index().expand(2)((1,1))
 
-
-def compute_partition():
-    G = gen_complete_pairwise_factor(7)
-    def potential(assgn):
-        # count potential: potential is # incoming nodes which are true
-        p = 0.0
-        for factor in G[1][1]:
-            connected = G[0].neighbors(factor)
-            for v in connected:
-                if v in assgn[0]:
-                    p += 1
-        return p
-    part = partition(G[0], G[1][0], [G[1][1]], potential)
-    print("partition: %d" % part)
-    # print("brute forced: %d" % (bruteforce_partition(G, partfun)))
-
 def find_representatives():
     # n=3
     pr = cProfile.Profile()
@@ -168,7 +128,7 @@ def find_representatives():
     # G = graphs.EmptyGraph()
     # for v in range(0, 300):
     #     G.add_vertex(v)
-    G = gen_complete_pairwise_factor(7)
+    G = gen_complete_pairwise_factor(30)
     # G = graphs.CompleteGraph(5)
     # G = gen_pigeonhole(5,5)
     # G = graphs.CycleGraph(20)
@@ -191,5 +151,4 @@ def find_representatives():
 
 
 if __name__ == "__main__":
-    compute_partition()
-    # find_representatives()
+    find_representatives()

query.py

@@ -0,0 +1,116 @@
+from sage.all import *
+from my_graphs import *
+import cProfile, pstats, StringIO
+from collections import deque
+import my_bliss
+from orbitgen import *
+
+### partition_function
+### computes the partition of a fully symmetric MLN the specified parameter
+###
+### potential: a function which maps variable truth assignments to
+### probabilities. This function *must* be invariant wrt. the automorphism
+### group of the graph.
+def partition(G, variables, factors, potential):
+    # folding function
+    print("G: %s, vars: %s, factors: %s" % (G, variables, factors))
+    def sum_prob(acc, g, color, A):
+        # TODO: avoid recomputing these two automorphism groups
+        g_order = g.automorphism_group().order()
+        aut_order = A.order()
+        orbit_sz = g_order / aut_order # yay orbit stabilizer theorem!
+        if len(color[0]) < len(g.vertices()) / 2.0:
+            return acc + (orbit_sz * potential(color)) + (orbit_sz * potential(color[::-1]))
+
+        return acc + (orbit_sz * potential(color))
+    return bfs_foldrep(G, variables, factors, 0.0, sum_prob)
+
+
+def compute_partition():
+    G = gen_complete_pairwise_factor(50)
+    def potential(assgn):
+        # count potential: potential is # incoming nodes which are true
+        p = 0.0
+        for factor in G[1][1]:
+            connected = G[0].neighbors(factor)
+            for v in connected:
+                if v in assgn[0]:
+                    p += 1
+        return p
+
+    pr = cProfile.Profile()
+    pr.enable()
+
+    part = partition(G[0], G[1][0], [G[1][1]], potential)
+    print("partition: %d" % part)
+    # 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 mpe(G, variables, factors, potential):
+    # folding function
+    print("G: %s, vars: %s, factors: %s" % (G, variables, factors))
+    def sum_prob(acc, g, color):
+        (cur_max_state, cur_max_value) = acc
+        if len(color[0]) < len(g.vertices()) / 2.0:
+            # check both possibilities
+            pot = potential(color[::-1])
+            if cur_max_value < pot:
+                cur_max_value = pot
+                cur_max_state = color[::-1]
+        pot = potential(color)
+        if cur_max_value < pot:
+            cur_max_value = pot
+            cur_max_state = color
+        return (cur_max_state, cur_max_value)
+    return bfs_foldrep(G, variables, factors, (None, 0.0), sum_prob)
+
+
+def compute_mpe_complete_2factor():
+    G = gen_complete_pairwise_factor(10)
+    # 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
+        # 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]:
+                p += 0.3
+            else:
+                p += 0.2
+            for v in connected:
+                if v in assgn[0]:
+                    p += 1
+        return 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]))
+    # 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_mpe_complete_2factor()