particlefilters.cc.svn-base

Probabilistic graphical models in matlab.

	}
	
	while (time ? time_used < max_steps : gibbs_step < max_steps )
	{
		++gibbs_step;
		
		if (!time)
		{
			if ( (gibbs_step % (max_steps/100) == 0) || (max_steps<100) )
			{ 
				cout << "." << flush; 
			}
			
			current_experiment->callback_monitoring(gibbs_step);
		}
		
		else
		{
			if (time_used > last_progress_dot_time)
			{ 
				cout << "." << flush;
				last_progress_dot_time += double (max_steps)/100.0;
			}
			
			current_experiment->callback_monitoring(time_used);
		}
		
		// Set up looping over trees (components over which we can infer probabilities). WARNING: HERE IT MUST BE SPECIALLY CRAFTED CHAINS !
		
		vector < iterator_property_map < vector< vector < double > >::iterator, property_map < subgraph<DiscreteGraph>, vertex_index_t>::type > >::iterator weights_property_maps_it = weights_property_maps_vec.begin();
		vector < iterator_property_map < vector< vector < unsigned int > >::iterator, property_map < subgraph<DiscreteGraph>, vertex_index_t>::type > >::iterator positions_property_maps_it = positions_property_maps_vec.begin();
		
		for ( tie(current_tree, current_tree_end) = subgraph_g.children(); current_tree != current_tree_end; ++current_tree)
		{
			SubGraphVertexDescriptor root_node = vertex(0, *current_tree);
			
			for ( tie(u, u_end) = vertices(*current_tree); u != u_end; ++u)
			{
				if (out_degree(*u,*current_tree) == 2) // we found the "root" of our chain 
				{
					root_node = *u;
					break;
				}
			}
			
			// We get a sample from the joint using FFBS with PFs
			
			SubGraphVertexDescriptor tail_node = forward_filtering_implementation (*current_tree,  root_node, *positions_property_maps_it, *weights_property_maps_it);
			
			// This gives us the weights, and the particles, now run the backward smoothing part
			
			backward_smoothing_implementation (*current_tree,  tail_node, *positions_property_maps_it, *weights_property_maps_it, backward_samples);
			
			//  display_variables_graph(*current_tree);
			
			//++current_tree_weights_map;
			//++current_tree_positions_map;
			
			++positions_property_maps_it;
			++weights_property_maps_it;
		}
		
		// cout << "*********************************" << endl;
		
		time_used = global_timer.elapsed();
	}
	
	if (time)
	{
		cout << "We were able to generate " << gibbs_step << " complete tree-Gibbs samples." << endl;
	}
	
	// Here the Gibbs loop is done. That's all.
	
	// Need to clean up memory !!
}


template <>
void non_parametric_tree_gibbs_sampler < DiscreteGraph> (DiscreteGraph & g,  const unsigned int gibbs_steps, const unsigned int number_particles, const unsigned int backward_samples, bool time)
{	
	//typedef typename vertex_bundle_type<Graph>::type VertexClassPointer;
	//typedef typename iterator_traits<VertexClassPointer>::value_type VertexClass;
	
	property_map<DiscreteGraph, vertex_index_t>::type vertex_index_map = get(vertex_index, g);
	
	typedef graph_traits<DiscreteGraph>::vertex_iterator  VertexIterator;
	
	VertexIterator u, u_end;
	
	vector < Potential * > original_potentials (num_vertices(g));
	
	for (tie(u,u_end) = vertices (g); u!= u_end; ++u)
	{
		if ( out_degree(*u,g) != 1)
		{
			original_potentials[ get(vertex_index_map, *u)] = g[*u]->get_potential();
			
			ChainedSingleDiscretePotential * p = new ChainedSingleDiscretePotential (* g[*u]->get_random_variable());
			
			p->add_potential(static_cast < DiscretePotential *> ( g[*u]->get_potential() ) );
			
			g[*u]->set_potential(p);
		}
	}
	
	non_parametric_tree_gibbs_sampler_implementation ( g,  gibbs_steps, number_particles, backward_samples, time);
	
	for (tie(u,u_end) = vertices (g); u!= u_end; ++u)
	{
		if ( out_degree(*u,g) != 1)
		{
			delete g[*u]->get_potential();
			g[*u]->get_random_variable()->obtain_inference_from_prior_samples();
			g[*u]->set_potential( original_potentials[ get(vertex_index_map, *u)]);
		}
	}
	
}

// No specialization needed for non_parametric_tree_gibbs_sampler_implementation (Graph & g, const unsigned int max_steps, const unsigned int number_particles);

//*********************** Generic Implementation ******************************//

//extern Experiment < NPGraph> * current_experiment;

