From c823f74aba5a2fb8970515afef6121e5f7143a80 Mon Sep 17 00:00:00 2001
From: Steven Holtzen <sholtzen@gmail.com>
Date: Fri, 8 Mar 2019 09:51:40 -0800
Subject: [PATCH] orbit gen exp

---
 factor.py | 104 +++++++++++---
 res.txt   | 407 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 491 insertions(+), 20 deletions(-)
 create mode 100644 res.txt

diff --git a/factor.py b/factor.py
index b2fe231..24c21f9 100644
--- a/factor.py
+++ b/factor.py
@@ -15,6 +15,12 @@ def findsubsets(S,m):
 def flip():
     return numpy.random.binomial(1, 0.5)
 
+def div_or_0(a, b):
+    if b == 0.0:
+        return 0
+    else:
+        return a / b
+
 ### returns a random element from a group g using product replacement
 def fast_random_element(g):
     stab = libgap(g)
@@ -231,10 +237,26 @@ def gibbs_transition(self):
                 st2.sort()
                 st2 = tuple(st2)
                 st2idx = state_to_idx[st2]
-                transition[idx, st1idx] += (1.0 / len(self.variables)) * (float(st1prob) / (st1prob + st2prob))
-                transition[idx, st2idx] += (1.0 / len(self.variables)) * (float(st2prob) / (st1prob + st2prob))
+                if st1prob + st2prob > 0:
+                    transition[idx, st1idx] += ((1.0/ len(self.variables)) * float(st1prob)/ (st1prob + st2prob))
+                    transition[idx, st2idx] += ((1.0 / len(self.variables)) * float(st2prob)/ (st1prob + st2prob))
+                else:
+                    # stay
+                    s_idx = state_to_idx[s]
+                    transition[s_idx, st1idx] += 1.0 / len(self.variables)
+        return transition
+
+    def uniform_transition(self):
+        states = self.gen_all_states()
+
+        transition = np.zeros([len(states),len(states)])
+        for x in range(0, len(states)):
+            for y in range(0, len(states)):
+                transition[x, y] = 1.0 / len(states)
+
         return transition
 
+
     def orbit_transition(self):
         states = self.gen_all_states()
         state_to_idx = dict()
@@ -272,14 +294,29 @@ def burnside_mh_transition(self, burnsidesteps):
             state_to_idx[st] = idx
             idx_to_state[idx] = st
 
+        aut_cache = dict()
         B = np.linalg.matrix_power(self.burnside_transition(), burnsidesteps)
         for x in range(0, len(states)):
             for y in range(0, len(states)):
                 st_x = idx_to_state[x]
-                orbx = self.graph_aut_order / self.graph.automorphism_group(partition=self.state_to_partition(dict(st_x))).order()
+                part_x = self.state_to_partition(dict(st_x))
+                orbx = None
+                if st_x in aut_cache:
+                    orbx = aut_cache[st_x]
+                else:
+                    orbx = self.graph_aut_order / self.graph.automorphism_group(partition=self.state_to_partition(dict(st_x))).order()
+                    aut_cache[st_x] = orbx
+
                 prx = self.potential(dict(st_x))
                 st_y = idx_to_state[y]
-                orby = self.graph_aut_order / self.graph.automorphism_group(partition=self.state_to_partition(dict(st_y))).order()
+                orby = None
+                if st_y in aut_cache:
+                    orby = aut_cache[st_y]
+                else:
+                    orby = self.graph_aut_order / self.graph.automorphism_group(partition=self.state_to_partition(dict(st_y))).order()
+                    aut_cache[st_y] = orby
+
+                # orby = self.graph_aut_order / self.graph.automorphism_group(partition=self.state_to_partition(dict(st_y))).order()
                 pry = self.potential(dict(st_y))
                 if prx != 0.0:
                     B[x,y] = B[x,y] * np.minimum(1, (pry * orby) / (prx * orbx))
@@ -308,10 +345,16 @@ def brute_force_prob_vector(self):
 
         return vec
 
