graphalgorithms.cc.svn-base

Probabilistic graphical models in matlab.

	typename subgraph<Graph>::children_iterator current_tree, tree_end;
	
	vector < vector < double> > sum_marginals (num_vertices(g));
	vector < Potential * > original_potentials (num_vertices(g));
	
	fibonnacci_number_generator.seed(std::time(NULL));
	
	initialize_graph_variables_random(g);
	
	// We now partition the graph. Warning: we send directly the subgraph type when we send to partition_graph_into_trees, not the graph
	
	//display_graph(g);
	
	subgraph <Graph> subgraph_g;
	copy_graph(g, subgraph_g);
	partition_graph_into_trees(subgraph_g);
	
	//display_subgraph(subgraph_g);
	
	// Here we setup stuff based on the partition. This is a specialized templated function.
	
	setup_tree_gibbs_sampler(subgraph_g, sum_marginals, original_potentials);
	
	for ( tie (current_tree, tree_end) = subgraph_g.children(); current_tree != tree_end; ++current_tree)
	{
		set_temp_observations(* current_tree);
	}
	
	
	cout << "Running MCMC Tree Gibbs Sampling with " << max_steps <<  (time ? " seconds."  : " steps." ) << endl;
	
	unsigned int current_gibbs_step (0);
	
	bool record = false;
	
	double last_progress_dot_time = double (max_steps)/100.0;
	
	double time_used(0.0);
	
	if (time)
	{
		global_timer.restart();
	}
	
	
	// Main loop after preparation stuff (Gibbs steps)
	
	while (time ? time_used < max_steps : current_gibbs_step < max_steps)
	{
		
		// Set up looping over trees (components over which we can infer probabilities)
		
		for ( tie (current_tree, tree_end) = subgraph_g.children(); current_tree != tree_end; ++current_tree)
		{
			//cout << "*************************************" << endl;
			
			// Initialization of the tree (we clear previous messages, etc...).
			
			initialize_graph_contents(* current_tree);
			
			// Must run inference
			
			belief_propagation_with_backward_sampling(* current_tree);
			
			//  display_variables_graph(*current_tree);
			
			// We have now obtained a sample from the joint. The following will also take care of RB estimates
			
			for (tie(u,u_end) = vertices (*current_tree); u!= u_end; ++u)
			{	
				// We could remove the following line to prevent use of RB, and only use Gibbs counting
				
				(*current_tree)[*u]->add_rao_blackwell_estimates(sum_marginals[current_tree->local_to_global(*u)]);
				
				(*current_tree)[*u]->set_variable_to_joint_sampled_value();
			}
			
			set_temp_observations((*current_tree));
		}
		
		time_used = global_timer.elapsed();
		++current_gibbs_step;
		
		if (!time)
		{
			if ( (current_gibbs_step % (max_steps/100) == 0) || (max_steps < 100) )
			{ 
				cout << "." << flush; 
			}
			
			if ( current_gibbs_step > burnoff_steps )
			{ 
				record = true; 
			}
			
			if (monitoring && record)
			{
				// Following is not portable...
				
				if (rao_blackwell)
				{
					for (tie (u, u_end) = vertices (g); u != u_end; ++u)
					{
						g[*u]->set_rao_blackwell_estimates(sum_marginals[*u]);
					}
				}
				else
				{
					for (tie (u, u_end) = vertices (g); u != u_end; ++u)
					{
						g[*u]->compute_marginals_from_samples();
					}
				}
				
				current_experiment->callback_monitoring(current_gibbs_step);
			}
			
			
		}
		
		else
		{
			if (time_used > last_progress_dot_time)
			{ 
				cout << "." << flush;
				last_progress_dot_time += double (max_steps)/100.0;
			}
			
			if ( time_used > burnoff_steps )
			{ 
				record = true; 
			}
			
			
			if (monitoring && record)
			{
				// Following is not portable...
				
				if (rao_blackwell)
				{
					for (tie (u, u_end) = vertices (g); u != u_end; ++u)
					{
						g[*u]->set_rao_blackwell_estimates(sum_marginals[*u]);
					}
				}
				else
				{
					for (tie (u, u_end) = vertices (g); u != u_end; ++u)
					{
						g[*u]->compute_marginals_from_samples();
					}
				}
				
				current_experiment->callback_monitoring(time_used);
				
				
			}
			
			
			
		}
	}
	
	cout << endl;
	
	// Following is not portable to other types of graphs than adjacency_list...
	
	if (rao_blackwell)
	{
		for (tie (u, u_end) = vertices (g); u != u_end; ++u)
		{
			g[*u]->set_rao_blackwell_estimates(sum_marginals[*u]);
		}
	}
	
	else
	{
		for (tie (u, u_end) = vertices (g); u != u_end; ++u)
		{
			g[*u]->compute_marginals_from_samples();
		}
	}
	
	if (time)
	{
		cout << "We were able to generate " << current_gibbs_step << " Tree Gibbs samples." << endl;
	}
	
	// Deallocating memory
	
	cleanup_tree_gibbs_sampler(g, original_potentials);
}

