inputfunctions.cc.svn-base

Probabilistic graphical models in matlab.

	
	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 );
	
	variate_generator <lagged_fibonacci607 &, normal_distribution < > > vg  (fibonnacci_number_generator, normal_distribution < > () );
	
	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
			}
		}

		//double strange = abs( vg() );
		
		double phys = exp( abs( vg() ) );
		
		//cout << phys << " " << strange << endl;
		
		for ( unsigned int j = 0 ; j < rv_size ; ++j)
		{
			double_vec[j*rv_size+j] = phys;
		}
		
		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;
		}
		
		
		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 randomize_non_parametric_pairwise_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;
	
	//unsigned int number_vertices = num_vertices(g);
	
	typename graph_traits <Graph>::edge_iterator e, e_end;
	
	RandomVariable * rv;
	GraphicalModelNode * node;
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		rv = new ContinuousRandomVariable(i);
		node = new GraphicalModelNode /* <RandomVariable, Potential, boost::adjacency_list_traits<vecS, vecS, undirectedS>::vertex_descriptor> */ (rv, new ConstantPotential () );
		g[*u] = node;
				
		++i;
	}
	
	for (tie (e, e_end) = edges (g); e != e_end; ++e)
	{
		g[*e] =  new GaussianPotential ( 0.1, *( g[source(*e,g)]->get_random_variable()) , *( g[target(*e,g)]->get_random_variable()) ); // first number is variance !!
	}
	
	//cout << "i: " << i << endl;
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		
		rv = new DiscreteRandomVariable(i, 2); // binary discrete RV
		
		node = new GraphicalModelNode (rv, new ConstantPotential ());
		
		v = add_vertex(node, g);
		
		ContinuousPotential * p = new BernouilliPotential( *(g[v]->get_random_variable() ), *(g[*u]->get_random_variable() ) );
		
		ed = add_edge(*u, v, g).first;

		g[ed] = p;
		
		++i;
	}
	
	//cout << "i (at end): " << i << endl;
}

template < class Graph>
void initialize_graph(Graph & g)
{
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		g[*u]->get_random_variable()->initialize();
	}
}

// Should be an "output" function...


template < class Graph>
void export_computed_means(string output_file_name, Graph & g)
{

	typename graph_traits<Graph>::vertex_iterator u, u_end;
	
	ofstream output_file (output_file_name.c_str());
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		if (out_degree (*u, g) != 1)
		{

			output_file << static_cast <ContinuousRandomVariable * > (g[*u]->get_random_variable())->get_index() << " " << static_cast <ContinuousRandomVariable * > (g[*u]->get_random_variable())->get_mean() << "\n";
			
		}
	}
	
	output_file << endl;

}

template <>
void export_computed_means <DiscretePairwiseGraph> (string output_file_name, DiscretePairwiseGraph & g)
{
	
	graph_traits<DiscretePairwiseGraph>::vertex_iterator u, u_end;
	
	ofstream output_file (output_file_name.c_str());
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		if (out_degree (*u, g) != 1)
		{
			
			//static_cast <DiscreteRandomVariable * > (g[*u]->get_random_variable())->compute_true_mean();
			
			static_cast <DiscreteRandomVariable * > (g[*u]->get_random_variable())->compute_mean();
			
			output_file << static_cast < DiscreteRandomVariable * > (g[*u]->get_random_variable())->get_index() << " " << static_cast <DiscreteRandomVariable * > (g[*u]->get_random_variable())->get_mean() << "\n";
		}
	}
	
	output_file << endl;
	
}

template < class Graph>
void read_reference_means (string input_file_name, Graph & g)
{
	ifstream input_file (input_file_name.c_str());
	
	if (!input_file.good())
	{
		cout << "Failed to open the file '" << input_file_name << "'" << endl;
		exit (0);
	}
	
	unsigned int variable_index(0);
	double reference_mean(0.0);
	
	while (input_file >> variable_index)
	{
		//input_file >> variable_index;
		input_file >> reference_mean;
		
		assert (g[vertex(variable_index-1,g)]->get_random_variable()->get_index() == variable_index );
		
		//cout << "With variable " << variable_index << " we associated mean " <<  reference_mean << endl;
		
		static_cast <ContinuousRandomVariable * > (g[vertex(variable_index-1,g)]->get_random_variable())->set_reference_mean (reference_mean);

	}
}

