inputfunctions.cc.svn-base

Probabilistic graphical models in matlab.

#include <InputFunctions.h>

#include <iostream>
#include <fstream>

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/connected_components.hpp>

#include <Generic.h>
#include <Potential.h>
#include <RandomVariable.h>
#include <Node.h>
#include <OutputGraph.h>

using namespace std;
using namespace boost;

// Explicit template instantiations //

typedef property< vertex_index_t, unsigned int> VertexProperty;
typedef property< edge_index_t, unsigned int> EdgeProperty;

typedef boost::property< boost::edge_index_t, unsigned int, DiscretePotential *> DiscretePairwiseEdgeProperty;
typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS, MessageNode<DiscreteRandomVariable, boost::adjacency_list_traits<boost::vecS, boost::vecS, boost::undirectedS>::vertex_descriptor> *, 	DiscretePairwiseEdgeProperty > DiscretePairwiseGraph;

typedef property< edge_index_t, unsigned int, ContinuousPotential *> NPEdgeProperty;
typedef adjacency_list<vecS, vecS, undirectedS, GraphicalModelNode *, NPEdgeProperty > NPGraph;

// End typedefs

//template DiscretePairwiseGraph * make_random_graph <DiscretePairwiseGraph> (const unsigned int, const double);
//template NPGraph * make_random_graph <NPGraph> (const unsigned int, const double);

template NPGraph * make_mrf_graph <NPGraph> (const unsigned int, const unsigned int);
template NPGraph * make_simple_chain_graph <NPGraph> (const unsigned int);
template DiscretePairwiseGraph * make_simple_chain_graph <DiscretePairwiseGraph> (const unsigned int);

//template void randomize_discrete_pairwise_graph <DiscretePairwiseGraph> (DiscretePairwiseGraph &);

//template void randomize_discrete_pairwise_symmetric_graph <DiscretePairwiseGraph> (DiscretePairwiseGraph &);
//template void fill_discrete_pairwise_symmetric_graph <DiscretePairwiseGraph>(DiscretePairwiseGraph &);

template void randomize_non_parametric_pairwise_graph <NPGraph> (NPGraph &);

template void export_computed_means < NPGraph> (string, NPGraph &);
template void read_reference_means < NPGraph> (string, NPGraph &);
template void set_observations < DiscretePairwiseGraph> (DiscretePairwiseGraph &);
template void set_observations < NPGraph> (NPGraph &);
template void initialize_random_variables < NPGraph> (NPGraph &);
template void initialize_random_variables < DiscretePairwiseGraph> (DiscretePairwiseGraph &);

template void set_special_mrf_observations < NPGraph> (NPGraph & g);

// End of Explicit template instantiations //

template void create_pairwise_square_lattice < DiscretePairwiseGraph> (const unsigned int, const unsigned int, DiscretePairwiseGraph &);
template void create_random_graph < DiscretePairwiseGraph> (const double, DiscretePairwiseGraph &);

//****************** set_vertex_index_map Implementation ********************* //

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

template <class Graph, class PropertyMapType>
set_vertex_index_map <Graph, PropertyMapType> ::set_vertex_index_map(Graph &)
{
}

template <class Graph, class PropertyMapType>
void set_vertex_index_map <Graph, PropertyMapType> ::operator()()
{
}

//*********** Partial Specialization of Class set_vertex_index_map ************//

template <class Graph>
set_vertex_index_map <Graph, boost::lvalue_property_map_tag>::set_vertex_index_map (Graph & a): g(a)
{
}

template <class Graph>
void set_vertex_index_map <Graph, boost::lvalue_property_map_tag>::operator () ()
{
	typename property_map<Graph, vertex_index_t>::type vertex_index_map = get(vertex_index, g);
	typename boost::graph_traits<Graph>::vertex_iterator v, v_end;
	
	unsigned int a = 0;
	
	for ( tie(v,v_end) = vertices (g); v != v_end; ++v)
	{	
		
		put(vertex_index_map , *v, a); // Only works if the property_map is mutable
		
		++a;
	}
	
}

//****************** End of set_vertex_index_map Implementation ********************* //

