switch to python 3 + blocks version of MLE implemented

5 files changed, 116 insertions, 12 deletions

diff --git a/simulation/main.py b/simulation/main.py
index decdf8f..1c0b3e8 100644
--- a/simulation/main.py
+++ b/simulation/main.py
@@ -7,6 +7,7 @@ from scipy.optimize import minimize
 import matplotlib.pyplot as plt
 import seaborn
 from random import random, randint
+from six.moves import range
 
 seaborn.set_style("white")
 
@@ -39,7 +40,7 @@ def uniform_source(graph, *args, **kwargs):
 
 
 def simulate_cascades(n, graph, source=uniform_source):
-    for t in xrange(n):
+    for t in range(n):
         x0 = source(graph, t)
         yield simulate_cascade(x0, graph)
 
diff --git a/simulation/mleNode.py b/simulation/mleNode.py
index ed32a12..c6b2e85 100644
--- a/simulation/mleNode.py
+++ b/simulation/mleNode.py
@@ -26,7 +26,7 @@ def test_gradient(x, y):
         f1 = likelihood(p, x, y)
         p[i] -= 2 * eps
         f2 = likelihood(p, x, y)
-        print g[i], (f1 - f2) / (2 * eps)
+        print(g[i], (f1 - f2) / (2 * eps))
 
 
 def infer(x, y):
diff --git a/simulation/mle_blocks.py b/simulation/mle_blocks.py
new file mode 100644
index 0000000..5acebab
--- /dev/null
+++ b/simulation/mle_blocks.py
@@ -0,0 +1,95 @@
+import main as mn
+import theano
+from theano import tensor as tsr
+import blocks
+import blocks.algorithms, blocks.main_loop, blocks.extensions.monitoring
+import theano.tensor.shared_randomstreams
+import picklable_itertools
+import numpy as np
+from six.moves import range
+import fuel
+import fuel.datasets
+import collections
+
+
+class JeaninuScheme(fuel.schemes.ShuffledScheme):
+    def get_request_iterator(self):
+        indices = list(self.indices)
+        start = np.random.randint(self.batch_size)
+        batches = list(map(
+            list,
+            picklable_itertools.extras.partition_all(self.batch_size,
+                                                     indices[start:])
+        ))
+        if indices[:start]:
+            batches.append(indices[:start])
+        batches = np.asarray(batches)
+        return iter(batches[np.random.permutation(len(batches))])
+
+
+def create_model(n_nodes):
+    x = tsr.matrix(name='x', dtype='int8')
+    s = tsr.matrix(name='s', dtype='int8')
+    params = theano.shared(
+        .5 + .01 *
+        np.random.normal(size=(n_nodes, n_nodes)).astype(theano.config.floatX),
+        name='params'
+    )
+    y = tsr.maximum(tsr.dot(x, params), 1e-5)
+    infect = tsr.log(1. - tsr.exp(-y[0:-1]))
+    lkl_pos = tsr.sum(infect * (x[1:] & s[1:]))
+    lkl_neg = tsr.sum(-y[0:-1] * (~x[1:] & s[1:]))
+    lkl_mle = lkl_pos + lkl_neg
+    lkl_mle.name = 'cost'
+    return x, s, params, lkl_mle
+
+
+def create_random_graph(n_nodes, p=.5):
+    graph = .5 * np.random.binomial(2, p=.5, size=(n_nodes, n_nodes))
+    for k in range(len(graph)):
+        graph[k, k] = 0
+    return np.log(1. / (1 - p * graph))
+
+
+def create_data_stream(n_cascades, graph, batch_size, shuffle=True):
+    """
+    shuffle (bool): shuffle minibatches but not within minibatch
+    """
+    cascades = mn.build_cascade_list(mn.simulate_cascades(n_cascades, graph),
+                                     collapse=True)
+    x_obs, s_obs = cascades[0], cascades[1]
+    data_set = fuel.datasets.base.IndexableDataset(collections.OrderedDict(
+        [('x', x_obs), ('s', s_obs)]
+    ))
+    if shuffle:
+        scheme = JeaninuScheme(len(x_obs), batch_size=batch_size)
+    else:
+        scheme = fuel.schemes.SequentialScheme(len(x_obs),
+                                               batch_size=batch_size)
+    return fuel.streams.DataStream(dataset=data_set, iteration_scheme=scheme)
+
+
+if __name__ == "__main__":
+    n_cascades = 10000
+    batch_size = 1000
+    graph = np.array([[0, 0, 1], [0, 0, 0.5], [0, 0, 0]])
+    graph = np.log(1. / (1 - .5 * graph))
+    print('GRAPH:\n', graph, '\n-------------\n')
+
+    x, s, params, cost = create_model(len(graph))
+
+    alg = blocks.algorithms.GradientDescent(
+       cost=-cost, parameters=[params], step_rule=blocks.algorithms.AdaDelta()
+    )
+    data_stream = create_data_stream(n_cascades, graph, batch_size,
+            shuffle=True)
+    loop = blocks.main_loop.MainLoop(
+        alg, data_stream,
+        extensions=[
+            blocks.extensions.FinishAfter(after_n_epochs = 1000),
+            blocks.extensions.monitoring.TrainingDataMonitoring([cost, params],
+                after_batch=True),
+            blocks.extensions.Printing()
+        ]
+    )
+    loop.run()
diff --git a/simulation/vi.py b/simulation/vi.py
index 1e45761..b4a83c1 100644
--- a/simulation/vi.py
+++ b/simulation/vi.py
@@ -73,7 +73,7 @@ if __name__ == '__main__':
     p = 0.5
     graph = np.log(1. / (1 - p * graph))
     print(graph)
