Extended examples/MCMC.py, improved MarkovChainSampler to use less memory...

Extended examples/MCMC.py, improved MarkovChainSampler to use less memory (lukas idea), started working on division

examples/MCMC.py

@@ -4,6 +4,7 @@ from  primo.core import BayesNet
 from  primo.reasoning import DiscreteNode
 from  primo.reasoning import MarkovChainSampler
 from primo.reasoning import GibbsTransitionModel
+from primo.reasoning.density import ProbabilityTable
 import numpy
 
 #Construct some simple BayesianNetwork
@@ -28,7 +29,22 @@ transition_model = GibbsTransitionModel()
 
 mcs = MarkovChainSampler()
 initial_state={burglary:"Safe",alarm:"Silent"}
-chain = mcs.generateMarkovChain(bn, 20, transition_model, initial_state)
+chain = mcs.generateMarkovChain(bn, 1000000, transition_model, initial_state)
+
+#for c in chain:
+#    print c
+
+
+pt = ProbabilityTable()
+pt.add_variable(burglary)
+pt.add_variable(alarm)
+pt.to_jpt_by_states(chain)
+print "----joint-probability----"
+print pt
+print "----burglary----"
+print pt.marginalization(alarm)
+print "----alarm----"
+print pt.division(burglary)
+
+bn.draw()
 
-for c in chain:
-    print c

primo/reasoning/MCMC.py

 import random
-
+import copy
 def weighted_random(weights):
     counter = random.random() * sum(weights)
     for i,w in enumerate(weights):
@@ -57,7 +57,7 @@ class MarkovChainSampler(object):
         
     def generateMarkovChain(self, network, time_steps, transition_model, initial_state):
         state=initial_state
-        for t in range(time_steps):
+        for t in xrange(time_steps):
             yield state
             state=transition_model.transition(network, state)
         

primo/reasoning/density/ProbabilityTable.py

@@ -10,16 +10,11 @@ class ProbabilityTable(Density):
 
     def __init__(self):
         super(ProbabilityTable, self).__init__()
-
-        #self.owner = owner
-        #self.variables = [owner]
-
-        #size_of_range = len(owner.value_range)
-        #self.table = numpy.ones(size_of_range) / size_of_range
-
         self.variables = []
         self.table = numpy.array(0)
 
+            
+
     def get_table(self):
         return self.table
 
@@ -43,9 +38,21 @@ class ProbabilityTable(Density):
         self.table = table
         self.variables = nodes
 
+    def to_jpt_by_states(self, samples):
+        '''This method uses a list of variable-instantiations to change this nodes internal
+            table to represent a joint probability table constructed from the given samples.
+            The Argument samples is a list of pairs (RandomNode, value).'''
+        
+        for state in samples:
+            index = self.get_cpt_index(state.items())
+            self.table[index]=self.table[index]+1
+
+        self.normalize_as_jpt()
+            
+
     def set_probability(self, value, node_value_pairs):
         index = self.get_cpt_index(node_value_pairs)
-        self.table[tuple(index)]=value
+        self.table[index]=value
 
     def get_cpt_index(self, node_value_pairs):
         nodes, values = zip(*node_value_pairs)
@@ -54,7 +61,7 @@ class ProbabilityTable(Density):
             index_in_values_list = nodes.index(node)
             value = values[index_in_values_list]
             index.append(node.value_range.index(value))
-        return index
+        return tuple(index)
 
 
     def is_normalized_as_cpt(self,owner):
@@ -67,9 +74,9 @@ class ProbabilityTable(Density):
     def is_normalized_as_jpt(self):
         return numpy.sum(self.table) == 1.0
 
-    def normalized_as_jpt(self):
-        '''This method returns a new ProbabilityTable which is a normalized version of this ProbabilityTable (As a joint probability).'''
-        return self.table/numpy.sum(self.table)
+    def normalize_as_jpt(self):
+        '''This method normalizes this ProbabilityTable so it represents a valid joint probability table'''
+        self.table=self.table*1.0/numpy.sum(self.table)
 
     def multiplication(self, inputFactor):
         '''This method returns a unified ProbabilityTable which contains the variables of both; the inputFactor
@@ -164,7 +171,9 @@ class ProbabilityTable(Density):
 
 
     def division(self, factor):
-        raise Exception("Called unimplemented function")
+        '''Returns a new ProbabilityTable which is the result of dividing this one by the one given
+            with the argument factor'''
+        
 
     def __str__(self):
         return str(self.table)