template <class Graph>
void create_random_graph ( const double density, Graph & g)
{
	typename property_map<Graph, vertex_index_t>::type vertex_index_map = get(vertex_index, g);
	
	typename boost::graph_traits<Graph>::vertex_descriptor t, u;
	
	double random_number(0.0);
	
	const unsigned int number_of_nodes = num_vertices(g);
	
	typename boost::graph_traits<Graph>::vertex_iterator v, v_end, w;
	
	for ( tie(v,v_end) = vertices (g); v != v_end; ++v)
	{
		w= v;
		
		for ( ++w; w != v_end; ++w)
		{
			random_number = uniform_distribution_0_1();
			
			if ( random_number < density)
			{
				add_edge(*v,*w,g);
			}
		}
	}

	// Below ensures that we only have ONE connected component in our Random Graph.
	
	std::vector<unsigned int> component(number_of_nodes);
	unsigned int number_of_components = connected_components(g, &component[0]);
		
	std::vector<unsigned int>::size_type i;
	
	vector <unsigned int> component_contents (number_of_components);
	
		for (i = 0; i != number_of_nodes; ++i)
		{
			component_contents[component[i]] = i;
		}
	
		for (i = 0; i < number_of_components -1 ; ++i)
		{
			t = vertex(component_contents[i],g);
			u = vertex(component_contents[i+1],g);
			add_edge (t, u, g);
		}
	
		assert (connected_components(g, &component[0]) == 1);
		
		//display_graph(g);
}



template <class Graph>
Graph * make_mrf_graph (const unsigned int rows, const unsigned int columns)
{
	unsigned int number_of_edges(0);
	
	typedef typename graph_traits<Graph>::vertex_descriptor VertexDescriptor;
	
	Graph * graph_ptr = new Graph(rows * columns);
	
	// Make "vertical" edges
	
	for (unsigned int i = 0 ; i < num_vertices(*graph_ptr); ++i)
	{
	
	if ( ((i+1) % rows) == 0 )
	{
		continue;
	}
		
	VertexDescriptor u = vertex(i, *graph_ptr);
	VertexDescriptor v = vertex(i+1, *graph_ptr);
	
	//cout << "edge from " << u+1 << " to " << v+1 << endl;
	
	add_edge(u,v,*graph_ptr);
	++number_of_edges;
	
	}
	
	// Make "horizontal" edges
	
	for (unsigned int i = 0 ; i < num_vertices(*graph_ptr); ++i)
	{
		
		if ( i >= rows*(columns-1) )
		{
			continue;
		}
		
		VertexDescriptor u = vertex(i, *graph_ptr);
		VertexDescriptor v = vertex(i+rows, *graph_ptr);
		
		add_edge(u,v,*graph_ptr);
		++number_of_edges;
		
		//cout << "edge from " << u+1 << " to " << v+1 << endl;
		
	}

	// cout << "We made " << number_of_edges << " edges and the graph has " << num_vertices(*graph_ptr) << " vertices." << endl;
	
	return graph_ptr;
}

template <class Graph>
Graph * make_simple_chain_graph (const unsigned int number_nodes)
{
	unsigned int number_of_edges(0);
	
	typedef typename graph_traits<Graph>::vertex_descriptor VertexDescriptor;
	
	Graph * graph_ptr = new Graph(number_nodes);
	
	// Make edges
	
	for (unsigned int i = 0 ; i < num_vertices(*graph_ptr) -1; ++i)
	{
		
		VertexDescriptor u = vertex(i, *graph_ptr);
		VertexDescriptor v = vertex(i+1, *graph_ptr);
		
		//cout << "edge from " << u+1 << " to " << v+1 << endl;
		
		add_edge(u,v,*graph_ptr);
		++number_of_edges;
	}
	
	cout << "We made " << number_of_edges << " edges and the graph has " << num_vertices(*graph_ptr) << " vertices." << endl;
	
	return graph_ptr;
}

// The following function uses hard coded constants for the number of possible variable values, etc...
/*
template <class Graph>
void randomize_discrete_pairwise_graph (Graph & g)
{
	
	typename graph_traits <Graph>::vertex_iterator u, u_end;
	
	unsigned int i = 1;
	
	//ConstantPotential * cp = new ConstantPotential ();
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		DiscreteRandomVariable * rv = new DiscreteRandomVariable(i,2);
		MRFMessageNode<DiscreteRandomVariable, DiscretePotential, boost::adjacency_list_traits<vecS, vecS, undirectedS>::vertex_descriptor> * node = new MRFMessageNode<DiscreteRandomVariable, DiscretePotential, boost::adjacency_list_traits<vecS, vecS, undirectedS>::vertex_descriptor>(rv, new SingleDiscretePotential (*rv));
		g[*u] = node;
		
		++i;
	}
	
	typename graph_traits <Graph>::edge_iterator e, e_end;
	
	vector <double> double_vec(4);
	
	for (tie (e, e_end) = edges (g); e != e_end; ++e)
	{
		g[*e] = new DiscretePotential();
		
		g[*e]->add_variable( *( g[source(*e,g)]->get_random_variable()));
		g[*e]->add_variable( *( g[target(*e,g)]->get_random_variable()));
		
		double_vec[0] = uniform_distribution_0_1();
		double_vec[1] = uniform_distribution_0_1();
		double_vec[2] = uniform_distribution_0_1();
		double_vec[3] = uniform_distribution_0_1();
		
		// Printing our values
		
		//cout << "Edge linking " << source(*e,g)+1 << " to " << target(*e,g)+1 << ": ";
		//cout << double_vec[0] << ", " << double_vec[1] << ", " << double_vec[2] << ", " << double_vec[3] << endl;
		
		g[*e]->setup_potential_values (double_vec);
	}
	
}
*/