-    def total_variation(self, transition_matrix, starting_point, num_steps):
+    def total_variation(self, transition_matrix, starting_point, num_steps, stagger=None, staggerlen=0):
         v = self.brute_force_prob_vector()
+        T = transition_matrix
         for i in range(1, num_steps):
-            p = starting_point.dot(np.linalg.matrix_power(transition_matrix, i))
+            if staggerlen > 0 and i % staggerlen == 0:
+                T = np.matmul(T,stagger)
+            else:
+                T = np.matmul(T, transition_matrix)
+
+            p = starting_point.dot(T)
             d = 0.0
             for idx in range(0, len(v)):
                 d += abs(v[idx] - p[idx])
@@ -648,24 +691,29 @@ def query(state):
 
 # n pigeons, m holes
 def mk_pigeonhole_fg(n, m):
-    w1 = 1000
-    w2 = 1000
+    w1 = 10000000
+    w2 = 100000
     g = gen_pigeonhole(n, m)
     def potential(state):
         total = 0.0
-        # to see every pigeon in some hole
+        # to see every pigeon in exactly one hole
         for p in range(0, n):
             # check the holes for the pigeons
             in_hole = False
             for h in range(0, m):
                 if state[(p, h)]:
-                    in_hole = True
+                    if in_hole:
+                        return 0.0
+                    else:
+                        in_hole = True
             if in_hole:
                 total += w1
-        # check to see no pigeon is in more than 1 hole
+
+
+        # check to see no no hole has 2 pigeons
         for h in range(0, m):
             for (p1, p2) in findsubsets(range(0, n), 2):
-                if state[(p1, h)] and state[(p2, h)]:
+                if not state[(p1, h)] or not state[(p2, h)]:
                     total += w2
 
         return total
@@ -693,24 +741,40 @@ def potential(state):
 
 
 def total_var():
-    # model = mk_pigeonhole_fg(2, 4)
-    model = mk_simple_complete_fg(6)
+    model = mk_pigeonhole_fg(2, 5)
+    # model = mk_simple_complete_fg(6)
     # print(np.linalg.matrix_power(graph.burnside_transition(), 10))
     # print("------")
     # print(graph.burnside_mh_transition(4))
     # M = graph.burnside_mh_transition(4)
     # M = graph.burnside_mh_transition(4)
 
-    M = model.gibbs_transition()
-    M2 = model.orbit_transition()
-    M = np.matmul(M2, M)
-    # M = model.burnside_mh_transition(4)
+    gibbs = model.gibbs_transition()
+    # print(gibbs)
+    # print(sum(gibbs))
+    within_orbit = model.orbit_transition()
+    orbitalmcmc = np.matmul(within_orbit, gibbs)
+    # unif = model.uniform_transition()
+    burnside = model.burnside_mh_transition(4)
+    M = orbitalmcmc
 
     pv = model.brute_force_prob_vector()
     start = np.zeros([2**len(model.variables)])
-    start[0] = 1
+    start[10] = 1
     # print(np.linalg.matrix_power(M, 5))
-    print(model.total_variation(M, start, 100))
+    print("-------------------")
+    print("pure gibbs")
+    model.total_variation(gibbs, start, 100)
+    print("------------------")
+    print("pure jump")
+    model.total_variation(burnside, start, 100, stagger=burnside, staggerlen=200)
+    print("------------------")
+    print("5-stagger")
+    print(model.total_variation(gibbs, start, 100, stagger=burnside, staggerlen=5))
+    print("------------------")
+    print("10-stagger")
+    print(model.total_variation(gibbs, start, 100, stagger=burnside, staggerlen=10))
+
 
 def pigeonhole_exact_experiment():
     fg = mk_pigeonhole_fg(2,34)