template <class Graph>
void non_parametric_tree_gibbs_sampler(Graph & g, const unsigned int gibbs_steps, const unsigned int number_particles, const unsigned int backward_samples = 1, bool time = false)
{	
	//typedef typename vertex_bundle_type<Graph>::type VertexClassPointer;
	//typedef typename iterator_traits<VertexClassPointer>::value_type VertexClass;
	
	typename property_map<Graph, vertex_index_t>::type vertex_index_map = get(vertex_index, g);
	
	typedef typename graph_traits<Graph>::vertex_iterator  VertexIterator;
	
	VertexIterator u, u_end;
	
	vector < Potential * > original_potentials (num_vertices(g));
	
	for (tie(u,u_end) = vertices (g); u!= u_end; ++u)
	{
		if ( out_degree(*u,g) != 1)
		{
			original_potentials[ get(vertex_index_map, *u)] = g[*u]->get_potential();
			
			ChainedSinglePotential * p = new ChainedSinglePotential ( g[*u]->get_random_variable()->get_index()); //g[*u]->get_potential()->clone();
			
			p->add_potential(g[*u]->get_potential() );
			
			g[*u]->set_potential(p);
		}
	}
	
	non_parametric_tree_gibbs_sampler_implementation ( g, gibbs_steps, number_particles, backward_samples, time);
	
	for (tie(u,u_end) = vertices (g); u!= u_end; ++u)
	{	
		if ( out_degree(*u,g) != 1)
		{
			delete g[*u]->get_potential();
			g[*u]->get_random_variable()->obtain_inference_from_prior_samples();
			g[*u]->set_potential( original_potentials[ get(vertex_index_map, *u)]);
		}
	}
	
}

