include/topology/vineyard.hpp
author Dmitriy Morozov <dmitriy@mrzv.org>
Sat, 11 Apr 2009 10:29:53 -0700
branchdev
changeset 126 3c3e77ac43d2
parent 87 2c2e2f3b5d15
child 179 d15c6d144645
permissions -rw-r--r--
Added Python bindings for the Rips complex + an example

/* Implementations */

#include <fstream>
#include <sstream>

#include "utilities/log.h"

#ifdef LOGGING
static rlog::RLogChannel* rlVineyard =			DEF_CHANNEL("topology/vineyard", rlog::Log_Debug);
#endif // LOGGING

template<class FS>
void
Vineyard<FS>::
start_vines(Index bg, Index end)
{
	AssertMsg(evaluator != 0, "Cannot start vines with a null evaluator");
	for (Index cur = bg; cur != end; ++cur)
	{
		if (!cur->sign()) continue;
		Dimension dim = cur->dimension();
		
		if (dim >= vines.size())
		{
			AssertMsg(dim == vines.size(), "New dimension has to be contiguous");
			vines.push_back(VineList());
            vines_vector.push_back(boost::prior(vines.end()));
		}

		start_vine(cur);
		record_knee(cur);
	}
}

template<class FS>
void
Vineyard<FS>::
switched(Index i, Index j)
{
	VineS* i_vine = i->vine();
	VineS* j_vine = j->vine();
	i->set_vine(j_vine);
	j->set_vine(i_vine);

	// Since the pairing has already been updated, the following assertions should be true
	AssertMsg(i->vine() == i->pair()->vine(), "i's vine must match the vine of its pair");
	AssertMsg(j->vine() == j->pair()->vine(), "j's vine must match the vine of its pair");

	if (!i->sign()) i = i->pair();
	if (!j->sign()) j = j->pair();
	record_knee(i);
	record_knee(j);
}

template<class FS>
void
Vineyard<FS>::
start_vine(Index i)
{
	rLog(rlVineyard, "Starting new vine");
	AssertMsg(i->sign(), "Can only start vines for positive simplices");
		
	Dimension dim = i->dimension();
	vines_vector[dim]->push_back(VineS());
	i->set_vine(&vines_vector[dim]->back());
	i->pair()->set_vine(i->vine());
}
	
/// Records the current diagram in the vineyard
template<class FS>
void 
Vineyard<FS>::
record_diagram(Index bg, Index end)
{
	rLog(rlVineyard, "Entered: record_diagram()");
	AssertMsg(evaluator != 0, "Cannot record diagram with a null evaluator");
	
	for (Index i = bg; i != end; ++i)
	{
		AssertMsg(i->vine() != 0, "Cannot process a null vine in record_diagram");
		if (!i->sign())		continue;
		record_knee(i);
	}
}


template<class FS>
void			
Vineyard<FS>::
save_edges(const std::string& filename) const
{
	for (unsigned int i = 0; i < vines_vector.size(); ++i)
	{
		std::ostringstream os; os << i;
		std::string fn = filename + os.str() + ".edg";
		std::ofstream out(fn.c_str());
		for (typename VineList::const_iterator vi = vines_vector[i]->begin(); vi != vines_vector[i]->end(); ++vi)
			for (typename VineS::const_iterator ki = vi->begin(), kiprev = ki++; ki != vi->end(); kiprev = ki++)
			{
				if (kiprev->is_infinite() || ki->is_infinite()) continue;
				out << kiprev->birth << ' ' << kiprev->death << ' ' << kiprev->time << std::endl;
				out << ki->birth << ' ' << ki->death << ' ' << ki->time << std::endl;
			}
		out.close();
	}
}

/// Records a knee for the given simplex
template<class FS>
void
Vineyard<FS>::
record_knee(Index i)
{
	rLog(rlVineyard, "Entered record_knee()");
	AssertMsg(evaluator != 0, "Cannot record knee with a null evaluator");
	AssertMsg(i->vine() != 0, "Cannot add a knee to a null vine");
	AssertMsg(i->sign(), "record_knee() must be called on a positive simplex");
	
	if (!i->is_paired())
		i->vine()->add(evaluator->value(*i), Infinity, evaluator->time());
	else
	{
		rLog(rlVineyard, "Creating knee");
		KneeS k(evaluator->value(*i), evaluator->value(*(i->pair())), evaluator->time());
		rLog(rlVineyard, "Knee created: %s", tostring(k).c_str());
        rLog(rlVineyard, "Vine: %s", tostring(*(i->vine())).c_str());

		if (!k.is_diagonal() || i->vine()->empty())			// non-diagonal k, or empty vine
		{
			rLog(rlVineyard, "Extending a vine");
			i->vine()->add(k);
		}
		else if (i->vine()->back().is_diagonal())			// last knee is diagonal
		{
			AssertMsg(i->vine()->size() == 1, "Only first knee may be diagonal for a live vine");
			rLog(rlVineyard, "Overwriting first diagonal knee");
			i->vine()->back() = k;
		} else												// finish this vine
		{
			rLog(rlVineyard, "Finishing a vine");
			i->vine()->add(k);
			start_vine(i);
			i->vine()->add(k);
		}
	}
	
	i->vine()->back().set_cycle(resolve_cycle(i));
	rLog(rlVineyard, "Leaving record_knee()");
}

template<class FS>
typename Vineyard<FS>::KneeS::SimplexList  
Vineyard<FS>::
resolve_cycle(Index i) const
{
	rLog(rlVineyard, "Entered resolve_cycle");
	const Cycle& cycle = i->cycle();
	
	// Resolve simplices
	typename KneeS::SimplexList lst;
	for (typename Cycle::const_iterator cur = cycle.begin(); cur != cycle.end(); ++cur)
		lst.push_back(**cur);

	return lst;
}