cplusplus_examples.svn-base

Probabilistic graphical models in matlab.

		u = current_tree->local_to_global(*v);
		
			for (tie(w, w_end)= adjacent_vertices ( u, g ); w != w_end; ++w)
			{
				if (!current_tree->find_vertex(*w).second)
				{
					// Add that edge
					
					cut_edges[u].push_back(edge(u, *w, g).first);
				}
			}
		
		}
	
	}
	
	// End 
	
	// Main loop after preparation stuff (Gibbs steps)
	
	for (unsigned int i = 0; i < max_steps; ++i)
	{
		
		cout << "Gibbs step " << i << endl;
	
	// Set up looping over trees (components over which we can infer probabilities)

		current_tree = tree_components.begin(); // hack because of bundled properties
	
		for (current_graph_tree = tree_graph_components.begin() ; current_graph_tree != tree_graph_components.end(); ++current_graph_tree)
		{
	
			// Must get the values of the adjacent edges and set them as potentials...
			// We will assume all such edges are stored somewhere... they are
			
			//Make product over ALL edges
			
			for (tie(v,v_end) = vertices (*current_graph_tree); v!= v_end; ++v)
			{

				for (unsigned int i = 0; i < ((*current_graph_tree)[*v])->get_random_variable()->get_number_possible_values(); ++i)
				{
					double a = 1.0;

					// Maybe we could somehow put the cut_edges list on the vertices. This would allow us to not all the time use the local_to_global
					
					for (typename list< EdgeDescriptor>::iterator it = cut_edges[current_tree->local_to_global(*v)].begin(); it != cut_edges[current_tree->local_to_global(*v)].end(); ++it)
					{
						a = a * g[*it].obtain_potential_value((*current_graph_tree)[*v]->get_random_variable()->get_index(), i);
					}
					
					// Change our SinglePotential so as to incorporate the changes
					
					((*current_graph_tree)[*v])->get_potential()->modify_potential_value(i, a);
				}
				
			}
			// Must run inference
			
			message_passing(*current_graph_tree); // this should be ok except initialization stuff
			
			// Must sample now on every variable, we can combine that with setting values
			// Must set values, this could (should) be done by a subroutine
			// Note: This will be done on the main graph
			
			for (tie(v,v_end) = vertices (*current_graph_tree); v!= v_end; ++v)
			{
			
			((*current_graph_tree)[*v])->get_random_variable()->sample();
			
				for ( typename list <EdgeDescriptor>::iterator it = cut_edges[current_tree->local_to_global(*v)].begin(); it != cut_edges[current_tree->local_to_global(*v)].end(); ++it)
				{
					g[*it].set_variable_value_by_index(  ( ((*current_graph_tree)[*v])->get_random_variable())->get_index(),  ((*current_graph_tree)[*v])->get_random_variable()->value );
				}
			}
			
			// Next line is a hack for no bundled properties
			
			current_tree++;
			
		}
		
		
	}
}