template <>
void read_reference_means <DiscretePairwiseGraph> (string input_file_name, DiscretePairwiseGraph & g)
{
	ifstream input_file (input_file_name.c_str());
	
	if (!input_file.good())
	{
		cout << "Failed to open the file '" << input_file_name << "'" << endl;
		exit (0);
	}
	
	unsigned int variable_index(0);
	double reference_mean(0.0);
	
	while (input_file >> variable_index)
	{
		//input_file >> variable_index;
		input_file >> reference_mean;
		
		assert (g[vertex(variable_index-1,g)]->get_random_variable()->get_index() == variable_index );
		
		//cout << "With variable " << variable_index << " we associated mean " <<  reference_mean << endl;
		
		static_cast <DiscreteRandomVariable * > (g[vertex(variable_index-1,g)]->get_random_variable())->set_reference_mean (reference_mean);
		
	}
}

template <class Graph>
void initialize_random_variables (Graph & g)
{
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	typename graph_traits<Graph>::adjacency_iterator v, v_end;
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		
		if (out_degree(*u,g) != 1)
		{
			
			g[*u]->get_random_variable()->initialize();
			
			for ( tie (v, v_end) = adjacent_vertices(*u, g); v != v_end; ++v)
			{
				g[edge(*u,*v,g).first]->set_variable_value( *(g[*u]->get_random_variable()));
			}
			
		}
		
	}
}

template <class Graph>
void set_observations (Graph & g)
{
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	typename graph_traits<Graph>::adjacency_iterator v, v_end;
	
	// We should get that in a better way
	
	unsigned int obs_size = 25;
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		
		if (out_degree(*u,g) == 1)
		{
			obs_size = static_cast <DiscreteRandomVariable *> (g[*u]->get_random_variable())->get_number_values();
			break;
		}
		
	}
	
	fibonnacci_number_generator.seed(15);
	
	variate_generator <lagged_fibonacci607 &, uniform_smallint < > > vg (fibonnacci_number_generator, uniform_smallint <  > (0,obs_size -1) );

	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		
		if (out_degree(*u,g) == 1)
		{
			unsigned int a = vg();
			
			//a = 0;
			//cout << "For variable " << g[*u]->get_random_variable()->get_index() << ": " << a << endl;
			
			tie (v, v_end) = adjacent_vertices(*u, g);
			
			static_cast < DiscreteRandomVariable * > (g[*u]->get_random_variable())->observe_variable_value(a);
			
			g[edge(*u,*v,g).first]->set_variable_value( *(g[*u]->get_random_variable()));
		}
		
	}
}


template <class Graph>
void set_special_mrf_observations (Graph & g)
{
	typename graph_traits<Graph>::vertex_iterator u, u_end;
	typename graph_traits<Graph>::adjacency_iterator v, v_end;
	
	fibonnacci_number_generator.seed(15);
	
	variate_generator <lagged_fibonacci607 &, uniform_smallint < > > vg (fibonnacci_number_generator, uniform_smallint <  > (0,1) );
	
	for (tie (u, u_end) = vertices (g); u != u_end; ++u)
	{
		
		if (out_degree(*u,g) == 1)
		{
			
			//cout << "For variable " << g[*u]->get_random_variable()->get_index() << ": " << a << endl;
			
			tie (v, v_end) = adjacent_vertices(*u, g);
			
			if (7 == (g[*u]->get_random_variable()->get_index() % mrf_rows) )
			{
			
				static_cast < DiscreteRandomVariable * > (g[*u]->get_random_variable())->observe_variable_value(0);
				
				dynamic_cast< BernouilliPotential * > (g[edge(*u,*v,g).first])->set_threshold (4.0);
			}
			
			else
				
			{
				static_cast < DiscreteRandomVariable * > (g[*u]->get_random_variable())->observe_variable_value(1);
			}
				
			g[edge(*u,*v,g).first]->set_variable_value( *(g[*u]->get_random_variable()));
			
			
		}
		
	}
}

// New functions

template <class Graph>
void create_pairwise_square_lattice (const unsigned int rows, const unsigned int columns, Graph & g)
{
	unsigned int number_of_edges(0);
	
	typedef typename graph_traits<Graph>::vertex_descriptor VertexDescriptor;

	// Make "vertical" edges
	
	for (unsigned int i = 0 ; i < num_vertices(g); ++i)
	{		
		if ( ((i+1) % rows) == 0 )
		{
			continue;
		}
		
		VertexDescriptor u = vertex(i, g);
		VertexDescriptor v = vertex(i+1, g);
		
		add_edge(u,v,g);
		++number_of_edges;
		
	}
	
	// Make "horizontal" edges
	
	for (unsigned int i = 0 ; i < num_vertices(g); ++i)
	{
		if ( i >= rows*(columns-1) )
		{
			continue;
		}
		
		VertexDescriptor u = vertex(i, g);
		VertexDescriptor v = vertex(i+rows, g);
		
		add_edge(u,v,g);
		++number_of_edges;
	}
}