randomvariable.cc.svn-base

Probabilistic graphical models in matlab.

#include <iostream>
#include <fstream>

#include <RandomVariable.h>
#include <Generic.h>

using namespace std;
using namespace boost;

//****************** RandomVariable Implementation ********************* //

RandomVariable::RandomVariable (const unsigned int a ) : index(a), conditionned(false)
{
}

RandomVariable::~RandomVariable ()
{
}

//****************** End of RandomVariable Implementation ********************* //

//****************** ContinuousRandomVariable Implementation ********************* //

// Constructor

ContinuousRandomVariable::ContinuousRandomVariable(const unsigned int a): RandomVariable(a), last_sampled_value(0.0), number_distribution_values(100), particle_mean (false)
{
}

// Inherited functions

void ContinuousRandomVariable::obtain_inference_from_prior_samples () 
{
	if ( conditionned)
	{
		return;
	}
	
	double max_sample_value, min_sample_value;
	vector <double>::iterator max_sample, min_sample;
	
	unsigned int number_samples = prior_samples.size();
	
	// First find the max and mean
	
	max_sample = prior_samples.begin();
	max_sample_value = *max_sample;
	min_sample = prior_samples.begin();
	min_sample_value = *min_sample;
	
	for ( vector <double>::iterator it = prior_samples.begin(); it != prior_samples.end(); ++it)
	{
		if ( *it < min_sample_value)
		{
			min_sample = it;
			min_sample_value = *it;			
		}
		
		if ( *it > max_sample_value)
		{
			max_sample = it;
			max_sample_value = *it;			
		}
	}
	
	// Our max and mean is found
	
	// unsigned int number_steps = (unsigned int) number_samples / 10;
	
	unsigned int number_steps = number_distribution_values;
	
	double step = (max_sample_value-min_sample_value) / number_steps;
	
	// cout << "max: " << max_sample_value << " min: " << min_sample_value <<endl;
	
	// unsigned int number_steps = (unsigned int) ( (max_sample_value-min_sample_value) / step);
	
	// cout << "number_steps: " << number_steps << endl;
	
	distribution_approximation.clear ();
	distribution_approximation.resize (number_steps );
	
	//distribution_approximation.resize (10, pair <double, double> (0.0, 0.0) );
	
	/*
	 for (unsigned int i = 0; i < number_steps; ++i)
	 {
		 distribution_approximation.push_back( pair <double, double> (0.0, 0.0) );
	 }
	 */
	
	for ( unsigned int i = 0; i < distribution_approximation.size(); ++i)
	{
		distribution_approximation[i] = pair <double, double> (min_sample_value + step/2 + step * i , 0.0) ;
	}
	
	for ( vector <double>::iterator it = prior_samples.begin(); it != prior_samples.end(); ++it)
	{
		unsigned int a = (unsigned int) ( (*it  - min_sample_value ) / step);
		
		if ( a >= number_steps)
		{
			a = number_steps -1;
		}
		
		distribution_approximation[a].second += (double) (1.0 / number_samples);
	}
	
	compute_mean();
}

void ContinuousRandomVariable::display_probabilities () const
{
	/*
	 cout << "Our approximation for RV " << get_index() << " consists in " << prior_samples.size() << " samples: " << endl;
	 
	 for ( vector <double>::const_iterator it = prior_samples.begin(); it != prior_samples.end(); ++it)
	 {
		 cout << *it << endl;
	 }
	 
	 
	 cout << "Our approximation for RV " << get_index() << " consists in the following distribution:"  << endl;
	 
	 for ( vector <pair <double, double> >::const_iterator it = distribution_approximation.begin(); it != distribution_approximation.end(); ++it)
	 {
		 cout << it->first << " : " << it->second << endl;
	 }
	 */
	cout << "Mean for RV " << get_index() << " is: "  << mean << endl;
}

void ContinuousRandomVariable::initialize()
{
	conditionned = false;
	last_sampled_value = 0.0;
	prior_samples.clear();
	distribution_approximation.clear();
	particle_mean = false;
	particle_means_vec.clear();
}

void ContinuousRandomVariable::export_inference_to_matlab () const
{
	ofstream output_file ("export.mat");
	
	cout << "Exporting RV " << get_index() << " to Matlab." << endl;
	
	for ( vector <pair <double, double> >::const_iterator it = distribution_approximation.begin(); it != distribution_approximation.end(); ++it)
	{
		output_file << it->first << " " << it->second << "\n";
	}
	
	output_file << endl;	
}

// Normal functions

void ContinuousRandomVariable::record_sample (const double a) 
{
	prior_samples.push_back(a);
	
	//last_sampled_value = a;
}

void ContinuousRandomVariable::compute_mean () 
{
	if (particle_mean)
	{
		double a(0.0);
		
		for ( unsigned int i = 0; i <  particle_means_vec.size();++i)
		{
			a += particle_means_vec[i];
		}
		
		mean = a / particle_means_vec.size();
		
		return;
	}
	
	double sum(0.0);
	
	for ( vector <double>::iterator it = prior_samples.begin(); it != prior_samples.end(); ++it)
	{
		sum += *it;
		
		/*
		if (index == 20)
		{
		
		//cout << "Sample: " << *it << endl;
		}
		*/
		}
	
	
	
	mean = sum /  prior_samples.size();
	/*
	if (index == 20)
	{
		
		//cout << "Mean: " << mean << endl;
		//cout << "ref: " << reference_mean << endl;
		//cout << "Mean - ref: " << mean -reference_mean << endl;
	}
	*/

}