template <class Graph>
void randomize_discrete_pairwise_symmetric_graph (Graph & g)
{
	typename graph_traits <Graph>::vertex_iterator u, u_end;
	typename graph_traits <Graph>::vertex_descriptor v;
	typename graph_traits <Graph>::edge_descriptor ed;
	
	unsigned int i = 1;
	
	const unsigned int rv_size = 2;
	const unsigned int obs_size = 2;
	
	//unsigned int number_vertices = num_vertices(g);
	
	typename graph_traits <Graph>::edge_iterator e, e_end;
	
	DiscreteRandomVariable * rv;
	MRFMessageNode<DiscreteRandomVariable, DiscretePotential, boost::adjacency_list_traits<vecS, vecS, undirectedS>::vertex_descriptor> * node;
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		rv = new DiscreteRandomVariable(i, rv_size );
		node = new MRFMessageNode<DiscreteRandomVariable, DiscretePotential, boost::adjacency_list_traits<vecS, vecS, undirectedS>::vertex_descriptor>(rv, new SingleDiscretePotential (*rv));
		
		g[*u] = node;
		++i;
	}
	
	vector <double> double_vec(rv_size * rv_size );
	
	for (tie (e, e_end) = edges (g); e != e_end; ++e)
	{
		g[*e] =  new DiscretePotential ();
		
		g[*e]->add_variable( *( g[source(*e,g)]->get_random_variable()));
		g[*e]->add_variable( *( g[target(*e,g)]->get_random_variable()));
		
		
		for ( unsigned int j = 0 ; j < rv_size ; ++j)
		{
			for ( unsigned int k = 0 ; k < rv_size ; ++k)
			{
				
				double_vec[j*rv_size+k] = 1.0; // We could randomize here as well
			}
		}
		
		for ( unsigned int j = 0 ; j < rv_size ; ++j)
		{
			double_vec[j*rv_size+j] = 10.0;
		}
		
		double_vec[0]= exp( 1.0 / 3.00);
		double_vec[1]= exp( -1.0 / 3.00);	
		double_vec[2]= exp( -1.0 / 3.00);
		double_vec[3]= exp( 1.0 / 3.00);
		
		g[*e]->setup_potential_values (double_vec);
		
	}
	
	vector <double> double_obs_vec(rv_size * obs_size );
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		
		rv = new DiscreteRandomVariable(i, obs_size ); // binary discrete RV
		
		node = new MRFMessageNode<DiscreteRandomVariable, DiscretePotential, boost::adjacency_list_traits<vecS, vecS, undirectedS>::vertex_descriptor> (rv, new ConstantPotential ());
		
		v = add_vertex(node, g);
		
		DiscretePotential * p = new DiscretePotential();
		
		p->add_variable( *( g[*u]->get_random_variable() ) );
		p->add_variable( *( g[v]->get_random_variable() ) ); // the second variable is the observation
		
		for ( unsigned int j = 0 ; j < obs_size ; ++j)
		{
			for ( unsigned int k = 0 ; k < rv_size ; ++k)
			{
				double_obs_vec[j*rv_size+k] = 1.0; // We could randomize here as well
			}
		}
		
		for ( unsigned int j = 0 ; j < obs_size ; ++j)
		{
			double_obs_vec[j*rv_size+j] = 5.0;
		}
		
		
		for ( unsigned int j = 0 ; j < obs_size ; ++j)
		{
			for ( unsigned int k = 0 ; k < rv_size ; ++k)
			{
				double_obs_vec[j*rv_size+k] = 1.0; // We could randomize here as well
			}
		}
		
		
		
		/*
		 
		 for (vector<double>::iterator it = double_obs_vec.begin(); it !=  double_obs_vec.end(); ++it)
		 {
			 cout << *it  << endl;
		 }
		 
		 */
		
		p->setup_potential_values (double_obs_vec);
		
		ed = add_edge(*u, v, g).first;
		
		g[ed] = p;
		
		++i;
	}
	
	//cout << "i (at end): " << i << endl;
}

/*
template <class Graph>
void fill_discrete_pairwise_symmetric_graph (Graph & g)
{
	typename graph_traits <Graph>::vertex_iterator u, u_end;
	typename graph_traits <Graph>::vertex_descriptor v;
	typename graph_traits <Graph>::edge_descriptor ed;
	
	unsigned int i = 1;
	
	const unsigned int rv_size = 10;
	const unsigned int obs_size = 10;
	
	//unsigned int number_vertices = num_vertices(g);