examples/grid/grid2Dvineyard.hpp
author Christos Mantoulidis <cmad@stanford.edu>
Tue, 04 Aug 2009 13:23:16 -0700
branchdev
changeset 156 f75fb57d2831
parent 30 6d4e450015e4
permissions -rw-r--r--
Changed implementation of WeightedRips to store simplex values (max distance between simplices' vertices) as an invisible layer on top of each simplex object, so that the data() field of WeightedRips has been freed for use by the users again.

#include <utilities/log.h>

/* Implementation */
	
Grid2DVineyard::
Grid2DVineyard(Grid2D* g): vertices_(g->size())
{
	grid_stack_.push_back(g); 
	for (CoordinateIndex i = 0; i < g->size(); ++i)
		vertices_[i].set_index(i);

	evaluator_ = new StaticEvaluator(grid(), 0);
	vineyard_ = new Vineyard(evaluator_);

	filtration_ = new LSFiltration(vertices_.begin(), vertices_.end(), VertexComparison(grid()), vineyard_);
	add_simplices();
}

Grid2DVineyard::
~Grid2DVineyard()
{
	delete filtration_;
	delete vineyard_;
	delete evaluator_;
}

void
Grid2DVineyard::
compute_pairing()
{
	filtration_->fill_simplex_index_map();
	filtration_->pair_simplices(filtration_->begin(), filtration_->end());
	vineyard_->start_vines(filtration_->begin(), filtration_->end());
}

void					
Grid2DVineyard::
compute_vineyard(Grid2D* g, bool explicit_events)
{
	AssertMsg(filtration_->is_paired(), "Simplices must be paired for a vineyard to be computed");
	
	typedef Traits::Kinetic_kernel::Point_1 								Point;
	typedef Traits::Kinetic_kernel::Function_kernel::Construct_function 	CF; 
	typedef Traits::Kinetic_kernel::Motion_function 						F; 
	
	Traits tr(0,1);
	Simulator::Handle sp = tr.simulator_handle();
	ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
	Sort sort(tr, SortVisitor(this));
	
	// Setup the (linear) trajectories
	std::cout << "Setting up trajectories" << std::endl;
	CF cf; 
	kinetic_map_.clear();
	for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
	{
		ValueType val0 = (*grid())(cur->index());
		ValueType val1 = (*g)(cur->index());
		F x = cf(F::NT(val0), F::NT(val1 - val0));			// x = val0 + (val1 - val0)*t
		Point p(x);
		cur->set_kinetic_key(apt->insert(p));
		kinetic_map_[cur->kinetic_key()] = cur->get_order();
		if (cur->index() % 10000 == 0)
			std::cout << "Added trajectory: " << cur->index() << " " << val0 << " " << val1 << std::endl;
	}
	
	// Process all the events (compute the vineyard in the process)
	change_evaluator(new KineticEvaluator(sp, apt, num_grids() - 1));
	if (explicit_events)
	{
		while (sp->next_event_time() < 1)
		{
			std::cout << "Next event time: " << sp->next_event_time() << std::endl;
			sp->set_current_event_number(sp->current_event_number() + 1);
			std::cout << "Processed event" << std::endl;
		}
	} else
		sp->set_current_time(1.0);
	std::cout << "Processed " << sp->current_event_number() << " events" << std::endl;
	
	// Add the grid to the stack
	grid_stack_.push_back(g); 
	change_evaluator(new StaticEvaluator(grid(), num_grids() - 1));
	vineyard_->record_diagram(filtration_->begin(), filtration_->end());
}
		
void 					
Grid2DVineyard::
swap(Key a, Key b)
{
	VertexOrderIndex ao = kinetic_map_[a], bo = kinetic_map_[b];
	AssertMsg(filtration_->get_vertex_cmp()(ao, bo), "In swap(a,b), a must precede b");
	filtration_->transpose_vertices(ao);
	AssertMsg(filtration_->get_vertex_cmp()(bo, ao), "In swap(a,b), b must precede a after the transposition");
}

void 
Grid2DVineyard::
add_simplices()
{
	// Takes advantage of LowerStarFiltration's smart append (which allows faces
	// to be inserted after cofaces, since everything is rearranged in the
	// proper lower star order anyway). Also note that vertices were added by
	// LowerStarFiltration's constructor
	for (CoordinateIndex x = 0; x < grid()->xsize() - 1; ++x)
		for (CoordinateIndex y = 0; y < grid()->ysize() - 1; ++y)
		{
			VertexIndex v(&vertices_[grid()->seq(x,y)]);
			VertexIndex vh(&vertices_[grid()->seq(x+1,y)]);
			VertexIndex vv(&vertices_[grid()->seq(x,y+1)]);
			VertexIndex vd(&vertices_[grid()->seq(x+1,y+1)]);

			Simplex sh(2, v);
			sh.add(vh);	filtration_->append(sh);		// Horizontal edge
			sh.add(vd);	filtration_->append(sh);		// "Horizontal" triangle
			
			Simplex sv(2, v);
			sv.add(vv);	filtration_->append(sv);		// Vertical edge
			sv.add(vd);	filtration_->append(sv);		// "Vertical" triangle
			
			Simplex sd(2, v);
			sd.add(vd); filtration_->append(sd);		// Diagonal edge

			if (y == grid()->ysize() - 2)
			{
				Simplex s(1, vv); 
				s.add(vd); filtration_->append(s);		// Top edge
			}
			if (x == grid()->xsize() - 2)
			{
				Simplex s(1, vh); 
				s.add(vd); filtration_->append(s);		// Right edge
			}
		}
}

void
Grid2DVineyard::
change_evaluator(Evaluator* eval)
{
	AssertMsg(evaluator_ != 0, "change_evaluator() assumes that existing evaluator is not null");
		
	delete evaluator_;
	evaluator_ = eval;
	vineyard_->set_evaluator(evaluator_);
}