-    cascades = mn.build_cascade_list(mn.simulate_cascades(100, graph))
+    cascades = mn.build_cascade_list(mn.simulate_cascades(1000, graph))
     mu0, sig0 = (1. + .2 * np.random.normal(size=graph.shape),
                  1 + .2 * np.random.normal(size=graph.shape))
     mu1, sig1 = (1. + .2 * np.random.normal(size=graph.shape),
diff --git a/simulation/vi_theano.py b/simulation/vi_theano.py
index 9e8fdb0..cfd1dcf 100644
--- a/simulation/vi_theano.py
+++ b/simulation/vi_theano.py
@@ -5,11 +5,11 @@ import theano.tensor.shared_randomstreams
 import numpy as np
 
 n_cascades = 1000
-n_nodes = 4
+n_nodes = 3
 n_samples = 100
 srng = tsr.shared_randomstreams.RandomStreams(seed=123)
-lr = 5*1e-3
-n_epochs = 20
+lr = 1e-1
+n_epochs = 1
 
 ###############Variational Inference####################
 
@@ -70,28 +70,36 @@ if __name__ == "__main__":
         graph[k, k] = 0
     p = 0.5
     graph = np.log(1. / (1 - p * graph))
+
+    graph = np.array([[0, 0, 1], [0, 0, 0.5], [0, 0, 0]])
+    p = 0.5
+    graph = np.log(1. / (1 - p * graph))
+
     cascades = mn.build_cascade_list(mn.simulate_cascades(n_cascades, graph),
                                      collapse=True)
     x_obs, s_obs = cascades[0], cascades[1]
 
     #mle
     lkl_plot = []
-    if 0:
+    if 1:
         for i in range(n_epochs):
-            for xt, st in zip(x_obs, s_obs):
+            for k in xrange(len(x_obs)/100):
+                xt = x_obs[k*100:(k+1)*100]
+                st = s_obs[k*100:(k+1)*100]
                 lkl = train_mle(xt, st)
                 lkl_plot.append(lkl)
+                print(params.get_value())
         print(graph)
         w = params.get_value()
         for k in range(len(w)):
             w[k, k] = 0
         print(w)
-        import matplotlib.pyplot as plt
-        plt.plot(lkl_plot)
-        plt.show()
+        #import matplotlib.pyplot as plt
+        #plt.plot(lkl_plot)
+        #plt.show()
 
     #variational inference
-    if 1:
+    if 0:
         for i in range(n_epochs):
             train_kl()
             for k in xrange(len(x_obs)/100):