old_code.svn-base

Probabilistic graphical models in matlab.

		
		// Computing new weights

		/* not needed anymore AND WRONG
			
		for (vector <double>::iterator current_particle_weight = current_weights.begin(); current_particle_weight != current_weights.end(); ++current_particle_weight)
		{		
			observation_potential->set_variable_value (current_variable_index, *current_particle_position);
			
			interaction_potential->set_variable_value (previous_variable_index, *previous_particle_position);
			
			interaction_potential->set_variable_value (current_variable_index, *current_particle_position);
			
			*current_particle_weight = *previous_particle_weight * observation_potential->get_potential_value() * interaction_potential->get_potential_value() / proposal->get_proposal_weight(*current_particle_position, *previous_particle_position, 1.0 ); // default variance is used
			
			++current_particle_position;
			++previous_particle_position;
			++previous_particle_weight;
		}
		
		*/
		
		
		/*
#include <ParticleFilters.h>
#include <Potential.h>
#include <Node.h>
#include <RandomVariable.h>
#include <GraphTreePartition.h>
#include <OutputGraph.h>

#include <boost/graph/subgraph.hpp>
#include <boost/graph/copy.hpp>

// Explicit Template Instantiations //

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

//template void non_parametric_tree_gibbs_sampler < NPGraph> (NPGraph &,  unsigned int, unsigned int);

// End of Explicit Template instantiations //


//*************** DiscreteParticleFilterProposal Declaration ******************//

class DiscreteParticleFilterProposal
{
	
private:
	
	vector <double> current_sampling_probabilities;
	vector <double> general_sampling_probabilities;

	
public:
		
	DiscreteParticleFilterProposal();
	
	void setup_discrete_proposal( const ChainedSingleDiscretePotential &, const DiscreteRandomVariable  &);
	void add_discrete_proposal( DiscretePotential &, const DiscreteRandomVariable  &);
	
	unsigned int get_proposal_position(DiscreteRandomVariable &);
	
	double get_proposal_weight(const unsigned int ) const;
	
};

//*************** End of DiscreteParticleFilterProposal Declaration ******************//

//*************** End of DiscreteParticleFilterProposal Implementation ******************//

void DiscreteParticleFilterProposal::setup_discrete_proposal( const ChainedSingleDiscretePotential & csdp, const DiscreteRandomVariable  & drv)
{
	general_sampling_probabilities.clear();
	general_sampling_probabilities.resize( drv.get_number_values() , 1.0);
	current_sampling_probabilities.clear();
	current_sampling_probabilities.resize( drv.get_number_values() , 1.0);
	
	unsigned int variable_index = drv.get_index();
	
		for (unsigned int sp = 0 ; sp < general_sampling_probabilities.size(); ++sp)
		{
			for (list <DiscretePotential * >::const_iterator it = csdp.potential_lst.begin(); it != csdp.potential_lst.end(); ++it)
			{
				general_sampling_probabilities[sp] = general_sampling_probabilities[sp] * (*it)->obtain_potential_value ( variable_index, sp);
			}
		}	
}

void DiscreteParticleFilterProposal::add_discrete_proposal(DiscretePotential & dp, const DiscreteRandomVariable  & drv)
{
	unsigned int variable_index = drv.get_index();
	
	for (unsigned int sp = 0 ; sp < current_sampling_probabilities.size(); ++sp)
	{
		current_sampling_probabilities[sp] = general_sampling_probabilities[sp] * dp.obtain_potential_value ( variable_index, sp);
	}	
}

unsigned int DiscreteParticleFilterProposal::get_proposal_position(DiscreteRandomVariable & rv)
{	
	double (*pf) (double) = &ReturnSelf <double>;
	void (*pf2) (double &,double) = &NaturalSetDouble <double>;
	
	normalize_over (current_sampling_probabilities.begin(), current_sampling_probabilities.end(), pf, pf2);

	return rv.sample_without_recording(current_sampling_probabilities);
}