template <class Graph>
void non_parametric_tree_gibbs_sampler_implementation(Graph & g, const unsigned int max_steps, const unsigned int number_particles, const unsigned int backward_samples, bool time)
{
	typedef typename graph_traits< Graph >::adjacency_iterator  BaseAdjacencyIterator;
	typedef typename graph_traits< Graph >::vertex_descriptor  BaseVertexDescriptor;
	
	typedef typename graph_traits< subgraph<Graph> >::vertex_iterator  SubGraphVertexIterator;
	typedef typename graph_traits< subgraph<Graph> >::vertex_descriptor  SubGraphVertexDescriptor;
	
	// typedef typename graph_traits<subgraph<Graph> >::edge_descriptor BaseEdgeDescriptor;
	
	typedef typename graph_traits<subgraph<Graph> >::edge_iterator SubGraphEdgeIterator;
	
	SubGraphVertexIterator u, u_end;
	BaseAdjacencyIterator v,v_end;
	BaseVertexDescriptor w;
	
	typename subgraph<Graph>::children_iterator current_tree, current_tree_end;
	
	subgraph <Graph> subgraph_g;
	
	copy_graph(g, subgraph_g);
	
	// If it is a simple chain, do this (no partitioning needed)
	
	//partition_simple_chain(subgraph_g);
	
	// Else do this for a MRF (with a PAIR and EQUAL number of rows / columns)
	
	partition_mrf_graph_into_chains(subgraph_g);
	
	//display_subgraph(subgraph_g);
	
	// Following block set up the correct internal potential for a variable given the edges linking it to the other variables (conditionned)
	
	list < vector < vector < double > > * > trees_weights_map;
	list < vector < vector < double > > * > trees_positions_map;
	
	for ( tie(current_tree, current_tree_end) = subgraph_g.children(); current_tree != current_tree_end; ++current_tree)
	{
		vector < double > particles_weights (number_particles);
		
		vector < vector < double > > * weights_vector = new  vector < vector < double > > ( num_vertices (*current_tree), particles_weights);
		
		trees_weights_map.push_back(weights_vector);
		
		vector < double > particles_positions (number_particles);
		
		vector < vector < double > > * positions_vector = new  vector < vector < double > > ( num_vertices (*current_tree), particles_positions);
		
		trees_positions_map.push_back(positions_vector);
		
		for (tie(u,u_end) = vertices (*current_tree); u!= u_end; ++u)
		{
			w = current_tree->local_to_global(*u);
			
			for (tie(v, v_end)= adjacent_vertices ( w, g ); v != v_end; ++v)
			{
				if (!current_tree->find_vertex(*v).second)
				{
					// Add that edge potential to the ChainedSinglePotential
					
					static_cast < ChainedSinglePotential * > (g[w]->get_potential())->add_potential( g[edge(w,*v,g).first] );
				}
			}
			
		}
		
	}
	
	// End
	
	
	list < vector < vector < double > > * >::iterator current_tree_weights_map = trees_weights_map.begin();
	list < vector < vector < double > > * >::iterator current_tree_positions_map = trees_positions_map.begin();
	
	vector < iterator_property_map < vector< vector < double > >::iterator, typename property_map < subgraph<Graph>, vertex_index_t>::type > > weights_property_maps_vec;
	vector < iterator_property_map < vector< vector < double > >::iterator, typename property_map < subgraph<Graph>, vertex_index_t>::type > > positions_property_maps_vec;
	
	
	for ( tie(current_tree, current_tree_end) = subgraph_g.children(); current_tree != current_tree_end; ++current_tree)
	{
		//iterator_property_map < vector< vector < double > >::iterator, typename property_map < subgraph<Graph>, vertex_index_t>::type > position_map = make_iterator_property_map( (*current_tree_positions_map)->begin(), get(vertex_index, *current_tree));
		
		//iterator_property_map < vector< vector < double > >::iterator, typename property_map < subgraph<Graph>, vertex_index_t>::type > weight_map = make_iterator_property_map( (*current_tree_weights_map)->begin(), get(vertex_index, *current_tree));
		
		positions_property_maps_vec.push_back(make_iterator_property_map( (*current_tree_positions_map)->begin(), get(vertex_index, *current_tree)));
		
		weights_property_maps_vec.push_back(make_iterator_property_map( (*current_tree_weights_map)->begin(), get(vertex_index, *current_tree)));
		
		++current_tree_weights_map;
		++current_tree_positions_map;
	}
	
	cout << "Running NP Chain Gibbs with " << max_steps << (time ? " seconds and " : " steps and ") << number_particles << " particles." << endl;
	
	// cout << "MapSize: " << trees_weights_map.size() << ", " << trees_positions_map.size() << endl;
	
	// Main loop after preparation stuff (Gibbs steps)
	
	fibonnacci_number_generator.seed(std::time(NULL));
	
	unsigned int gibbs_step(0);
	
	double time_used(0.0);
	
	double last_progress_dot_time = double (max_steps)/100.0;
	
	if (time)
	{
		global_timer.restart();
	}
	
	while (time ? time_used < max_steps : gibbs_step < max_steps )
	{
		++gibbs_step;
		
		if (!time)
		{
			if ( (gibbs_step % (max_steps/100) == 0) || (max_steps<100) )
			{ 
				cout << "." << flush; 
			}
		}
		
		else
		{
			if (time_used > last_progress_dot_time)
			{ 
				cout << "." << flush;
				last_progress_dot_time += double (max_steps)/100.0;
			}
			
			current_experiment->callback_monitoring(time_used);
		}
		
		// Set up looping over trees (components over which we can infer probabilities). WARNING: HERE IT MUST BE SPECIALLY CRAFTED CHAINS !
		
		typename vector < iterator_property_map < vector< vector < double > >::iterator, typename property_map < subgraph<Graph>, vertex_index_t>::type > >::iterator weights_property_maps_it = weights_property_maps_vec.begin();
		typename vector < iterator_property_map < vector< vector < double > >::iterator, typename property_map < subgraph<Graph>, vertex_index_t>::type > >::iterator positions_property_maps_it = positions_property_maps_vec.begin();
		
		for ( tie(current_tree, current_tree_end) = subgraph_g.children(); current_tree != current_tree_end; ++current_tree)
		{
			/*
			 SubGraphEdgeIterator e, e_end;
			 
			 tie (e,e_end) = edges(g);
			 
			 for (tie (e,e_end) = edges(g); e != e_end; ++e)
			 {
				 g[*e]->debug_display_state();
			 }
			 */
			
			SubGraphVertexDescriptor root_node = vertex(0, *current_tree);
			
			for ( tie(u, u_end) = vertices(*current_tree); u != u_end; ++u)
			{
				if (out_degree(*u,*current_tree) == 2) // we found the "root" of our chain 
				{
					root_node = *u;
					break;
				}
			}
			
			//cout << "Root node of FF is: " << g[root_node]->get_random_variable()->get_index() << endl;
			
			SubGraphVertexDescriptor tail_node = forward_filtering_implementation (*current_tree,  root_node, *positions_property_maps_it, *weights_property_maps_it);
			
			// This gives us the weights, and the particles, now run the backward smoothing part
			
			backward_smoothing_implementation (*current_tree,  tail_node, *positions_property_maps_it, *weights_property_maps_it, backward_samples);
			
			//  display_variables_graph(*current_tree);
			
			//++current_tree_weights_map;
			//++current_tree_positions_map;
			
			++positions_property_maps_it;
			++weights_property_maps_it;
		}
		
		// cout << "*********************************" << endl;