diff --git a/orbitgen.py b/orbitgen.py index 9c59cf6..48a7d41 100644 --- a/orbitgen.py +++ b/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() diff --git a/query.py b/query.py new file mode 100644 index 0000000..dfb52cb --- /dev/null +++ b/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()