diff --git a/res.txt b/res.txt
new file mode 100644
index 0000000..7d1e5ce
--- /dev/null
+++ b/res.txt
@@ -0,0 +1,407 @@
+[(0, 0), (0, 1), (0, 2), (0, 3), (0, 4), (1, 0), (1, 1), (1, 2), (1, 3), (1, 4)]
+-------------------
+pure gibbs
+1	1.5981501538664418
+2	1.5004375195247162
+3	1.4496907228487779
+4	1.4182910989845776
+5	1.3981186925395104
+6	1.3849631390089514
+7	1.377777799999998
+8	1.3777777800000028
+9	1.3777777780000016
+10	1.377777777799998
+11	1.3777777777799998
+12	1.3777777777779998
+13	1.3777777777777995
+14	1.3777777777777809
+15	1.3777777777777789
+16	1.3777777777777798
+17	1.3777777777777802
+18	1.3777777777777769
+19	1.3777777777777762
+20	1.3777777777777773
+21	1.3777777777777758
+22	1.3777777777777795
+23	1.3777777777777764
+24	1.3777777777777769
+25	1.3777777777777758
+26	1.3777777777777755
+27	1.3777777777777773
+28	1.377777777777776
+29	1.3777777777777782
+30	1.377777777777775
+31	1.3777777777777787
+32	1.3777777777777782
+33	1.3777777777777778
+34	1.3777777777777744
+35	1.377777777777775
+36	1.3777777777777764
+37	1.3777777777777747
+38	1.377777777777776
+39	1.3777777777777753
+40	1.3777777777777764
+41	1.3777777777777755
+42	1.377777777777775
+43	1.3777777777777775
+44	1.3777777777777747
+45	1.377777777777777
+46	1.3777777777777778
+47	1.3777777777777747
+48	1.3777777777777733
+49	1.3777777777777735
+50	1.377777777777775
+51	1.377777777777776
+52	1.377777777777771
+53	1.3777777777777729
+54	1.3777777777777764
+55	1.3777777777777733
+56	1.3777777777777727
+57	1.3777777777777742
+58	1.3777777777777738
+59	1.377777777777775
+60	1.377777777777776
+61	1.3777777777777747
+62	1.3777777777777742
+63	1.3777777777777749
+64	1.377777777777772
+65	1.3777777777777733
+66	1.3777777777777749
+67	1.3777777777777722
+68	1.3777777777777713
+69	1.3777777777777716
+70	1.3777777777777753
+71	1.3777777777777707
+72	1.3777777777777733
+73	1.3777777777777702
+74	1.3777777777777707
+75	1.3777777777777702
+76	1.3777777777777727
+77	1.3777777777777733
+78	1.3777777777777707
+79	1.3777777777777707
+80	1.377777777777771
+81	1.3777777777777716
+82	1.3777777777777713
+83	1.377777777777771
+84	1.3777777777777716
+85	1.3777777777777693
+86	1.3777777777777724
+87	1.3777777777777698
+88	1.3777777777777702
+89	1.3777777777777702
+90	1.3777777777777707
+91	1.3777777777777698
+92	1.377777777777771
+93	1.3777777777777698
+94	1.377777777777771
+95	1.3777777777777704
+96	1.3777777777777698
+97	1.3777777777777698
+98	1.3777777777777684
+99	1.3777777777777707
+------------------
+pure jump
+1	1.6635003723115502
+2	1.5166433359182248
+3	1.3825953516110996
+4	1.2602591788325461
+5	1.1486284933971642
+6	1.0467813252904412
+7	0.9538732652636945
+8	0.883764596250683
+9	0.8444326655038324
+10	0.8068360858839644
+11	0.7709011737698167
+12	0.7365569337568991
+13	0.7037350244697836
+14	0.6723697144309758
+15	0.6423978299006284
+16	0.613758696332558
+17	0.586394074857604
+18	0.5602480950009882
+19	0.5352671846627567
+20	0.5113999982358108
+21	0.48859734360152307
+22	0.46681210862697753
+23	0.44599918768640534
+24	0.42611540864241554
+25	0.4071194606469051
+26	0.388971823056626
+27	0.3716346957022396
+28	0.3550719307012587
+29	0.3392489659645463
+30	0.3241327605100748
+31	0.3096917316681709
+32	0.2958956942361696
+33	0.28271580161945054
+34	0.27012448897707475
+35	0.25809541837515876
+36	0.2466034259383597
+37	0.23562447097922168
+38	0.2251355870764453
+39	0.21511483506641446
+40	0.20554125790627836
+41	0.19639483736286117
+42	0.18765645247807033
+43	0.1793078397590062
+44	0.17133155503937333
+45	0.16371093695743477
+46	0.156430071995329
+47	0.1494737610243802
+48	0.1428274873011619
+49	0.13647738585976046
+50	0.13041021424628954
+51	0.12461332454282445
+52	0.1190746366289855
+53	0.11378261263069665
+54	0.10872623250699996
+55	0.10389497072729391
+56	0.09927877399274633
+57	0.09486803995731417
+58	0.09065359690516432
+59	0.08662668434298745
+60	0.082778934467164
+61	0.07910235446723689
+62	0.07558930962875184
+63	0.07223250719987348
+64	0.06902498098774269
+65	0.0659600766518591
+66	0.06303143766320209
+67	0.06023299189910032
+68	0.05755893884514074
+69	0.05500373737668573
+70	0.052562094093695035
+71	0.05022895218380866
+72	0.047999480789651446
+73	0.045869064857448225
+74	0.04383329544505926
+75	0.041887960468506444
+76	0.04002903586703544
+77	0.03825267716762546
+78	0.036555211430775995
+79	0.03493312956015295
+80	0.03338307895956135
+81	0.03190185652137959
+82	0.03048640193137863
+83	0.029133791275510492
+84	0.027841230934919785
+85	0.026606051756052396
+86	0.025425703483365293
+87	0.02429774944267017
+88	0.023219861463748784
+89	0.022189815031333902
+90	0.02120548465412679
+91	0.02026483944192541
+92	0.019365938881447088
+93	0.018506928801818073
+94	0.01768603752114526
+95	0.01690157216599194
+96	0.016151915155893613
+97	0.01543552084548374
+98	0.014750912317087809
+99	0.014096678317008365
+------------------
+5-stagger
+1	1.5981501538664418
+2	1.5004375195247162
+3	1.4496907228487779
+4	1.4182910989845776
+5	1.2976802849144122
+6	1.2252678114941664
+7	1.203277682519918
+8	1.1926152826590317
+9	1.1875231764482665
+10	1.084517289429584
+11	1.033354184418927
+12	1.0227497535327723
+13	1.0222585334638412
+14	1.0222094114569509
+15	0.9316458902375643
+16	0.8856564490286666
+17	0.8803638712259338
+18	0.8799410525277808
+19	0.8798987706579686
+20	0.8019432900106778
+21	0.7620523908135743
+22	0.7578007010379936
+23	0.7574367467550058
+24	0.7574003513267072
+25	0.6902986213160989
+26	0.6559611611459036
+27	0.6523006234481656
+28	0.6519873384318713
+29	0.6519560099302396
+30	0.5941968409295196
+31	0.5646396685123674
+32	0.5614881363466384
+33	0.5612184664833659
+34	0.5611914994970408
+35	0.5114740291415684
+36	0.48603167153918225
+37	0.4833184515315103
+38	0.4830863247735172
+39	0.483063112097719
+40	0.44026762925490565
+41	0.4183672550833013
+42	0.4160314536843645
+43	0.41583164332690337
+44	0.41581166229115646
+45	0.3789743846578366
+46	0.3601229038038671
+47	0.35811206856736694
+48	0.35794007556107066
+49	0.35792287626044283
+50	0.3262142306193522
+51	0.3099871944590654
+52	0.30825614870671225
+53	0.3081081003543189
+54	0.30809329551907766
+55	0.28079922397414786
+56	0.26683127085247876
+57	0.2653411084290785
+58	0.26521367118543904
+59	0.265200927461077
+60	0.24170679878104068
+61	0.22968342923068424
+62	0.2284006470520044
+63	0.2282909514640227
+64	0.22827998190522183
+65	0.20805674466849794
+66	0.19770724131879178
+67	0.1966029910805106
+68	0.19650856717679213
+69	0.19649912478642004
+70	0.1790913896694809
+71	0.1701827225802727
+72	0.1692321654982083
+73	0.1691508871737393
+74	0.16914275934129158
+75	0.154158541569524
+76	0.14649012187469573
+77	0.14567187244615817
+78	0.14560190958962954
+79	0.14559491330397778
+80	0.13269680384886698
+81	0.12609596766208842
+82	0.12539161429210824
+83	0.12533139157805123
+84	0.12532536930664562
+85	0.1142229357153165
+86	0.10854105666956562
+87	0.10793474862892571
+88	0.10788291004673449
+89	0.10787772618851539
+90	0.09832097171222598
+91	0.09343011464219139
+92	0.09290820624137146
+93	0.09286358456256318
+94	0.09285912239468265
+95	0.08463285597825658
+96	0.08042289642725103
+97	0.0799736405248435
+98	0.07993523102104463
+99	0.07993139007066502
+None
+------------------
+10-stagger
+1	1.5981501538664418
+2	1.5004375195247162
+3	1.4496907228487779
+4	1.4182910989845776
+5	1.3981186925395104
+6	1.3849631390089514
+7	1.377777799999998
+8	1.3777777800000028
+9	1.3777777780000016
+10	1.2558275791818083
+11	1.1936392370877034
+12	1.187046322469131
+13	1.1864755310148118
+14	1.186418451869384
+15	1.1864127439548382
+16	1.1864121731633877
+17	1.186412116084244
+18	1.1864121103763288
+19	1.1864121098055278
+20	1.0805432175357892
+21	1.0260910695080114
+22	1.0203998900894466
+23	1.0199118043442146
+24	1.0198629957696865
+25	1.0198581149122372
+26	1.0198576268264945
+27	1.0198575780179182
+28	1.0198575731370632
+29	1.0198575726489745
+30	0.9288473584406167
+31	0.8820342216852579
+32	0.8771411373813263
+33	0.8767215926664986
+34	0.8766796381950113
+35	0.8766754427478628
+36	0.8766750232031527
+37	0.8766749812486786
+38	0.8766749770532347
+39	0.8766749766336897
+40	0.7984425162873132
+41	0.7582016107627431
+42	0.7539951149359403
+43	0.7536344722050614
+44	0.7535984079319726
+45	0.7535948015046651
+46	0.7535944408619337
+47	0.753594404797659
+48	0.7535944011912354
+49	0.7535944008305923
+50	0.6863455815826366
+51	0.6517542606754316
+52	0.6481381489387478
+53	0.6478281385594645
+54	0.6477971375215369
+55	0.647794037417747
+56	0.6477937274073667
+57	0.647793696406327
+58	0.6477936933062212
+59	0.6477936929962174
+60	0.5899863620776037
+61	0.5602514677532395
+62	0.557142946668039
+63	0.5568764601227036
+64	0.5568498114681747
+65	0.5568471466027209
+66	0.5568468801161734
+67	0.5568468534675185
+68	0.5568468508026574
+69	0.5568468505361686
+70	0.5071554190516555
+71	0.4815951413518653
+72	0.47892299431872126
+73	0.4786939210889535
+74	0.47867101376597354
+75	0.4786687230336783
+76	0.47866849396044814
+77	0.4786684710531288
+78	0.47866846876239244
+79	0.47866846853332023
+80	0.43595348013957225
+81	0.4139817281850074
+82	0.41168471369677767
+83	0.4114878011489741
+84	0.4114681098941938
+85	0.41146614076871757
+86	0.4114659438561663
+87	0.41146592416491073
+88	0.4114659221957896
+89	0.411465921998875
+90	0.3747479073996548
+91	0.35586087267570043
+92	0.3538863355949171
+93	0.35371706853773643
+94	0.35370014183201626
+95	0.3536984491614457
+96	0.353698279894387
+97	0.3536982629676809
+98	0.35369826127501086
+99	0.3536982611057453
+None