old_code.svn-base

Probabilistic graphical models in matlab.

	l.push_back(25);
	
	combine_nodes(*g, l, 35);
	
	cout << "Exact Inference " << endl;
	
	message_passing(*g);
	
	uncombine_node (*g, 31);
	uncombine_node (*g, 32);
	uncombine_node (*g, 33);
	uncombine_node (*g, 34);
	uncombine_node (*g, 35);
	*/
	
	/*
	cout << "Gibbs Approx." << endl;
	
	gibbs_sampler( *g, 100000 , 10000);
	*/
	
	/*
	cout << "Tree Gibbs Approx." << endl;
	
	tree_gibbs_sampler( *g, 10000 , 10);
	
	
	display_variables_graph(*g);
	*/
	
	//(*g)[new_node]->get_potential()->display_table();
	
	/*
	 
	 EdgeIterator e, e_end;
	 
	 for (tie (e, e_end) = edges (*g); e != e_end; ++e)
	 {
		 cout << "Edge linking " << source(*e,*g)+1 << " to " << target(*e,*g)+1 << ": ";
		 (*g) [*e].display_table();
	 }
	 
	 // seems to be ok
	 
	 */
	
	//cout << "Combine Node " << new_node << endl;
	
	//display_graph(*g);
	
	//cout << "Out " << endl;
	
	//message_passing(*g);
	
	/*
	 tree_gibbs_sampler( *g, 50 , 10);
	 
	 cout << "End Tree Gibbs " << endl;
	 
	 
	 */


// End Gibbs sampler

/* Algorithm not used yet

template <class Graph>
double obtain_probability (Graph & G)
{	
	typedef typename Graph::edge_bundled  EdgeClass;
	typedef typename Graph::vertex_bundled  VertexClass;
	
	typename graph_traits<Graph>::vertex_iterator  begin_vertex, end_vertex;
	
	tie (begin_vertex, end_vertex) = vertices(G);
	
	typename graph_traits<Graph>::edge_iterator  begin_edge, end_edge;
	
	tie (begin_edge, end_edge) = edges(G);
	
	// Next line makes a product over all the edge's potentials
	
	double a = make_product_over ( begin_edge, end_edge, adapt_graph_edges <double> (G, mem_fun_ref(&EdgeClass::get_potential_value) ) );
	
	// Next line makes a product over all the potentials contained in the Vertexes
	
	a = a * make_product_over ( begin_vertex, end_vertex, adapt_graph_vertexes <double> (G, boost::bind(&Potential::get_potential_value, bind(&VertexClass::get_potential, _1) ) ) );
	
	return a;
}

*/
// Obsolete
/*
template < class Graph >
void forward_filtering_backward_smoothing(Graph & g, const unsigned int number_particles)
{
	// Initialize particles, create an external map for all weights and particles
	
	typedef typename graph_traits<Graph>::vertex_descriptor VertexDescriptor;
	
	typename property_map <Graph, vertex_index_t>::type vertex_index_map = get(vertex_index, g);
	
	vector < double > particles_weights (number_particles);
	
	vector < vector < double > > weights_vector ( num_vertices (g), particles_weights);
	
	vector < double > particles_positions (number_particles);
	
	vector < vector < double > > positions_vector ( num_vertices (g), particles_positions);
	
	iterator_property_map < vector< vector < double > >::iterator, typename property_map <Graph, vertex_index_t>::type > position_map = make_iterator_property_map(positions_vector.begin(), get(vertex_index,g));
	
	VertexDescriptor root_node;
	
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	
	for ( tie(u, u_end) = vertices(g); u != u_end; ++u)
	{
		if (out_degree(*u,g) == 2) // we found the "root" of our chain 
		{
			root_node = *u;
			break;
		}
	}
	
	//cout << "Root node of FF is: " << root_node << " (" << g[root_node]->get_random_variable()->get_index() << ")" << endl;
	
	//cout << "Root node of FF is: " << g[root_node]->get_random_variable()->get_index() << endl;
	
	VertexDescriptor tail_node = forward_filtering_implementation (g,  root_node, position_map, make_iterator_property_map(weights_vector.begin(), get(vertex_index,g)));
	
	//cout << "Finished FF" << endl;
	
	//tail_node = tail_node + 1;
	
	// This gives us the weights, and the particles, now run the backward smoothing part
	
	backward_smoothing_implementation (g,  tail_node, make_iterator_property_map(positions_vector.begin(), get(vertex_index,g)), make_iterator_property_map(weights_vector.begin(), get(vertex_index,g)));
	
	//cout << "Finished BS" << endl;
}
*/