double ContinuousRandomVariable::sample_from_particles (const vector <double> & positions, const vector <double> & weights)
{
	assert ( positions.size() == weights.size());
	
	unsigned int i;
	
	double a = 0;
	double random_number = uniform_distribution_0_1();
	
	for (i = 0; i < weights.size(); ++i )
	{
		a+= weights[i];
		
		if (random_number < a)
		{
			break;
		}
		
	}
	
	assert ( i < weights.size());
	
	prior_samples.push_back(positions[i]);
	
	last_sampled_value = positions[i];
	
	return positions[i];
}

// Testing purposes

void  ContinuousRandomVariable::set_mean_bw_particles(vector <double> & particle_vec,  vector <double> & weights_vec)
{
double a(0.0);

for ( unsigned int i = 0; i <  particle_vec.size();++i)
{
	a += particle_vec[i] * weights_vec[i];
}

//mean = a / particle_vec.size();

particle_mean = true;

particle_means_vec.push_back(a);

}

#ifndef NDEBUG

void ContinuousRandomVariable::debug_display_state()
{
	cout << "Number of samples: " << prior_samples.size() << endl;
	
}

#endif

//****************** End of ContinuousRandomVariable Implementation ********************* //

//****************** DiscreteRandomVariable Implementation ********************* //

// Constructors, destructor

DiscreteRandomVariable::DiscreteRandomVariable(unsigned int a): RandomVariable(a)
{
}

DiscreteRandomVariable::DiscreteRandomVariable(unsigned int a, unsigned int b) : RandomVariable(a), start_loop (true), end_loop (false), number_values(b), last_sampled_value(0), sampling_probabilities(b, 1.0 / b), prior_samples(b,0), particle_mean(false)
{
}

DiscreteRandomVariable::~DiscreteRandomVariable()
{	
}

// Inherited virtual methods

void DiscreteRandomVariable::initialize()
{
	cout << "Performing INIT" << endl;
	
	//mean = 0.0;
	//conditionned = false;
	last_sampled_value = 0;
	prior_samples.clear();
	prior_samples.resize(number_values, 0);
	sampling_probabilities.clear();
	sampling_probabilities.resize(number_values, 1.0 / number_values) ; // not really needed
	particle_mean = false;
	particle_means_vec.clear();
}


void DiscreteRandomVariable::remove_prior_samples()
{
	prior_samples.clear();
	prior_samples.resize(number_values, 0);
}


void DiscreteRandomVariable::obtain_inference_from_prior_samples()
{
	if (  conditionned)
	{
		return;
	}
	
	assert ( prior_samples.size() == sampling_probabilities.size() );
	
	vector <unsigned int>::iterator current_index = prior_samples.begin();
	
	for (vector <double>::iterator it = sampling_probabilities.begin(); it != sampling_probabilities.end(); ++it)
	{
		*it = *current_index;
		++current_index;
	}
			 
	normalize_probabilities ();
}

void DiscreteRandomVariable::display_probabilities() const
{
	if ( conditionned)
	{
		cout << "Variable " << get_index() << " is observed and set to " << last_sampled_value << endl;
		return;
	}
	
	cout << "This is RV " << get_index() << ".\n" ;
	
	for (unsigned int i =0; i < number_values; ++i)
	{
		cout << "Value " <<  i << ": " << sampling_probabilities[i] << endl;
	}	
}

void DiscreteRandomVariable::export_inference_to_matlab () const
{
	cout << "Feature not implemented yet (export_inference_to_matlab () )." << endl;
}

// Normal methods

double DiscreteRandomVariable::compute_error_to_reference() const
{
	vector<double>::const_iterator jt = sampling_probabilities.begin();
	
	double error(0.0);
	
	for (vector<double>::const_iterator it = reference_marginals.begin(); it != reference_marginals.end(); ++it)
	{
		error += (*it - *jt) * (*it - *jt);		
		++jt;
	}
	
	return error;
}

void  DiscreteRandomVariable::set_reference_marginals (const vector <double> & rm_vec)
{
	assert (number_values == rm_vec.size());
	
	reference_marginals.resize(number_values);
	
	copy(rm_vec.begin(), rm_vec.end(), reference_marginals.begin());
}

void DiscreteRandomVariable::normalize_probabilities ()
{	
	/*
	if (conditionned)
	{
		cout << "Sampling probs" << endl;
		
		for (vector <double>::iterator it = sampling_probabilities.begin(); it != sampling_probabilities.end(); ++it)
		{
			cout << *it << endl;
		}
	
		cout << "Sampling probs - end" << endl;
	}
	 */
	
	 if (!fast_double_normalize_over (sampling_probabilities.begin(),sampling_probabilities.end()) )
	{
		 cout << "ERROR [FATAL]: while normalizing over Discrete RV " << index << " sampling probabilities." << endl;
		 exit (0);
	}
}

void DiscreteRandomVariable::observe_variable_value(const unsigned int a)
{
	conditionned = true;
	//conditionned_value = a;
	last_sampled_value = a;
	
	for (vector <double>::iterator it = sampling_probabilities.begin(); it != sampling_probabilities.end(); ++it)
	{
		*it = 0.0;
	}
	
	sampling_probabilities[a] = 1.0;