First alpha wighted sampling

examples/ParticleFilterDBN.py

@@ -30,24 +30,30 @@ twoTBN.add_edge(weather, weather, True);
 
 weather0_cpt = numpy.array([.4, .6])
 weather0.set_probability_table(weather0_cpt, [weather0])
-ice_cream_eaten0_cpt = numpy.array([[.9, .1],
-                                    [.2, .8]])
+ice_cream_eaten0_cpt = numpy.array([[.2, .8],
+                                    [.9, .1]])
 ice_cream_eaten0.set_probability_table(ice_cream_eaten0_cpt, [weather0, ice_cream_eaten0])
 
-weather_cpt=numpy.array([[.7, .3],
-                         [.5, .5]])
+weather_cpt=numpy.array([[.5, .5],
+                         [.7, .3]])
 weather.set_probability_table(weather_cpt, [weather, weather])
-ice_cream_eaten_cpt = numpy.array([[.9, .1],
-                                   [.2, .8]])
+ice_cream_eaten_cpt = numpy.array([[.2, .8],
+                                   [.9, .1]])
 ice_cream_eaten.set_probability_table(ice_cream_eaten_cpt, [weather, ice_cream_eaten])
 
 dbn.set_B0(B0)
 dbn.set_TwoTBN(twoTBN)
 
-N = 20000
+#N = 20000
+#T = 2
+#evidence = [{weather:"Rain"}, {weather:"Sun"}]
+#samples = fs.sample_DBN(dbn, N, T, evidence)
+#print("P(W_1 = R, W_2 = S) = " + str(len(samples) * 1.0 / N))
+
+N = 1
 T = 2
-evidence = [{weather:"Rain"}, {weather:"Sun"}]
-samples = fs.sample_DBN(dbn, N, T, evidence)
+evidence = [{weather: "Rain"}, {ice_cream_eaten: True}]
+samples = fs.weighted_sample_DBN(dbn, N, T, evidence)
 print("P(W_1 = R, W_2 = S) = " + str(len(samples) * 1.0 / N))
 
 #pt = ProbabilityTable()

primo/core/TwoTBN.py

@@ -6,10 +6,10 @@ from primo.core import BayesNet
 class TwoTBN(BayesNet):
     ''' This is the implementation of a 2-time-slice Bayesian network (2-TBN).
     '''
-    
+
     def __init__(self):
-        BayesNet.__init__(self)            
-        
+        BayesNet.__init__(self)
+
 
     def add_edge(self, node_from, node_to, inter=False):
         '''Add an directed edge to the graph.
@@ -17,7 +17,7 @@ class TwoTBN(BayesNet):
         Keyword arguments:
         node_from -- from node
         node_to -- to node
-        inter -- is this edge a temporal (inter-slice) conditional dependency 
+        inter -- is this edge a temporal (inter-slice) conditional dependency
         (default: False)
         '''
         super(TwoTBN, self).add_edge(node_from, node_to)
@@ -51,4 +51,17 @@ class TwoTBN(BayesNet):
             if successor == origin:
                 return True
             else:
-                return self.has_loop(successor, origin)
+                return self.has_loop(successor, origin)
\ No newline at end of file
+
+    def get_nodes_in_topological_sort(self):
+        # Remove temporal/inter-slice edges for toplogical sort
+        inter_edges = []
+        for edge in self.graph.edges(data = True):
+            if edge[2]['inter'] == True:
+                inter_edges.append(edge)
+
+        self.graph.remove_edges_from(inter_edges)
+        ts = super(TwoTBN, self).get_nodes_in_topological_sort()
+        # Add temporal/inter-slice edges
+        self.graph.add_edges_from(inter_edges)
+        return ts
\ No newline at end of file

primo/reasoning/particlebased/ForwardSampling.py

 # -*- coding: utf-8 -*-
 from primo.core import BayesNet
 from primo.core import DynamicBayesNet
-from primo.core import TwoTBN
 import random
 import copy
 
@@ -16,11 +15,7 @@ def sample(network, state={}):
     if not isinstance(network, BayesNet):
         raise Exception("The given network is not an instance of BayesNet.")
 
-    if not isinstance(network, TwoTBN):
-        nodes = network.get_nodes_in_topological_sort()
-    else:
-        nodes = network.get_nodes()
-    #print("#### Sampling ####")
+    nodes = network.get_nodes_in_topological_sort()
     for node in nodes:
         # reduce this node's cpd
         parents = network.get_parents(node)
@@ -32,10 +27,38 @@ def sample(network, state={}):
 
         new_state = weighted_random(reduced_cpd.get_table())
         state[node] = node.get_value_range()[new_state]
