--- a/examples/triangle/CMakeLists.txt Sat Dec 27 14:51:38 2008 -0800
+++ b/examples/triangle/CMakeLists.txt Wed Dec 31 11:51:14 2008 -0800
@@ -1,5 +1,6 @@
set (targets
- triangle)
+ triangle
+ triangle-zigzag)
set (libraries ${libraries} ${Boost_SERIALIZATION_LIBRARY})
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/triangle/triangle-zigzag.cpp Wed Dec 31 11:51:14 2008 -0800
@@ -0,0 +1,87 @@
+#include <vector>
+#include <cassert>
+#include <iostream>
+
+#include <topology/simplex.h>
+#include <topology/zigzag-persistence.h>
+#include <boost/tuple/tuple.hpp>
+
+typedef ZigzagPersistence<unsigned> Zigzag;
+typedef Zigzag::SimplexIndex Index;
+typedef Zigzag::Death Death;
+typedef Zigzag::ZColumn Boundary;
+typedef std::vector<Index> Complex;
+
+int main(int argc, char** argv)
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ //stdoutLog.subscribeTo(RLOG_CHANNEL("topology/persistence"));
+#endif
+
+
+ Zigzag zz;
+ Complex c;
+ Index i; Death d;
+ unsigned birth = 0;
+
+ // Adding the triangle
+ std::cout << birth << ": adding 0" << std::endl;
+ boost::tie(i, d) = zz.add(Boundary(), birth++); // A
+ c.push_back(i);
+ assert(!d); // birth
+
+ std::cout << birth << ": adding 1" << std::endl;
+ boost::tie(i, d) = zz.add(Boundary(), birth++); // B
+ c.push_back(i);
+ assert(!d); // birth
+
+ std::cout << birth << ": adding 2" << std::endl;
+ boost::tie(i, d) = zz.add(Boundary(), birth++); // C
+ c.push_back(i);
+ assert(!d); // birth
+
+ std::cout << birth << ": adding 3" << std::endl;
+ boost::tie(i, d) = zz.add(Boundary(c.begin(),
+ boost::next(c.begin(),2)),
+ birth++); // AB
+ c.push_back(i);
+ assert(d); // death
+ if (d) std::cout << "Death of: " << *d << std::endl;
+
+ std::cout << birth << ": adding 4" << std::endl;
+ boost::tie(i, d) = zz.add(Boundary(boost::next(c.begin()),
+ boost::next(c.begin(),3)),
+ birth++); // BC
+ c.push_back(i);
+ assert(d); // death
+ if (d) std::cout << "Death of: " << *d << std::endl;
+
+ std::cout << birth << ": adding 5" << std::endl;
+ {
+ Boundary bdry; bdry.append(*c.begin(), zz.cmp); bdry.append(*boost::next(c.begin(), 2), zz.cmp);
+ boost::tie(i, d) = zz.add(bdry, birth++); // AC
+ }
+ c.push_back(i);
+ assert(!d); // birth
+
+ std::cout << birth << ": adding 6" << std::endl;
+ boost::tie(i, d) = zz.add(Boundary(boost::next(c.begin(), 3),
+ boost::next(c.begin(), 6)),
+ birth++); // ABC
+ c.push_back(i);
+ assert(d); // death
+ if (d) std::cout << "Death of: " << *d << std::endl;
+
+ //zz.show_all();
+
+ // Removing the triangle in reverse order
+ for (Complex::reverse_iterator cur = c.rbegin(); cur != c.rend(); ++cur)
+ {
+ std::cout << birth << ": removing " << (*cur)->order << std::endl;
+ d = zz.remove(*cur, birth++);
+ if (d) std::cout << "Death of: " << *d << std::endl;
+ //zz.show_all();
+ }
+}
--- a/include/topology/chain.h Sat Dec 27 14:51:38 2008 -0800
+++ b/include/topology/chain.h Wed Dec 31 11:51:14 2008 -0800
@@ -35,6 +35,7 @@
/// \name Template Parameters
/// @{
typedef Container_ Container;
+ typedef SizeStorage<Container> Size;
typedef typename boost::iterator_value<typename Container::iterator>::type Item;
/// @}
@@ -55,6 +56,7 @@
public:
ChainWrapper();
ChainWrapper(const ChainWrapper& c);
+ template<class Iterator> ChainWrapper(Iterator bg, Iterator end);
/// \name Whole Chain operations
/// @{
@@ -67,7 +69,7 @@
template<class ConsistencyComparison>
void sort(const ConsistencyComparison& cmp); ///< Sort elements to enforce consistency
- size_t size() const { return SizeStorage<Container>::size(*this); }
+ size_t size() const { return Size::size(*this); }
using ChainRepresentation::empty;
using ChainRepresentation::clear;
@@ -75,8 +77,8 @@
/// \name Modifiers
/// @{
- using ChainRepresentation::erase;
-
+ void remove(const_reference x);
+
template<class ConsistencyComparison>
void append(const_reference x, const ConsistencyComparison& cmp);
/// @}
@@ -85,6 +87,8 @@
/// @{
using ChainRepresentation::begin;
using ChainRepresentation::end;
+ using ChainRepresentation::front;
+ using ChainRepresentation::back;
template<class OrderComparison>
const_reference top(const OrderComparison& cmp) const; ///< First element in cmp order
@@ -105,6 +109,7 @@
private:
using ChainRepresentation::push_back;
using ChainRepresentation::insert;
+ using ChainRepresentation::erase;
private:
// Serialization
--- a/include/topology/chain.hpp Sat Dec 27 14:51:38 2008 -0800
+++ b/include/topology/chain.hpp Wed Dec 31 11:51:14 2008 -0800
@@ -29,25 +29,42 @@
template<class C>
ChainWrapper<C>::
-ChainWrapper(const ChainWrapper& c): ChainRepresentation(c)
+ChainWrapper(const ChainWrapper& c):
+ ChainRepresentation(c), Size(c)
+{}
+
+template<class C>
+template<class Iterator>
+ChainWrapper<C>::
+ChainWrapper(Iterator bg, Iterator end):
+ ChainRepresentation(bg, end), Size(ChainRepresentation::size())
{}
template<class C>
+void
+ChainWrapper<C>::
+remove(const_reference x)
+{
+ erase(std::find(begin(), end(), x));
+ Size::operator--();
+}
+
+template<class C>
template<class ConsistencyCmp>
void
ChainWrapper<C>::
append(const_reference x, const ConsistencyCmp& cmp)
{
- SizeStorage<Container>::operator++();
+ Size::operator++();
// First try the special cases that x goes at the end
- if (empty() || cmp(ChainRepresentation::back(), x))
- {
+ if (empty() || cmp(back(), x))
push_back(x);
- return;
- }
-
- insert(std::upper_bound(begin(), end(), x, cmp), x);
+ // Then try the special case that x goes at the front
+ else if (cmp(x, front()))
+ ContainerTraits<C,ConsistencyCmp>::push_front(*this, x);
+ else
+ insert(std::upper_bound(begin(), end(), x, cmp), x);
}
template<class C>
@@ -66,7 +83,7 @@
swap(ChainWrapper& c)
{
ChainRepresentation::swap(c);
- SizeStorage<Container>::swap(c);
+ Size::swap(c);
}
template<class C>
@@ -127,7 +144,7 @@
std::set_symmetric_difference(begin(), end(), c.begin(), c.end(), bi, cmp);
static_cast<ChainRepresentation*>(this)->swap(tmp);
- static_cast<SizeStorage<Container>*>(this)->swap(bi);
+ static_cast<Size*>(this)->swap(bi);
return *this;
}
--- a/include/topology/filtration.hpp Sat Dec 27 14:51:38 2008 -0800
+++ b/include/topology/filtration.hpp Wed Dec 31 11:51:14 2008 -0800
@@ -11,7 +11,7 @@
complex_order_map_(bg, reverse_order_.begin()),
simplex_index_map_(bg, end)
{
- std::sort(order_.begin(), order_.end(), IndirectComparison<ComplexIndex, Comparison>(cmp));
+ std::sort(order_.begin(), order_.end(), make_indirect_comparison(cmp));
for (Index obg = order_.begin(), cur = obg; cur != order_.end(); ++cur)
reverse_order_[*cur - bg] = cur - obg;
}
--- a/include/topology/order.h Sat Dec 27 14:51:38 2008 -0800
+++ b/include/topology/order.h Wed Dec 31 11:51:14 2008 -0800
@@ -79,68 +79,34 @@
};
#if 0
-template<class StoragePolicy_ = ListRV<> >
-struct ListOrderDescriptor: public StoragePolicy_::template RebindOrder<typename std::list<ListOrderDescriptor<StoragePolicy_> >::iterator>::other
+template<class Chains_ = VectorChains<>,
+ class Data_ = Empty>
+struct ListOrderDescriptor:
+ public Chains_::template rebding<typename OrderTraits<ListOrderDescriptor<Chains_, Data_> >::Index>::other,
+ public Data_
{
- typedef ListOrderDescriptor<StoragePolicy_> Self;
- typedef typename std::list<Self>::iterator OrderIndex;
- OrderIndex pair;
+ typedef ListOrderDescriptor<Chains_, Data_> Self;
+
+ typedef typename OrderTraits<Self>::Index OrderIndex;
+ typedef typename Chains_::template rebind<OrderIndex>::other Chains;
+ typedef Data_ Data;
+
+ template<class OtherData_> struct RebindData
+ { typedef ListOrderDescriptor<Chains_, OtherData_> other; };
- typedef typename StoragePolicy_::template RebindOrder<OrderIndex>::other StoragePolicy;
- typedef typename StoragePolicy::ComplexIndex ComplexIndex;
-
- ListOrderDescriptor(ComplexIndex i):
- StoragePolicy(i) {}
-
- // Acts as a rebind
- template<class OtherStoragePolicy_> struct Order
- { typedef std::list<ListOrderDescriptor<OtherStoragePolicy_> > type; };
+ OrderIndex pair;
};
-template<class T>
-struct ConsistencyCmp
-{
- int compare(T a, T b) const { if (a < b) return -1; if (a == b) return 0; return 1; }
- bool operator()(T a, T b) const { return compare(a,b) == -1; }
-};
-
-// Traits
-template<class Order_>
-struct OrderTraits
-{
- typedef Order_ Order;
- typedef void StoragePolicy;
- typedef void Index;
- typedef void Element;
-};
-
-// Specialization
-template<class T>
-struct OrderTraits<std::vector<T> >
+// Specialization for ListOrderDescriptor
+template<class Chains, class Data>
+struct OrderTraits<ListOrderDescriptor<Chains, Data> >
{
- typedef std::vector<T> Order;
- typedef typename T::StoragePolicy StoragePolicy;
- typedef typename T::OrderIndex Index;
- typedef T Element;
-
- typedef std::less<Index> OrderComparison;
- typedef ConsistencyCmp<Index> ConsistencyComparison;
-
- static Index begin(Order& order) { return order.begin(); }
- static Index end(Order& order) { return order.end(); }
+ typedef ListOrderDescriptor<Chains, Data> Descriptor;
- template<class Comparison>
- static void sort(Order& order, const Comparison& cmp) { std::sort(order.begin(), order.end(), cmp); }
-};
-
-template<class T>
-struct OrderTraits<std::list<T> >
-{
- typedef std::list<T> Order;
- typedef typename T::StoragePolicy StoragePolicy;
- typedef typename T::OrderIndex Index;
- typedef T Element;
+ typedef std::list<Descriptor> Order;
+ typedef Descriptor Element;
+ typedef Order::iterator Index;
// FIXME
typedef std::less<Index> OrderComparison;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/zigzag-persistence.h Wed Dec 31 11:51:14 2008 -0800
@@ -0,0 +1,138 @@
+#ifndef __ZIGZAG_PERSISTENCE_H__
+#define __ZIGZAG_PERSISTENCE_H__
+
+#include <list>
+#include "cycles.h"
+#include "utilities/types.h"
+#include <sstream>
+
+/**
+ * Class: ZigzagPersistence
+ * TODO: this should probably be parametrized by Chain or Field
+ */
+template<class BirthID_ = Empty<> >
+class ZigzagPersistence
+{
+ public:
+ typedef BirthID_ BirthID;
+
+ struct ZNode;
+ struct BNode;
+ struct SimplexNode;
+
+ typedef std::list<ZNode> ZList;
+ typedef typename ZList::iterator ZIndex;
+ typedef std::list<BNode> BList;
+ typedef typename BList::iterator BIndex;
+ typedef std::list<SimplexNode> SimplexList;
+ typedef typename SimplexList::iterator SimplexIndex;
+
+ // TODO: deques might make a lot of sense, at least for BColumns and ZRows
+ typedef typename VectorChains<ZIndex>::Chain ZRow;
+ typedef typename VectorChains<ZIndex>::Chain BColumn;
+ typedef typename VectorChains<BIndex>::Chain BRow;
+ typedef typename VectorChains<BIndex>::Chain CRow;
+ typedef typename VectorChains<SimplexIndex>::Chain ZColumn;
+ typedef typename VectorChains<SimplexIndex>::Chain CColumn;
+
+ typedef boost::optional<BirthID> Death;
+ typedef std::pair<SimplexIndex, Death> IndexDeathPair;
+
+ // TODO: probably should store something to identify the birth to the outside world; this should probably
+ // be a template parameter (perhaps a template parameter to ZigzagPersistence)
+ struct ZNode
+ {
+ ZNode(int o, const BirthID& b):
+ order(o), birth(b) {}
+
+ int order;
+ ZColumn z_column;
+ BRow b_row;
+ BIndex low; // which BColumn has this ZIndex as low
+
+ BirthID birth; // TODO: need to do empty-member optimization
+ };
+
+ struct BNode
+ {
+ BNode(unsigned o): order(o) {}
+
+ unsigned order;
+ BColumn b_column;
+ CColumn c_column;
+ };
+
+ struct SimplexNode
+ {
+ SimplexNode(unsigned o): order(o) {}
+
+ unsigned order;
+ ZRow z_row;
+ CRow c_row;
+ ZIndex low; // which ZColumn has this SimplexNode as low
+ };
+
+ // Constructor: ZigzagPersistence()
+ ZigzagPersistence() {}
+
+ // Function: add(s)
+ IndexDeathPair add(ZColumn bdry, const BirthID& birth = BirthID());
+
+ // Function: remove(s)
+ Death remove(SimplexIndex s, const BirthID& birth = BirthID());
+
+ // Debug
+ void show_all();
+
+ private:
+ ZList z_list;
+ BList b_list;
+ SimplexList s_list;
+
+ /* Helper functors */
+ template<class Member, class Element> struct Appender;
+ template<class Member, class Element> struct Remover;
+ template<class Member, class Chain> struct Adder;
+
+ template<class Member, class Element>
+ Appender<Member, Element> make_appender(Member m, Element e) const { return Appender<Member, Element>(m,e); }
+ template<class Member, class Element>
+ Remover<Member, Element> make_remover(Member m, Element e) const { return Remover<Member, Element>(m,e); }
+ template<class Member, class Chain>
+ Adder<Member, Chain> make_adder(Member m, Chain& c) const { return Adder<Member, Chain>(m, c); }
+
+ template<class Index, class IndexFrom, class PrimaryMember, class SecondaryMember>
+ void add_chains(Index bg, Index end, IndexFrom j, PrimaryMember pm, SecondaryMember sm);
+ template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
+ void add_chains(IndexTo bg, IndexTo end, IndexFrom j,
+ PrimaryMemberTo pmt, SecondaryMemberTo smt,
+ PrimaryMemberFrom pmf);
+ template<class Index, class PrimaryMember, class SecondaryMember>
+ void add_chain(Index to, Index from,
+ PrimaryMember pmt, SecondaryMember smt);
+ template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
+ void add_chain(IndexTo to, IndexFrom from,
+ PrimaryMemberTo pmt, SecondaryMemberTo smt,
+ PrimaryMemberFrom pmf);
+ template<class IndexTo, class IndexFrom, class PrimaryMember, class SecondaryMember, class DualPrimaryMember, class DualSecondaryMember>
+ void change_basis(IndexTo bg, IndexTo end, IndexFrom j,
+ PrimaryMember pm, SecondaryMember sm,
+ DualPrimaryMember dpm, DualSecondaryMember dsm);
+
+ public:
+ struct OrderComparison
+ {
+ template<class T>
+ bool operator()(T a, T b) const { return a->order < b->order; }
+ } cmp;
+
+ struct OrderOutput
+ {
+ template<class T>
+ std::string operator()(T a) const { std::stringstream s; s << a->order; return s.str(); }
+ } out;
+};
+
+#include "zigzag-persistence.hpp"
+
+#endif // __ZIGZAG_PERSISTENCE_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/zigzag-persistence.hpp Wed Dec 31 11:51:14 2008 -0800
@@ -0,0 +1,403 @@
+#include <utilities/log.h>
+#include <boost/utility.hpp>
+#include <algorithm>
+#include <utilities/indirect.h>
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlZigzagAdd = DEF_CHANNEL("topology/persistence/zigzag/add", rlog::Log_Debug);
+static rlog::RLogChannel* rlZigzagRemove = DEF_CHANNEL("topology/persistence/zigzag/remove", rlog::Log_Debug);
+#endif // LOGGING
+
+
+template<class BID>
+typename ZigzagPersistence<BID>::IndexDeathPair
+ZigzagPersistence<BID>::
+add(ZColumn bdry, const BirthID& birth)
+{
+ rLog(rlZigzagAdd, "Entered ZigzagPersistence::add()");
+ rLog(rlZigzagAdd, " Boundary: %s", bdry.tostring(out).c_str());
+ rLog(rlZigzagAdd, " Boundary size: %d", bdry.size());
+
+ { // scoping to not pollute with the name order
+ unsigned order = s_list.empty() ? 0 : boost::prior(s_list.end())->order + 1;
+ s_list.push_back(SimplexNode(order));
+ }
+ SimplexIndex last_s = boost::prior(s_list.end());
+ last_s->low = z_list.end();
+
+ // Reduce bdry among the cycles
+ BColumn v; // representation of the boundary in the cycle basis
+ while (!bdry.empty())
+ {
+ SimplexIndex l = bdry.back();
+ ZIndex k = l->low;
+ v.append(k, cmp);
+ bdry.add(k->z_column, cmp);
+ }
+
+ // Reduce v among boundaries
+ BRow u;
+ while (!(v.empty()))
+ {
+ ZIndex l = v.back();
+ BIndex k = l->low;
+
+ if (k == b_list.end())
+ break;
+
+ u.append(k, cmp);
+ v.add(k->b_column, cmp);
+ }
+
+ if (v.empty())
+ {
+ // Birth
+ int order = z_list.empty() ? 0 : boost::prior(z_list.end())->order + 1;
+ z_list.push_back(ZNode(order, birth));
+ ZIndex last_z = boost::prior(z_list.end());
+
+ // Set z_column
+ ZColumn& z = last_z->z_column;
+ std::for_each(u.begin(), u.end(), make_adder(&BNode::c_column, z));
+ z.append(last_s, cmp);
+
+ // Set s_row
+ std::for_each(z.begin(), z.end(), make_appender(&SimplexNode::z_row, last_z));
+
+ // Set low
+ last_z->low = b_list.end();
+ last_s->low = last_z;
+
+ return std::make_pair(last_s, Death());
+ } else
+ {
+ // Death
+ unsigned order = b_list.empty() ? 0 : boost::prior(b_list.end())->order + 1;
+ b_list.push_back(BNode(order));
+ BIndex last_b = boost::prior(b_list.end());
+
+ // Set b_column and low
+ last_b->b_column.swap(v);
+ last_b->b_column.back()->low = last_b;
+
+ // Set b_row
+ std::for_each(last_b->b_column.begin(), last_b->b_column.end(), make_appender(&ZNode::b_row, last_b));
+
+ // Set c_column
+ CColumn& c = last_b->c_column;
+ std::for_each(u.begin(), u.end(), make_adder(&BNode::c_column, c));
+ c.append(last_s, cmp);
+
+ // Set c_row
+ std::for_each(c.begin(), c.end(), make_appender(&SimplexNode::c_row, last_b));
+
+ return std::make_pair(last_s, Death(last_b->b_column.back()->birth));
+ }
+}
+
+
+template<class BID>
+typename ZigzagPersistence<BID>::Death
+ZigzagPersistence<BID>::
+remove(SimplexIndex s, const BirthID& birth)
+{
+ rLog(rlZigzagRemove, "Entered ZigzagPersistence::remove(%d)", s->order);
+ rLog(rlZigzagRemove, " s->z_row: %s", s->z_row.tostring(out).c_str());
+
+ if (s->z_row.empty())
+ {
+ AssertMsg(!(s->c_row.empty()), "Birth after removal, row in C must be non-empty");
+
+ // Birth
+ int order = z_list.empty() ? 0 : z_list.begin()->order - 1;
+ z_list.push_front(ZNode(order, birth));
+ ZIndex first_z = z_list.begin();
+ ZColumn& z = first_z->z_column;
+
+ // Prepend DC[j] = ZB[j] to Z
+ BIndex j = s->c_row.front();
+ std::for_each(j->b_column.begin(), j->b_column.end(), make_adder(&ZNode::z_column, z));
+ std::for_each(z.begin(), z.end(), make_appender(&SimplexNode::z_row, first_z));
+
+ // Prepend row of s in C to B
+ first_z->b_row = s->c_row; // copying instead of swapping is inefficient,
+ // but it simplifies logic when subtracting chains later
+ std::for_each(first_z->b_row.begin(), first_z->b_row.end(), make_appender(&BNode::b_column, first_z));
+
+ // Subtract C[j] from every column of C that contains s
+ add_chains(first_z->b_row.begin(), first_z->b_row.end(), j, &BNode::c_column, &SimplexNode::c_row);
+
+ // Subtract B[j] from every other column of B that has l
+ ZIndex l = j->b_column.back();
+ BRow l_row = l->b_row; // again, copying is inefficient, but it simplifies the logic
+ add_chains(l_row.begin(), l_row.end(), j, &BNode::b_column, &ZNode::b_row);
+
+ // Drop j, l, and s
+ // TODO: add some assertions, like everything being empty
+ //
+ // l->z_column is the only non-empty thing, but we drop it,
+ // the basis is preserved because we added first_z
+ l->z_column.back()->low = z_list.end();
+ b_list.erase(j);
+ std::for_each(l->z_column.begin(), l->z_column.end(), make_remover(&SimplexNode::z_row, l));
+ z_list.erase(l);
+ s_list.erase(s);
+
+ // Reduce Z
+ SimplexIndex ls = first_z->z_column.back();
+ while(ls->low != first_z)
+ {
+ if (ls->low == z_list.end()) { ls->low = first_z; break; }
+
+ // if ls->low precedes first_z, swap them
+ if (cmp(ls->low, first_z)) std::swap(ls->low, first_z);
+
+ add_chain(first_z, ls->low, &ZNode::z_column, &SimplexNode::z_row);
+ std::swap(ls->low, first_z);
+ }
+
+ return Death();
+ } else
+ {
+ // Death
+ ZIndex j = s->z_row.front();
+ CRow c_row = s->c_row;
+
+ // Subtract Z[j] from every chain in C that contains s
+ add_chains(c_row.begin(), c_row.end(), j, &BNode::c_column, &SimplexNode::c_row, &ZNode::z_column);
+
+ /**
+ * Change basis to remove s from Z
+ * All the processing is done from back to front, so that when elements are removed from s->z_row,
+ * it doesn't invalidate the rest of the elements (iterators) in the row.
+ */
+ // Compute reducers --- columns that we will be adding to other columns
+ typedef typename ZRow::reverse_iterator ZRowReverseIterator;
+ typedef std::list<ZRowReverseIterator> ReducersContainer;
+ ReducersContainer reducers; // list of ZColumns that we will be adding to other columns
+ reducers.push_back(boost::prior(s->z_row.rend())); // j is the first reducer
+ SimplexIndex low = j->z_column.back();
+ for (typename ZRow::iterator cur = s->z_row.begin(); cur != s->z_row.end(); ++cur)
+ if (cmp((*cur)->z_column.back(), low))
+ {
+ reducers.push_back(ZRowReverseIterator(cur));
+ low = (*(reducers.back()))->z_column.back();
+ }
+ reducers.push_back(s->z_row.rbegin());
+ //rLog(rlZigzagRemove, " Reducers size: %d", reducers.size());
+
+ // Add each reducer to the columns that follow them until the next reducer
+ for (typename ReducersContainer::reverse_iterator cur = boost::next(reducers.rbegin()); cur != reducers.rend(); ++cur)
+ change_basis(*boost::prior(cur), *cur, **cur,
+ &ZNode::z_column, &SimplexNode::z_row,
+ &ZNode::b_row, &BNode::b_column);
+
+ // Drop j and s
+ // TODO: add some assertions
+ BirthID birth = j->birth;
+ std::for_each(j->z_column.begin(), j->z_column.end(), make_remover(&SimplexNode::z_row, j));
+ z_list.erase(j);
+ s_list.erase(s);
+
+ return Death(birth);
+ }
+}
+
+template<class BID>
+void
+ZigzagPersistence<BID>::
+show_all()
+{
+ std::cout << "s_list:" << std::endl;
+ for (SimplexIndex cur = s_list.begin(); cur != s_list.end(); ++cur)
+ {
+ std::cout << " " << cur->order << ":" << std::endl;
+
+ std::cout << " z_row: ";
+ for (typename ZRow::const_iterator zcur = cur->z_row.begin(); zcur != cur->z_row.end(); ++zcur)
+ std::cout << (*zcur)->order << " ";
+ std::cout << std::endl;
+
+ std::cout << " c_row: ";
+ for (typename CRow::const_iterator ccur = cur->c_row.begin(); ccur != cur->c_row.end(); ++ccur)
+ std::cout << (*ccur)->order << " ";
+ std::cout << std::endl;
+ }
+
+ std::cout << "z_list:" << std::endl;
+ for (ZIndex cur = z_list.begin(); cur != z_list.end(); ++cur)
+ {
+ std::cout << " " << cur->order << ":" << std::endl;
+
+ std::cout << " birth: " << cur->birth << std::endl;
+
+ std::cout << " z_column: ";
+ for (typename ZColumn::const_iterator zcur = cur->z_column.begin(); zcur != cur->z_column.end(); ++zcur)
+ std::cout << (*zcur)->order << " ";
+ std::cout << std::endl;
+
+ std::cout << " b_row: ";
+ for (typename BRow::const_iterator bcur = cur->b_row.begin(); bcur != cur->b_row.end(); ++bcur)
+ std::cout << (*bcur)->order << " ";
+ std::cout << std::endl;
+
+ }
+
+ std::cout << "b_list:" << std::endl;
+ for (BIndex cur = b_list.begin(); cur != b_list.end(); ++cur)
+ {
+ std::cout << " " << cur->order << ":" << std::endl;
+
+ std::cout << " b_column: ";
+ for (typename BColumn::const_iterator bcur = cur->b_column.begin(); bcur != cur->b_column.end(); ++bcur)
+ std::cout << (*bcur)->order << " ";
+ std::cout << std::endl;
+
+ std::cout << " c_column: ";
+ for (typename CColumn::const_iterator ccur = cur->c_column.begin(); ccur != cur->c_column.end(); ++ccur)
+ std::cout << (*ccur)->order << " ";
+ std::cout << std::endl;
+ }
+}
+
+/* Private */
+
+// Class: Appender
+//
+// Functor that appends given element to the given member of whatever parameter it is invoked with
+template<class BID>
+template<class Member, class Element>
+struct ZigzagPersistence<BID>::Appender
+{
+ Appender(Member mm, Element ee):
+ m(mm), e(ee) {}
+
+ template<class T>
+ void operator()(T& a) { ((*a).*m).append(e, cmp); }
+
+ Member m;
+ Element e;
+ OrderComparison cmp;
+};
+
+// Class: Remover
+//
+// Functor that removes given element from the given member of whatever parameter it is invoked with
+template<class BID>
+template<class Member, class Element>
+struct ZigzagPersistence<BID>::Remover
+{
+ Remover(Member mm, Element ee):
+ m(mm), e(ee) {}
+
+ template<class T>
+ void operator()(T& a) { ((*a).*m).remove(e); }
+
+ Member m;
+ Element e;
+};
+
+// Class: Adder
+//
+// Functor that adds the given member of whatever it is invoked with to the given chain
+template<class BID>
+template<class Member, class Chain>
+struct ZigzagPersistence<BID>::Adder
+{
+ Adder(Member mm, Chain& cc):
+ m(mm), c(cc) {}
+
+ template<class T>
+ void operator()(T& a) { c.add((*a).*m, cmp); }
+
+ Member m;
+ Chain& c;
+ OrderComparison cmp;
+};
+
+
+// Function: add_chains()
+//
+// Special case of add_chains where all Indexes are the same, and
+// therefore PrimaryMemberFrom and PrimaryMemberTo are the same
+template<class BID>
+template<class Index, class IndexFrom, class PrimaryMember, class SecondaryMember>
+void
+ZigzagPersistence<BID>::
+add_chains(Index bg, Index end, IndexFrom j, PrimaryMember pm, SecondaryMember sm)
+{
+ add_chains(bg, end, j, pm, sm, pm);
+}
+
+// Function: add_chains()
+//
+// Adds PrimaryMember pm of j to pm of every element in the range [bg,end)
+// Fixes SecondaryMembers by adding and removing the corresponding elements.
+// For example, if we add a column to a number of other columns, then PrimaryMember is that
+// column member, and SecondaryMember is the corresponding row member.
+template<class BID>
+template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
+void
+ZigzagPersistence<BID>::
+add_chains(IndexTo bg, IndexTo end, IndexFrom j, PrimaryMemberTo pmt, SecondaryMemberTo smt, PrimaryMemberFrom pmf)
+{
+ for (IndexTo cur = bg; cur != end; ++cur)
+ add_chain(*cur, j, pmt, smt, pmf);
+}
+
+// Function: change_basis()
+//
+// Changes basis by adding PrimaryMember pm of j to pm of every element in range [bg, end).
+// In parallel it performs the reverse (complementary) update on the dual members, i.e.
+// column and row operations are performed in sync, so that the product of the two matrices doesn't change
+template<class BID>
+template<class IndexTo, class IndexFrom, class PrimaryMember, class SecondaryMember, class DualPrimaryMember, class DualSecondaryMember>
+void
+ZigzagPersistence<BID>::
+change_basis(IndexTo bg, IndexTo end, IndexFrom j, PrimaryMember pm, SecondaryMember sm, DualPrimaryMember dpm, DualSecondaryMember dsm)
+{
+ for (IndexTo cur = bg; cur != end; ++cur)
+ {
+ add_chain(*cur, j, pm, sm, pm);
+ add_chain(j, *cur, dpm, dsm, dpm);
+ }
+}
+
+template<class BID>
+template<class Index, class PrimaryMember, class SecondaryMember>
+void
+ZigzagPersistence<BID>::
+add_chain(Index to, Index from, PrimaryMember pm, SecondaryMember sm)
+{
+ add_chain(to, from, pm, sm, pm);
+}
+
+// Function: add_chain()
+//
+// Adds PrimaryMemberFrom pmf of `from` to PrimaryMemberTo pmt of `to`.
+// Fixes SecondaryMemberTos. See add_chains().
+template<class BID>
+template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
+void
+ZigzagPersistence<BID>::
+add_chain(IndexTo to, IndexFrom from, PrimaryMemberTo pmt, SecondaryMemberTo smt, PrimaryMemberFrom pmf)
+{
+ // Fix secondaries
+ std::for_each(make_intersection_iterator(((*from).*pmf).begin(), ((*from).*pmf).end(),
+ ((*to).*pmt).begin(), ((*to).*pmt).end(),
+ cmp),
+ make_intersection_iterator(((*from).*pmf).end(), ((*from).*pmf).end(),
+ ((*to).*pmt).end(), ((*to).*pmt).end(),
+ cmp),
+ make_remover(smt, to));
+ std::for_each(make_difference_iterator(((*from).*pmf).begin(), ((*from).*pmf).end(),
+ ((*to).*pmt).begin(), ((*to).*pmt).end(),
+ cmp),
+ make_difference_iterator(((*from).*pmf).end(), ((*from).*pmf).end(),
+ ((*to).*pmt).end(), ((*to).*pmt).end(),
+ cmp),
+ make_appender(smt, to));
+
+ // Add primaries
+ ((*to).*pmt).add((*from).*pmf, cmp);
+}
--- a/include/utilities/containers.h Sat Dec 27 14:51:38 2008 -0800
+++ b/include/utilities/containers.h Wed Dec 31 11:51:14 2008 -0800
@@ -11,10 +11,12 @@
{
typedef Container_ Container;
typedef typename Container::value_type Item;
+ typedef typename Container::const_reference const_reference;
typedef Comparison_ Comparison;
static void reserve(Container& c, size_t sz) {}
static void sort(Container& c, const Comparison& cmp = Comparison()) { c.sort(cmp); }
+ static void push_front(Container& c, const_reference x) { c.push_front(x); }
};
/**
@@ -75,10 +77,12 @@
{
typedef T Item;
typedef std::vector<T> Container;
+ typedef typename Container::const_reference const_reference;
typedef Comparison_ Comparison;
static void reserve(Container& c, size_t sz) { c.reserve(sz); }
static void sort(Container& c, const Comparison& cmp = Comparison()) { std::sort(c.begin(), c.end(), cmp); }
+ static void push_front(Container& c, const_reference x) { c.insert(c.begin(), x); }
};
template<class T, class Comparison_>
@@ -86,6 +90,7 @@
{
typedef T Item;
typedef List<T> Container;
+ typedef typename Container::const_reference const_reference;
typedef Comparison_ Comparison;
static void reserve(Container& c, size_t sz) { }
@@ -95,6 +100,7 @@
std::sort(tmp.begin(), tmp.end(), cmp);
std::copy(tmp.begin(), tmp.end(), c.begin());
}
+ static void push_front(Container& c, const_reference x) { c.push_front(x); }
};
// TODO: specialize for List (singly-linked list)
--- a/include/utilities/indirect.h Sat Dec 27 14:51:38 2008 -0800
+++ b/include/utilities/indirect.h Wed Dec 31 11:51:14 2008 -0800
@@ -2,39 +2,139 @@
#define __INDIRECT_H__
#include <boost/iterator/iterator_adaptor.hpp>
+#include <boost/iterator/iterator_facade.hpp>
+#include <iterator>
// TODO: write documentation
-template<class Iterator_, class Comparison_>
+template<class Comparison_>
struct IndirectComparison
{
- typedef Iterator_ Iterator;
typedef Comparison_ Comparison;
IndirectComparison(const Comparison& cmp):
cmp_(cmp)
{}
+ template<class Iterator>
bool operator()(Iterator a, Iterator b) const
{ return cmp_(*a, *b); }
const Comparison& cmp_;
};
+template<class Comparison>
+IndirectComparison<Comparison> make_indirect_comparison(const Comparison& cmp)
+{ return IndirectComparison<Comparison>(cmp); }
+
template<class Comparison>
struct ThreeOutcomeCompare: public Comparison
{
- typedef typename Comparison::first_argument_type first_argument_type;
- typedef typename Comparison::second_argument_type second_argument_type;
+ typedef typename Comparison::first_argument_type first_argument_type;
+ typedef typename Comparison::second_argument_type second_argument_type;
- ThreeOutcomeCompare(const Comparison& cmp = Comparison()): Comparison(cmp) {}
+ ThreeOutcomeCompare(const Comparison& cmp = Comparison()): Comparison(cmp) {}
- int compare(const first_argument_type& a, const second_argument_type& b) const
+ int compare(const first_argument_type& a, const second_argument_type& b) const
{ if (operator()(a,b)) return -1;
else if (operator()(b,a)) return 1;
else return 0;
}
};
+// Iterates over the difference of the two sorted sequences, dereferencing into the first sequence
+template<class Iterator1, class Iterator2, class StrictWeakOrdering>
+class difference_iterator: public boost::iterator_facade<difference_iterator<Iterator1, Iterator2, StrictWeakOrdering>,
+ typename std::iterator_traits<Iterator1>::value_type,
+ boost::forward_traversal_tag>
+{
+ public:
+ typedef typename std::iterator_traits<Iterator1>::reference reference;
+
+ difference_iterator(Iterator1 cur1, Iterator1 end1,
+ Iterator2 cur2, Iterator2 end2,
+ const StrictWeakOrdering& cmp = StrictWeakOrdering()):
+ cur1_(cur1), end1_(end1), cur2_(cur2), end2_(end2), cmp_(cmp) { catchup(); }
+
+ private:
+ friend class boost::iterator_core_access;
+
+ void increment() { ++cur1_; catchup(); }
+ bool equal(const difference_iterator& other) const { return (cur1_ == other.cur1_) && (cur2_ == other.cur2_); }
+ reference dereference() const { return *cur1_; }
+
+ private:
+ Iterator1 cur1_, end1_;
+ Iterator2 cur2_, end2_;
+ const StrictWeakOrdering& cmp_;
+
+ private:
+ void catchup()
+ {
+ while ((cur1_ != end1_) && (cur2_ != end2_))
+ {
+ if (cmp_(*cur1_, *cur2_)) break;
+ else if (cmp_(*cur2_, *cur1_)) ++cur2_;
+ else { ++cur1_; ++cur2_; }
+ }
+
+ if (cur1_ == end1_) cur2_ = end2_;
+ }
+};
+
+template<class Iterator1, class Iterator2, class StrictWeakOrdering>
+difference_iterator<Iterator1, Iterator2, StrictWeakOrdering>
+make_difference_iterator(Iterator1 cur1, Iterator1 end1, Iterator2 cur2, Iterator2 end2, const StrictWeakOrdering& cmp)
+{ return difference_iterator<Iterator1, Iterator2, StrictWeakOrdering>(cur1, end1, cur2, end2, cmp); }
+
+// Iterates over the intersection of the two sorted sequences, dereferencing into the first sequence
+template<class Iterator1, class Iterator2, class StrictWeakOrdering>
+class intersection_iterator: public boost::iterator_facade<intersection_iterator<Iterator1, Iterator2, StrictWeakOrdering>,
+ typename std::iterator_traits<Iterator1>::value_type,
+ boost::forward_traversal_tag>
+{
+ public:
+ typedef typename std::iterator_traits<Iterator1>::reference reference;
+
+ intersection_iterator(Iterator1 cur1, Iterator1 end1,
+ Iterator2 cur2, Iterator2 end2,
+ const StrictWeakOrdering& cmp = StrictWeakOrdering()):
+ cur1_(cur1), end1_(end1), cur2_(cur2), end2_(end2), cmp_(cmp) { catchup(); }
+
+ private:
+ friend class boost::iterator_core_access;
+
+ void increment() { ++cur1_; ++cur2_; catchup(); }
+ bool equal(const intersection_iterator& other) const { return (cur1_ == other.cur1_) && (cur2_ == other.cur2_); }
+ reference dereference() const { return *cur1_; }
+
+ private:
+ Iterator1 cur1_, end1_;
+ Iterator2 cur2_, end2_;
+ const StrictWeakOrdering& cmp_;
+
+ private:
+ void catchup()
+ {
+ while ((cur1_ != end1_) && (cur2_ != end2_))
+ {
+ if (cmp_(*cur1_, *cur2_)) ++cur1_;
+ else if (cmp_(*cur2_, *cur1_)) ++cur2_;
+ else break;
+ }
+
+ if ((cur1_ == end1_) || (cur2_ == end2_))
+ {
+ cur1_ = end1_;
+ cur2_ = end2_;
+ }
+ }
+};
+
+template<class Iterator1, class Iterator2, class StrictWeakOrdering>
+intersection_iterator<Iterator1, Iterator2, StrictWeakOrdering>
+make_intersection_iterator(Iterator1 cur1, Iterator1 end1, Iterator2 cur2, Iterator2 end2, const StrictWeakOrdering& cmp)
+{ return intersection_iterator<Iterator1, Iterator2, StrictWeakOrdering>(cur1, end1, cur2, end2, cmp); }
+
#endif // __INDIRECT_H__
--- a/tests/utilities/CMakeLists.txt Sat Dec 27 14:51:38 2008 -0800
+++ b/tests/utilities/CMakeLists.txt Wed Dec 31 11:51:14 2008 -0800
@@ -1,4 +1,5 @@
set (targets
+ test-set-iterators
test-consistencylist
test-orderlist
test-counters)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/utilities/test-set-iterators.cpp Wed Dec 31 11:51:14 2008 -0800
@@ -0,0 +1,59 @@
+#include <vector>
+#include <list>
+#include <algorithm>
+#include <iostream>
+
+#include <utilities/indirect.h>
+
+#include <boost/lambda/lambda.hpp>
+using boost::lambda::_1;
+
+int main()
+{
+ std::vector<int> v;
+ v.push_back(1);
+ v.push_back(3);
+ v.push_back(5);
+ v.push_back(7);
+ v.push_back(9);
+ std::cout << "v: ";
+ std::for_each(v.begin(), v.end(), std::cout << _1 << ' ');
+ std::cout << std::endl;
+
+ std::list<int> l;
+ l.push_back(2);
+ l.push_back(3);
+ l.push_back(4);
+ l.push_back(5);
+ l.push_back(6);
+ l.push_back(8);
+ std::cout << "l: ";
+ std::for_each(l.begin(), l.end(), std::cout << _1 << ' ');
+ std::cout << std::endl;
+ std::cout << std::endl;
+
+ std::cout << "v \\cap l: ";
+ std::for_each(make_intersection_iterator(v.begin(), v.end(), l.begin(), l.end(), std::less<int>()),
+ make_intersection_iterator(v.end(), v.end(), l.end(), l.end(), std::less<int>()),
+ std::cout << _1 << ' ');
+ std::cout << std::endl;
+
+ std::cout << "v - l: ";
+ std::for_each(make_difference_iterator(v.begin(), v.end(), l.begin(), l.end(), std::less<int>()),
+ make_difference_iterator(v.end(), v.end(), l.end(), l.end(), std::less<int>()),
+ std::cout << _1 << ' ');
+ std::cout << std::endl;
+ std::cout << std::endl;
+
+ std::cout << "l \\cap v: ";
+ std::for_each(make_intersection_iterator(l.begin(), l.end(), v.begin(), v.end(), std::less<int>()),
+ make_intersection_iterator(l.end(), l.end(), v.end(), v.end(), std::less<int>()),
+ std::cout << _1 << ' ');
+ std::cout << std::endl;
+
+ std::cout << "l - v: ";
+ std::for_each(make_difference_iterator(l.begin(), l.end(), v.begin(), v.end(), std::less<int>()),
+ make_difference_iterator(l.end(), l.end(), v.end(), v.end(), std::less<int>()),
+ std::cout << _1 << ' ');
+ std::cout << std::endl;
+}