template <class Graph>
void set_7_7_mrf_observations (Graph & g)
{
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	typename graph_traits<Graph>::adjacency_iterator v, v_end;
	
	const unsigned int increment(49);
	
	static_cast < DiscreteRandomVariable * > (g[0+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[1+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[2+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[3+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[4+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[5+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[6+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[7+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[8+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[9+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[10+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[11+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[12+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[13+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[14+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[15+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[16+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[17+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[18+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[19+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[20+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[21+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[22+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[23+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[24+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[25+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[26+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[27+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[28+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[29+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[30+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[31+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[32+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[33+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[34+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[35+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[36+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[37+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[38+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[39+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[40+increment]->get_random_variable())->observe_variable_value( 0);
	static_cast < DiscreteRandomVariable * > (g[41+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[42+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[43+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[44+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[45+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[46+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[47+increment]->get_random_variable())->observe_variable_value( 1);
	static_cast < DiscreteRandomVariable * > (g[48+increment]->get_random_variable())->observe_variable_value( 1);
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		
		if (out_degree(*u,g) == 1)
		{
			tie (v, v_end) = adjacent_vertices(*u, g);
			
			g[edge(*u,*v,g).first]->set_variable_value( *(g[*u]->get_random_variable()));
		}
		
	}
}







// Not used anymore

template < >
void forward_filtering_backward_smoothing < subgraph<DiscreteGraph> > (subgraph<DiscreteGraph> & g, const unsigned int number_particles)
{
	// Initialize particles, create an external map for all weights and particles
	
	typedef graph_traits< subgraph<DiscreteGraph> >::vertex_descriptor VertexDescriptor;
	
	property_map < subgraph<DiscreteGraph> , vertex_index_t>::type vertex_index_map = get(vertex_index, g);
	
	vector < double > particles_weights (number_particles);
	
	vector < vector < double > > weights_vector ( num_vertices (g), particles_weights);
	
	vector < unsigned int > particles_positions (number_particles);
	
	vector < vector < unsigned int > > positions_vector ( num_vertices (g), particles_positions);
	
	iterator_property_map < vector< vector < unsigned int > >::iterator, property_map < subgraph<DiscreteGraph>, vertex_index_t>::type > position_map = make_iterator_property_map(positions_vector.begin(), get(vertex_index,g));
	iterator_property_map < vector< vector < double > >::iterator, property_map < subgraph<DiscreteGraph> , vertex_index_t>::type > weight_map = make_iterator_property_map(weights_vector.begin(), get(vertex_index,g));
	
	VertexDescriptor root_node = vertex(0,g);
	
	graph_traits<subgraph<DiscreteGraph> >::vertex_iterator u, u_end;
	
	for ( tie(u, u_end) = vertices(g); u != u_end; ++u)
	{
		if (out_degree(*u,g) == 2) // we found the "root" of our chain 
		{
			root_node = *u;
			break;
		}
	}
	
	//cout << "Root node of FF is: " << root_node << " (" << g[root_node]->get_random_variable()->get_index() << ")" << endl;
	
	//cout << "Root node of FF is: " << g[root_node]->get_random_variable()->get_index() << endl;
	
	VertexDescriptor tail_node = forward_filtering_implementation (g,  root_node, position_map, weight_map);
	
	//cout << "Finished FF" << endl;
	
	//tail_node = tail_node + 1;
	
	// This gives us the weights, and the particles, now run the backward smoothing part
	
	backward_smoothing_implementation (g,  tail_node, position_map, weight_map);
	
	//cout << "Finished BS" << endl;
}

	// We sample the tail node variable
	
	vector <double> & last_weights = get (weights, tail_node);
	
	// We normalize the last weights, it is not really needed though as if resampling always take place, their sum will be equal to 1 already
	
	normalize_over (last_weights.begin(), last_weights.end(), pf, pf2);
	
	
	
		if ( argc < 3 )
	{
		cout << "Give more args" << endl;
		return 0;
	}
	
	istringstream st (argv[1]);
	
	st >> mrf_rows;
	
	st.str (argv[2]);
	
	st.clear();
	
	st >> mrf_columns;
	
	//cout << "Rows, Columns: " << mrf_rows << " " << mrf_columns << endl;
	 
	string gr_tr = "discrete";
	
	if ( gr_tr.compare( argv[ 3] ) == 0)
	{
		//DiscreteGraph * g = make_mrf_graph < DiscreteGraph > ( mrf_rows, mrf_columns);
		//DiscreteGraph * g = make_simple_chain_graph <DiscreteGraph > (20);
		
		//current_discrete_experiment = new Experiment <DiscreteGraph> ( * g, 1, 300, true, true);
		
		current_discrete_experiment = new Experiment <DiscreteGraph> ();
		
		string gr_tr = "ground";
		
		if ( gr_tr.compare( argv[ 4] ) == 0)
		{

			if ( argc < 9 )
			{
				cout << "Give more args for ground: method time particles export_file_name" << endl;
				return 0;
			}
			

			current_discrete_experiment->set_ground_truth_method(argv[ 5]);

			
			unsigned int time,  particles;
			
			st.str (argv[6]);
			st.clear();
			st >>time;
			
			st.str (argv[7]);
			st.clear();
			st >> particles;
			
			//current_discrete_experiment->set_running_time(time);
			current_discrete_experiment->set_running_particles(particles);
			
			current_discrete_experiment->obtain_ground_truth(argv[8], time);
			
		}
		
		string run_str = "run";
		
		if ( run_str.compare( argv[ 4] ) == 0)
		{
			if ( argc < 10 )
			{
				cout << "Give more args for run: steps/timer time burnoff particles true_means_file produced_files_extension" << endl;
				return 0;
			}
			
			 run_str = "steps";
			
			if ( run_str.compare( argv[ 5] ) == 0)
			{
				current_discrete_experiment->set_timer_method(false);
			}
			
			
			unsigned int time, burnoff, particles;
			
			st.str (argv[6]);
			st.clear();
			st >>time;
			
			st.str (argv[7]);
			st.clear();