inline double DiscreteParticleFilterProposal::get_proposal_weight(const unsigned int a) const
{
	return current_sampling_probabilities[a];
	
}

//*************** End of DiscreteParticleFilterProposal Implementation ******************//

typedef property< edge_index_t, unsigned int, DiscretePotential *> DiscreteEdgeProperty;

typedef adjacency_list<vecS, vecS, undirectedS, MessageNode<DiscreteRandomVariable, boost::adjacency_list_traits<vecS, vecS, undirectedS>::vertex_descriptor> *, DiscreteEdgeProperty > DiscreteGraph;

template <>
void non_parametric_tree_gibbs_sampler < DiscreteGraph> (DiscreteGraph & g,  const unsigned int gibbs_steps, const unsigned int number_particles)
{	
	//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()->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);
	
	for (tie(u,u_end) = vertices (g); u!= u_end; ++u)
	{		
		delete g[*u]->get_potential();
		
		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);

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


template < class Graph,  class VertexDescriptor, class PositionMap , class WeightMap>
VertexDescriptor forward_filtering_implementation (Graph & g, const VertexDescriptor root_node, PositionMap positions, WeightMap weights)
{
	// Variable declarations
	
	DiscreteRandomVariable * current_random_variable;
	
	unsigned int current_variable_index, previous_variable_index ;
	
	VertexDescriptor double_previous_node, previous_node, current_node;
	
	previous_node = root_node;
	
	double_previous_node = root_node;
	
	double (*pf) (double) = &ReturnSelf <double>;
	void (*pf2) (double &,double) = &NaturalSetDouble <double>;
	
	typename graph_traits<Graph>::adjacency_iterator u, u_end;
	
	ParticleFilterProposal * proposal = new ParticleFilterProposal();
	
	unsigned int number_particles = get (positions, root_node).size();
	
	initialize_particles (g, get (positions, root_node), get (weights, root_node), root_node, proposal);
	
	DiscretePotential * internal_potential;
	DiscretePotential * observation_potential;
	DiscretePotential * interaction_potential;
	
	vector <double> cumulative_distribution (number_particles);
	vector <double> stratified_weights (number_particles);
	vector <unsigned int> original_particles (number_particles);
	
	// We got our initial distribution covered, resample it ?? Currently, NO resampling is performed.
	
	// Initialization done.
	// Beginning of main loop.
	
	while (true)
	{
		
		// * We must obtain:
		
		// * potential coming from the edges linking the node to the other nodes not in the chain, call that internal_potential (because we arranged it as an internal potential)
		// * potential between node and observation, call that observation_potential;
		// * potential between the two nodes, call that interaction potential;
		
		for (tie (u, u_end) = adjacent_vertices(previous_node, g); u != u_end; ++u)
		{
			if (out_degree(*u,g) != 1)
			{
				
				if ( *u != double_previous_node)
				{
					current_node = *u;
					interaction_potential = g[ edge(previous_node, current_node, g).first];
					
					internal_potential = g[current_node]->get_potential();
				}
			}
		}
		
		for (tie (u, u_end) = adjacent_vertices(current_node, g); u != u_end; ++u)
		{
			if (out_degree(*u,g) == 1)
			{
				observation_potential =  g[edge(current_node, *u, g).first];
			}
		}
		
		//cout << "Current node is " << current_node << " (" << g[current_node]->get_random_variable()->get_index() << ")" << ", previous was " << previous_node << " (" << g[previous_node]->get_random_variable()->get_index() << ")" << endl;
		
		//cout << "Current node (FF) is " << g[current_node]->get_random_variable()->get_index()  << endl;
		
		previous_variable_index = g[previous_node]->get_random_variable()->get_index();
		
		current_random_variable = g[current_node]->get_random_variable();
		current_variable_index = current_random_variable->get_index();
		
		// Sampling particle from proposal
		
		vector <double> & current_weights = get (weights, current_node);		
		vector <double> & previous_weights = get (weights, previous_node);
		vector <unsigned int> & current_positions = get (positions, current_node);		
		vector <unsigned int> & previous_positions = get (positions, previous_node);
		
		vector <unsigned int>::iterator previous_particle_position = previous_positions.begin();
		vector <double>::iterator previous_particle_weight = previous_weights.begin();
		vector <double>::iterator current_particle_weight = current_weights.begin();
		
		//cout << "First weight (before update): " << *current_weights.begin() << endl;
		
		proposal->setup_discrete_proposal(static_cast<ChainedSingleDiscretePotential *>(internal_potential), current_variable_index);
		
		for (vector <double>::iterator current_particle_position = current_positions.begin(); current_particle_position != current_positions.end(); ++current_particle_position)
		{
			// Sampling the new particle
			
			proposal->add_discrete_proposal(*previous_particle_position, *interaction_potential);
			
			*current_particle_position = proposal->get_proposal_position(*current_random_variable); // This also sets the current_random_variable.last_sampled_value
			
			// Compute the new weight
			
			observation_potential->set_variable_value (*current_random_variable); // All of these are discrete potentials	
			
			interaction_potential->set_variable_value (*current_random_variable);
			interaction_potential->set_variable_value_by_index (previous_variable_index, *previous_particle_position);
			
			// Setting Internal Potential
			
			internal_potential->set_variable_value (*current_random_variable);
			
			// cout << "New particle, weight of old particle, obs, int, internal, prop: " << *current_particle_position << ", " << *previous_particle_weight << ", " << observation_potential->get_potential_value() << ", " << interaction_potential->get_potential_value() << ", " <<internal_potential->get_potential_value() << ", " << proposal->get_proposal_weight(*current_particle_position) << endl;
			
			*current_particle_weight = ( *previous_particle_weight * observation_potential->get_potential_value() * interaction_potential->get_potential_value() * internal_potential->get_potential_value() )/ proposal->get_proposal_weight(*current_particle_position);
			
			// cout << "Weight of new particle: " << *current_particle_weight << endl;
			
			++previous_particle_position;
			++previous_particle_weight;
			++current_particle_weight;
		}
		
		// Normalize the weights 
		
		normalize_over (current_weights.begin(), current_weights.end(), pf, pf2);
		
		// Resample part
		
		resample ( current_positions, current_weights, cumulative_distribution, stratified_weights, original_particles);
		if ( out_degree(current_node,g) == 2)
		{
			// We are finished

			return current_node;
		}
		
		double_previous_node = previous_node;
		previous_node = current_node;
		
	}
	
	delete proposal;
	
}
*/

	
	cout << endl << "*******************" << endl;
	
	//partition_graph_into_trees(*g);
	
	//create_factor_graph_test();
	
	//create_mrf_graph_test();
	
	typedef graph_traits <MRFGraph>::vertex_descriptor VertexDescriptor;
	typedef graph_traits <MRFGraph>::edge_iterator EdgeIterator;
	
	/*
	list < unsigned int > l;
	
	l.push_back(1);
	l.push_back(2);
	l.push_back(3);
	l.push_back(4);
	l.push_back(5);

	combine_nodes(*g, l, 31);
	
	l.clear();
	l.push_back(6);
	l.push_back(7);
	l.push_back(8);
	l.push_back(9);
	l.push_back(10);
	
	combine_nodes(*g, l, 32);
	
	l.clear();
	
	l.push_back(11);
	l.push_back(12);
	l.push_back(13);
	l.push_back(14);
	l.push_back(15);
	
	combine_nodes(*g, l, 33);
	
	l.clear();
	
	l.push_back(16);
	l.push_back(17);
	l.push_back(18);
	l.push_back(19);
	l.push_back(20);
	
	combine_nodes(*g, l, 34);
	
	l.clear();
	
	l.push_back(21);
	l.push_back(22);
	l.push_back(23);
	l.push_back(24);