template <class Graph>
void tree_gibbs_sampler_bundled (Graph & g, unsigned int max_steps, unsigned int burnoff)
{
	typedef typename graph_traits<Graph>::vertex_iterator  VertexIterator;
//	typedef typename graph_traits< subgraph <Graph> >::vertex_iterator  SubGraphVertexIterator;
	typedef typename graph_traits<Graph>::adjacency_iterator  AdjacencyIterator;
	typedef typename graph_traits<Graph>::vertex_descriptor  VertexDescriptor;
	typedef typename graph_traits<Graph>::edge_descriptor EdgeDescriptor;
	typedef typename Graph::vertex_bundled  VertexClassPointer;
	typedef typename iterator_traits<VertexClassPointer>::value_type VertexClass;
	typedef typename VertexClass::random_variable_type RandomVariableClass;
	
	AdjacencyIterator w,w_end;
	
	//SubGraphVertexIterator v, v_end;
	
	VertexIterator v, v_end;
	
	VertexDescriptor u;
	vector < list < EdgeDescriptor> > cut_edges;
	
	typename vector < subgraph <Graph> >::iterator current_tree;

	
	//(*current_tree)


	//Get tree components, create cut_edges
	
	vector < subgraph <Graph> > tree_components = partition_graph_into_trees(g);
	
		
	current_tree = tree_components.begin();
	
	
	for (current_tree = tree_components.begin();current_tree !=tree_components.begin(); ++current_tree)
	{
		for (tie(v,v_end) = vertices (*current_tree); v!= v_end; ++v)
		{

		u = current_tree->local_to_global(*v);
		
			for (tie(w, w_end)= adjacent_vertices ( u, g ); w != w_end; ++w)
			{
				if (!current_tree->find_vertex(*w).second)
				{
					// Add that edge
					
					cut_edges[u].push_back(edge(u, *w, g).first);
				}
			}
		
		}
	
	}
	
	// Set up looping over steps
	
	for (unsigned int i = 0; i < max_steps; ++i) // main loop (Gibbs steps)
	{
	
	
	// Set up looping over trees (components over which we can infer probabilities)
		
		for (current_tree = tree_components.begin();current_tree !=tree_components.begin(); ++current_tree)
		{
	
			// Must get the values of the adjacent edges and set them as potentials...
			// We will assume all such edges are stored somewhere... they are
			
			//Make product over ALL edges
			
			for (tie(v,v_end) = vertices (*current_tree); v!= v_end; ++v)
			{

			
				for (unsigned int i = 0; i < ((*current_tree)[*v])->get_random_variable()->get_number_possible_values(); ++i)
				{
					double a = 1.0;
					
					for (typename list< EdgeDescriptor>::iterator it = cut_edges[current_tree->local_to_global(*v)].begin(); it != cut_edges[current_tree->local_to_global(*v)].end(); it++
					)
					{
						a = a * g[*it].obtain_potential_value(g[*v]->get_random_variable()->get_index(), i);
					}
						
					//modify_potential_value(i, a);  // HAS TO BE UNCOMMENTED
				}
				
			}
			// Must run inference
			
			message_passing(*current_tree); // this should be ok except initialization stuff
			
			//Now we got to store the results somewhere...
			
			
			
			// Must sample now on every variable, we can combine that with setting
			
			// Must set values
			
			for (tie(v,v_end) = vertices (*current_tree); v!= v_end; ++v)
			{
			
			(*current_tree[*v])->get_random_variable()->sample();
			
				for ( typename list <EdgeDescriptor>::iterator it = cut_edges[current_tree->local_to_global(*v)].begin(); it != cut_edges[current_tree->local_to_global(*v)].end(); ++it)
				{
					g[*it].set_variable_value_by_index(  ( ((*current_tree)[*v])->get_random_variable())->get_index(), ((*current_tree[*v])->get_random_variable())->value);
				}
			}
			
		}
		
		
	}
}


	
// Following function needs testing...
/*
double DiscretePotential::internal_marginalize_over_variable(unsigned int index)
{
	
	unsigned int left_block_size,right_block_size,right_blocks, fixed_value;
	
	
	// Right_blocks is the number of blocks we have to "the right", it is in fact the number of times we must sum in the outer loop
	// Right_block_size is the size of the blocks to the right... we add that once we are over in the inner loop
	
	if (index == number_of_variable_values.size()-1)
	{
		right_blocks = 1;
		right_block_size = 0; // Not used in this case
	}
	else
	{
		right_blocks = number_of_variable_values[index+1];
		right_block_size = help_vector[index+1];
	}
	
	// The fixed value represents how much we must add constantly to read the point where "the sum is frozen"
	// Left block size represents the size of the blocks to the left, ie, 1 if we are not summing over the first variable
	
	left_block_size = help_vector[index];
	fixed_value = variable_values[index]*left_block_size;
	
	double sum = 0;
	unsigned int a = fixed_value;
	
	for (unsigned int i = 0; i < right_blocks; i++)
	{
		
		for (unsigned int j = 0; j < left_block_size; j++)
		{
			sum += potential_table[a+j];
		}
		
		a += right_block_size;
		
	}
	
	return sum;
	
}
*/

/* untested function below

double internal_sum_over_variable(int index)
{
	// We mustn't disturb the state of the Potential, so we save the variable_value 
	
	int backup_value = variable_values[index];
	variable_values[index] = 0;
	
	int a = inner_product(variable_values.begin(), variable_values.end(), help_vector.begin(), 0);
	
	int increment = help_vector[index];
	
	double sum = 0;
	
	for (int i = 0; i < ; i++)
	{
		a += i*increment;
		
	sum += potential_table[a];
	}
	
	variable_values[index] = backup_value;
	
	return a;
}

*/