// Implementation of a Simple Gibbs Sampler algorithm, which samples nodes one at a time.

template <class Graph>
void gibbs_sampler(Graph & g, const unsigned int max_steps, const unsigned int burnoff_steps, const bool time = false)
{
	typename graph_traits<Graph>::vertex_iterator  u, u_end;
	
	cout << "Running Gibbs Sampler with " << max_steps <<  (time ? " seconds and "  : " steps and " ) << burnoff_steps << " burnoff steps." << endl;
	
	assert ( burnoff_steps < max_steps );
	
	fibonnacci_number_generator.seed(std::time(NULL));
	
	initialize_graph_variables_random(g);
	
	// Now we setup the proposals, these should not change
	
	for ( tie (u, u_end) = vertices(g) ; u != u_end; ++u)
	{	
		// Here we do some graph dependent setup for setting the proposal
		
		setup_gibbs_proposal(g, *u);
	}
	
	unsigned int current_gibbs_step = 0;
	
	bool record = false;
	
	double last_progress_dot_time = double (max_steps)/100.0;
	
	double time_used(0.0);
	
	if (time)
	{
		global_timer.restart();
	}
	
	// Begin of the main Gibbs loop
	
	while (time ? time_used < max_steps : current_gibbs_step <= max_steps)
	{
		
		for ( tie (u, u_end) = vertices(g) ; u != u_end; ++u)
		{	
			// Here come the generic virtual method
			
			g[*u]->sample_from_proposal(record);
		}
		
		++current_gibbs_step;
		time_used = global_timer.elapsed();
		
		if (!time)
		{
			if ( (current_gibbs_step % (max_steps/100) == 0) || (max_steps < 100) )
			{ 
				cout << "." << flush; 
			}
			
			
			if (monitoring && record)
			{
				for (tie (u, u_end) = vertices (g); u != u_end; ++u)
				{
					g[*u]->compute_marginals_from_samples();
				}
				
				current_experiment->callback_monitoring(current_gibbs_step);
			}
			
			
			if ( current_gibbs_step > burnoff_steps )
			{ 
				record = true; 
			}
		}
		
		else
		{
			if (time_used > last_progress_dot_time)
			{ 
				cout << "." << flush;
				last_progress_dot_time += double (max_steps)/100.0;
			}
			
			if (monitoring && record)
			{
				for (tie (u, u_end) = vertices (g); u != u_end; ++u)
				{
					g[*u]->compute_marginals_from_samples();
				}
				
				current_experiment->callback_monitoring(time_used);
			}
			
			if ( time_used > burnoff_steps )
			{ 
				record = true; 
			}
		}
		
	}
	
	// End of the Gibbs loop
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		g[*u]->compute_marginals_from_samples();
	}
	
	cout << endl;
	
	if (time)
	{
		cout << "We were able to generate " << current_gibbs_step << " Gibbs samples." << endl;
	}
}

// The message_passing algorithm consists in two depth first search passes.