-        #print("Node: " + node.name + " (" + str(state[node]) + ")")
-    #print("#### Sampling done ####")
     return state
 
+def weighted_sample(network, events={}, w=1.0, state={}):
+    if not isinstance(network, BayesNet):
+        raise Exception("The given network is not an instance of BayesNet.")
+
+    nodes = network.get_nodes_in_topological_sort()
+    print "---------" + str(w)
+    for node in nodes:
+        print node
+        # reduce this node's cpd
+        parents = network.get_parents(node)
+        if parents:
+            evidence = [(parent, state[parent]) for parent in parents]
+            reduced_cpd = node.get_cpd_reduced(evidence)
+            print reduced_cpd
+        else:
+            reduced_cpd = node.get_cpd()
+
+        if node in events:
+            print("w before: " + str(w))
+            print "foreced state: " + str(events[node])
+            w *= reduced_cpd.get_table()[node.get_value_range().index(events[node])]
+            state[node] = events[node];
+            print("w after: " + str(w))
+        else:
+            new_state = weighted_random(reduced_cpd.get_table())
+            state[node] = node.get_value_range()[new_state]
+        print state[node]
+    return (state, w)
+
+
 def sample_DBN(network, N, T, evidence=[]):
     if not isinstance(network, DynamicBayesNet):
         raise Exception("The given network is not an instance of DynamicBayesNet.")
@@ -44,40 +67,82 @@ def sample_DBN(network, N, T, evidence=[]):
         raise Exception("The list of evidences must have as many entries as " +
         "the number of time slices (T).")
 
-    samples = {}
+    samples = []
     for n in xrange(N):
-        #print("###### Sample " + str(n) + " ######")
         # Sample from initial distribution
         init_state = sample(network.get_B0())
 
-        # Copy nodes from 2TBN in state but keep values
+        # Take nodes from 2TBN in state but keep values
         twoTBN = network.get_TwoTBN()
         state = {}
-        #print("### Initial state ###")
         for node in init_state:
-            #print(str(node.name) + ": " + str(init_state[node]))
             state[twoTBN.get_node(node.name)] = init_state[node]
 
         # Sample N Particles
         tmp_sample = []
         tmp_sample.append(copy.copy(state))
         for t in xrange(T - 1):
-            #yield state
             state = sample(twoTBN, state)
             tmp_sample.append(copy.copy(state))
 
-        #for i in xrange(len(tmp_sample)):
-        #    state = tmp_sample[i]
-        #    for node in state:
-        #        print("Node: " + node.name + " (" + str(state[node]) + ")")
+        # Reject samples not consistent with evidence
+        accept_sample = True
+        if len(evidence) != 0:
+            for t in xrange(T):
+                for node in evidence[t]:
+                    if not evidence[t][node] == tmp_sample[t][node]:
+                        accept_sample = False
+
+        if accept_sample:
+             samples.append(tmp_sample)
+
+    return samples
+
+def weighted_sample_DBN(network, N, T, evidence=[]):
+    if not isinstance(network, DynamicBayesNet):
+        raise Exception("The given network is not an instance of DynamicBayesNet.")
+
+    if len(evidence) != 0 and not len(evidence) == T:
+        raise Exception("The list of evidences must have as many entries as " +
+        "the number of time slices (T).")
+
+    samples = []
+    w_match = 0.0
+    w_all = 0.0
+    for n in xrange(N):
+        ###
+        ### TODO: Solve problem withe evindence for initial and other nodes!
+        ###
+
+
+        # Sample from initial distribution
+        (init_state, w) = weighted_sample(network.get_B0(), evidence[0])
+
+        # Take nodes from 2TBN in state but keep values
+        twoTBN = network.get_TwoTBN()
+        state = {}
+        for node in init_state:
+            state[twoTBN.get_node(node.name)] = init_state[node]
+
+        # Sample N Particles
+        tmp_sample = []
+        tmp_sample.append(copy.copy(state))
+        for t in xrange(T - 1):
+            (state, w) = weighted_sample(twoTBN, evidence[t + 1], w, state)
+            tmp_sample.append(copy.copy(state))
 
+        # Reject samples not consistent with evidence
         accept_sample = True
-        for t in xrange(T):
-            for node in evidence[t]:
-                if not evidence[t][node] == tmp_sample[t][node]:
-                    accept_sample = False
+        if len(evidence) != 0:
+            for t in xrange(T):
+                for node in evidence[t]:
+                    if not evidence[t][node] == tmp_sample[t][node]:
+                        accept_sample = False
 
         if accept_sample:
-             samples[n] = tmp_sample
+             samples.append(tmp_sample)
+             w_match += w
+        w_all += w
 
-    return samples
+    print(str(w_match/w_all))
\ No newline at end of file
+    return samples