Implemented reduction for Discrete Node, started implementation of mcmc

examples/MCMC.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+from  primo.core import BayesNet
+from  primo.reasoning import DiscreteNode
+from  primo.reasoning import MarkovChainSampler
+from primo.reasoning import GibbsTransitionModel
+import numpy
+
+#Construct some simple BayesianNetwork
+bn = BayesNet()
+burglary = DiscreteNode("Burglary", ["Intruder","Safe"])
+alarm = DiscreteNode("Alarm", ["Ringing", "Silent","Kaputt"])
+
+bn.add_node(burglary)
+bn.add_node(alarm)
+
+bn.add_edge(burglary,alarm)
+
+burglary_cpt=numpy.array([0.2,0.8])
+burglary.set_probability_table(burglary_cpt, [burglary])
+
+alarm_cpt=numpy.array([[0.8,0.15,0.05],[0.05,0.9,0.05]])
+alarm.set_probability_table(alarm_cpt, [burglary,alarm])
+
+#Construct a Markov Chain by sampling states from this Network
+
+transition_model = GibbsTransitionModel()
+
+mcs = MarkovChainSampler()
+initial_state={burglary:"Safe",alarm:"Silent"}
+chain = mcs.generateMarkovChain(bn, 20, transition_model, initial_state)
+
+for c in chain:
+    print c

primo/reasoning/DiscreteNode.py

@@ -25,6 +25,12 @@ class DiscreteNode(RandomNode):
 
     def set_probability_table(self, table, nodes):
         self.cpd.set_probability_table(table, nodes)
+        
+    def get_cpd_reduced(self, evidence):
+        return self.cpd.reduction(evidence)
+        
+    def get_cpd(self):
+        return self.cpd
 
     def is_valid(self):
         return self.cpd.is_normalized_as_cpt(self)

primo/reasoning/MCMC.py

+
+
+class GibbsTransitionModel(object):
+    def __init__(self):
+        pass
+        
+    def transition(self, network, state):
+        nodes = network.get_nodes([])
+        for node in nodes:
+            #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)
+            else:
+                reduced_cpd = node.get_cpd()
+                
+            
+                
+            #reduce the children's cpds
+            children = network.get_children(node)
+            for child in children:
+                
+                #reduce this node's cpd
+                parents=network.get_parents(child)
+                if parents:
+                    evidence=[(parent,state[parent]) for parent in parents]
+                    reduced_child_cpd = child.get_cpd_reduced(evidence)
+                else:
+                    reduced_child_cpd = child.get_cpd()
+                reduced_cpd = reduced_cpd.multiplication(reduced_child_cpd)
+                
+            
+        return state
+            
+        
+
+class MarkovChainSampler(object):
+    def __init__(self):
+        pass
+        
+    def generateMarkovChain(self, network, time_steps, transition_model, initial_state):
+        state=initial_state
+        for t in range(time_steps):
+            yield state
+            state=transition_model.transition(network, state)
+        

primo/reasoning/__init__.py

 from RandomNode import RandomNode
 from DiscreteNode import DiscreteNode
 from GaussNode import GaussNode
+from MCMC import MarkovChainSampler
+from MCMC import GibbsTransitionModel

primo/reasoning/density/ProbabilityTable.py

 # -*- coding: utf-8 -*-
 import numpy
+import copy
 from primo.reasoning.density import Density
 
 class ProbabilityTable(Density):
@@ -69,8 +70,24 @@ class ProbabilityTable(Density):
     def marginalization(self, variable):
         raise Exception("Called unimplemented function")       
         
-    def reduction(self):
-        raise Exception("Called unimplemented function")
+    def reduction(self, evidence):
+        '''Returns a reduced version of this ProbabilityTable, evidence is a list of pairs.
+            Important: This node is not being changed!'''
+        reduced = ProbabilityTable()
+        reduced.variables = copy.copy(self.variables)
+        reduced.table = self.table
+        for node,value in evidence:
+
+            axis=reduced.variables.index(node)
+            position=node.value_range.index(value)
+            reduced.table = numpy.take(reduced.table,[position],axis=axis)
+            
+            reduced.table=reduced.table.squeeze()
+            reduced.variables.remove(node)
+            
+        return reduced
+            
+            
 
     def division(self, factor):
         raise Exception("Called unimplemented function")