template <class Graph>
void belief_propagation(Graph & g)
{	
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	
	// We record every node as its own parent
	
	for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
	{
		g[*u]->parent_node = *u;
	}
	
	depth_first_search(g, visitor(make_dfs_visitor(std::make_pair(Collect < Graph > (g),  RecordParent  < Graph > (g) ) ) ));	
	depth_first_search(g, visitor(make_dfs_visitor(Distribute< Graph > (g) ) ) );
	
	for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
	{
		g[*u]->compute_marginals_from_messages();
	}
}

template <class Graph>
void belief_propagation_with_backward_sampling(Graph & g)
{	
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	
	// We record every node as its own parent
	
	for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
	{
		g[*u]->parent_node = *u;
	}
	
	//cout << "Starting Collect" << endl;
	
	depth_first_search(g, visitor(make_dfs_visitor(std::make_pair(Collect < Graph > (g),  RecordParent  < Graph > (g) ) ) ));	
	
	//cout << "Starting B Sampling" << endl;
	
	depth_first_search(g, visitor(make_dfs_visitor(BackwardSampling< Graph > (g) ) ) );
	
	//cout << "Starting Distribute" << endl;
	
	remove_temp_observations(g);
	
	depth_first_search(g, visitor(make_dfs_visitor(Distribute< Graph > (g) ) ) );
	
	for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
	{
		g[*u]->compute_marginals_from_messages();
	}
	
	/*
	 Attempt at NOT using FFBS for debug purposes
	 
	for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
	{
		g[*u]->sample();
	}
	 */
}

template <class Graph>
void loopy_belief_propagation(Graph & g, const unsigned int max_steps, const bool time = false)
{		
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	typename graph_traits<Graph>::adjacency_iterator v, v_end;
	
	// Main loop, (loopy steps)
	
	cout << "Running Loopy Belief Propagation with " << max_steps <<  (time ? " seconds."  : " steps." ) << endl;
	
	double last_progress_dot_time = double (max_steps)/100.0;
	
	double time_used(0.0);
	
	unsigned int current_loopy_step(0);
	
	if (time)
	{
		global_timer.restart();
	}
	
	while (time ? time_used < max_steps : current_loopy_step < max_steps)
	{
		
		for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
		{
			for (tie (v, v_end) = adjacent_vertices(*u, g); v != v_end; ++v)
			{
				// Here we do some graph dependent setup for sending the message
				
				setup_message_computation (g, *u , *v);
				
				g[*v]->update_message ( g[*u]->send_message(* g[*v] ) );
			}
		}
		
		time_used = global_timer.elapsed();
		++current_loopy_step;
		
		if (!time)
		{
			if ( (current_loopy_step % (max_steps/100) == 0) || (max_steps < 100) )
			{ 
				cout << "." << flush; 
			}
			
			if (monitoring)
			{
				for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
				{
					g[*u]->compute_marginals_from_messages();
				}
				
				current_experiment->callback_monitoring(current_loopy_step);
			}
		}
		
		else
		{
			if (time_used > last_progress_dot_time)
			{ 
				cout << "." << flush;
				last_progress_dot_time += double (max_steps)/100.0;
			}
			
			
			if (monitoring)
			{
				for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
				{
					g[*u]->compute_marginals_from_messages();
				}
				
				current_experiment->callback_monitoring(time_used);
			}
		}
		
	}
	
	// We are finished, let's get the marginals (beliefs)
	
	if (time)
	{
		cout << "\nWe were able to generate " << current_loopy_step << " Loopy steps." << endl;
	}
	
	for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
	{
		g[*u]->compute_marginals_from_messages();
	}
}


template <class Graph>
void forward_filtering_backward_sampling(Graph & g)
{	
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	
	// We record every node as its own parent
	
	for ( tie (u, u_end) = vertices(g); u != u_end; ++u)
	{
		g[*u]->parent_node = *u;
	}
	
	depth_first_search(g, visitor(make_dfs_visitor(std::make_pair(Collect < Graph > (g),  RecordParent  < Graph > (g) ) ) ));
	
	depth_first_search(g, visitor(make_dfs_visitor(BackwardSampling< Graph > (g) ) ) );
}