Added field arithmetic to CohomologyPersistence (Zp + Q) dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Thu, 09 Jul 2009 00:59:32 -0700
branchdev
changeset 137 069596c71902
parent 136 beff535b53ff
child 138 96030f8d2f2c
Added field arithmetic to CohomologyPersistence (Zp + Q)
examples/cohomology/rips-pairwise-cohomology.cpp
examples/cohomology/triangle-cohomology.cpp
include/topology/cohomology-persistence.h
include/topology/cohomology-persistence.hpp
include/topology/field-arithmetic.h
include/utilities/indirect.h
--- a/examples/cohomology/rips-pairwise-cohomology.cpp	Thu May 14 14:04:43 2009 -0700
+++ b/examples/cohomology/rips-pairwise-cohomology.cpp	Thu Jul 09 00:59:32 2009 -0700
@@ -6,6 +6,7 @@
 
 #include <utilities/containers.h>           // for BackInsertFunctor
 #include <utilities/timer.h>
+#include <utilities/log.h>
 
 #include <string>
 
@@ -15,28 +16,36 @@
 typedef     PairwiseDistances<PointContainer, L2Distance>           PairDistances;
 typedef     PairDistances::DistanceType                             DistanceType;
 typedef     PairDistances::IndexType                                Vertex;
-    
-typedef     CohomologyPersistence<DistanceType>                     Persistence;
-typedef     Persistence::SimplexIndex                               Index;
-typedef     Persistence::Death                                      Death;
-
+ 
 typedef     Rips<PairDistances>                                     Generator;
 typedef     Generator::Simplex                                      Smplx;
+typedef     std::vector<Smplx>                                      SimplexVector;
+typedef     SimplexVector::const_iterator                           SV_const_iterator;
 
-typedef     std::map<Smplx, Index, 
-                     Smplx::VertexComparison>                       Complex;
-typedef     std::vector<Smplx>                                      SimplexVector;
+typedef     boost::tuple<Dimension, DistanceType>                   BirthInfo;
+typedef     CohomologyPersistence<BirthInfo, SV_const_iterator>     Persistence;
+typedef     Persistence::SimplexIndex                               Index;
+typedef     Persistence::Death                                      Death;
+typedef     std::map<Smplx, Index, Smplx::VertexComparison>         Complex;
 
-void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance);
+void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime);
 
 int main(int argc, char* argv[])
 {
+#ifdef LOGGING
+    rlog::RLogInit(argc, argv);
+
+    stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
     Dimension               skeleton;
     DistanceType            max_distance;
+    ZpField::Element        prime;
     std::string             infilename;
 
-    program_options(argc, argv, infilename, skeleton, max_distance);
+    program_options(argc, argv, infilename, skeleton, max_distance, prime);
 
+    Timer total_timer; total_timer.start();
     PointContainer          points;
     read_points(infilename, points);
 
@@ -46,30 +55,62 @@
     SimplexVector           v;
     Complex                 c;
     
+    Timer rips_timer; rips_timer.start();
     rips.generate(skeleton, max_distance, make_push_back_functor(v));
     std::sort(v.begin(), v.end(), Generator::Comparison(distances));
+    rips_timer.stop();
     std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
 
     Timer persistence_timer; persistence_timer.start();
-    Persistence p;
+    ZpField                 zp(prime);
+    Persistence             p(zp);
     for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
     {
         std::vector<Index>      boundary;
-        for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); 
-                                     bcur != cur->boundary_end();       ++bcur)
+        for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
             boundary.push_back(c[*bcur]);
         
         Index idx; Death d;
-        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), size(*cur));
-        c[*cur] = idx;
-        if (d && (size(*cur) - *d) > 0)
-            std::cout << (cur->dimension() - 1) << " " << *d << " " << size(*cur) << std::endl;
+        bool store = cur->dimension() < skeleton;
+        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, cur);
+        
+        // c[*cur] = idx;
+        if (store)
+            c.insert(std::make_pair(*cur, idx));
+
+        if (d && (size(*cur) - d->get<1>()) > 0)
+        {
+            AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
+            std::cout << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
+        }
     }
+    // output infinte persistence cocycles
+    for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
+        std::cout << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
     persistence_timer.stop();
+
+
+    // p.show_cocycles();
+    // Output alive cocycles
+    for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
+    {
+        std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
+        for (Persistence::ZColumn::const_iterator zcur = cur->zcolumn.begin(); zcur != cur->zcolumn.end(); ++zcur)
+        {
+            const Smplx& s = **(zcur->si);
+            std::cout << zcur->coefficient << " ";
+            for (Smplx::VertexContainer::const_iterator vcur = s.vertices().begin(); vcur != s.vertices().end(); ++vcur)
+                std::cout << *vcur << " ";
+            std::cout << std::endl;
+        }
+    }
+    total_timer.stop();
+    rips_timer.check("Rips timer");
     persistence_timer.check("Persistence timer");
+    total_timer.check("Total timer");
 }
 
-void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance)
+void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime)
 {
     namespace po = boost::program_options;
 
@@ -81,7 +122,13 @@
     visible.add_options()
         ("help,h",                                                                                  "produce help message")
         ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2),                  "Dimension of the Rips complex we want to compute")
+        ("prime,p",             po::value<ZpField::Element>(&prime)->default_value(2),              "Prime p for the field F_p")
         ("max-distance,m",      po::value<DistanceType>(&max_distance)->default_value(Infinity),    "Maximum value for the Rips complex construction");
+#if LOGGING
+    std::vector<std::string>    log_channels;
+    visible.add_options()
+        ("log,l",               po::value< std::vector<std::string> >(&log_channels),           "log channels to turn on (info, debug, etc)");
+#endif
 
     po::positional_options_description pos;
     pos.add("input-file", 1);
@@ -94,6 +141,11 @@
                   options(all).positional(pos).run(), vm);
     po::notify(vm);
 
+#if LOGGING
+    for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
+        stderrLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
+#endif
+
     if (vm.count("help") || !vm.count("input-file"))
     { 
         std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl;
--- a/examples/cohomology/triangle-cohomology.cpp	Thu May 14 14:04:43 2009 -0700
+++ b/examples/cohomology/triangle-cohomology.cpp	Thu Jul 09 00:59:32 2009 -0700
@@ -54,7 +54,7 @@
 
     stderrLog.subscribeTo(RLOG_CHANNEL("info"));
     stderrLog.subscribeTo(RLOG_CHANNEL("error"));
-    //stderrLog.subscribeTo(RLOG_CHANNEL("topology"));
+    stderrLog.subscribeTo(RLOG_CHANNEL("topology"));
 #endif
 
     SimplexVector v;
@@ -66,10 +66,13 @@
 
     // Compute persistence
     Complex         c;
-    Persistence     p;
+    ZpField         zp(11);
+    Persistence     p(zp);
     unsigned i = 0;
     for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
     {
+        std::cout << "-------" << std::endl;
+
         std::vector<Index>      boundary;
         for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); 
                                      bcur != cur->boundary_end();       ++bcur)
@@ -84,7 +87,7 @@
             // (i.e. when a 1 class kills 0, it's really that in cohomology forward 0 kills 1,
             //  in cohomology backward 1 kills 0, and in homology 1 kills 0)
 
-        //p.show_cycles();
+        p.show_cocycles();
     }
 }
 
--- a/include/topology/cohomology-persistence.h	Thu May 14 14:04:43 2009 -0700
+++ b/include/topology/cohomology-persistence.h	Thu Jul 09 00:59:32 2009 -0700
@@ -16,6 +16,7 @@
 #include <list>
 #include <utility>
 
+#include <topology/field-arithmetic.h>
 #include "utilities/types.h"
 
 #include <boost/optional.hpp>
@@ -23,86 +24,105 @@
 namespace bi = boost::intrusive;
 
 
-template<class BirthInfo_, class SimplexData_ = Empty<> >
+template<class BirthInfo_, class SimplexData_ = Empty<>, class Field_ = ZpField>
 class CohomologyPersistence
 {
     public:
         typedef             BirthInfo_                                                  BirthInfo;
         typedef             SimplexData_                                                SimplexData;
+        typedef             Field_                                                      Field;
+
+        typedef             typename Field::Element                                     FieldElement;
 
 
-        struct SNode;
-        typedef             bi::list<SNode, bi::constant_time_size<false> >             ZRow;
+                            CohomologyPersistence(const Field& field = Field()):
+                                field_(field)                                           {}
+
 
-        // Simplex representation
-        struct SHead: public SimplexData
-        {
-                            SHead(const SHead& other):
-                                SimplexData(other), order(other.order)                  {}  // don't copy row since we can't
-                            SHead(const SimplexData& sd, unsigned o): 
-                                SimplexData(sd), order(o)                               {}
+        // An entry in a cocycle column
+        struct  SNode;      // members: si, coefficient, ci
+        typedef             s::vector<SNode>                                            ZColumn;
+        typedef             bi::list<SNode, bi::constant_time_size<false> >             ZRow;
+        class   CompareSNode;
 
-            // intrusive list corresponding to row of s in Z^*, not ordered in any particular order
-            ZRow            row;
-            unsigned        order;
-        };
-
+        struct  SHead;      // members: row, order
         typedef             s::list<SHead>                                              Simplices;
         typedef             typename Simplices::iterator                                SimplexIndex;
 
-        struct Cocycle;
+        struct  Cocycle;    // members: zcolumn, birth, order
         typedef             s::list<Cocycle>                                            Cocycles;
         typedef             typename Cocycles::iterator                                 CocycleIndex;
-        
-        // An entry in a cocycle column; it's also an element in an intrusive list, hence the list_base_hook<>
-        typedef             bi::list_base_hook<bi::link_mode<bi::auto_unlink> >         auto_unlink_hook;
-        struct SNode: public auto_unlink_hook
-        {
-                            SNode()                                                     {}
-                            SNode(SimplexIndex sidx): si(sidx)                          {}
-
-            // eventually store a field element
-
-            SimplexIndex    si;
-            CocycleIndex    cocycle;                    // TODO: is there no way to get rid of this overhead?
-
-            void            unlink()                    { auto_unlink_hook::unlink(); }
-        };
-        class CompareSNode;
+        typedef             std::pair<CocycleIndex, FieldElement>                       CocycleCoefficientPair;
 
-        typedef             s::vector<SNode>                                            ZColumn;
-        struct Cocycle
-        {
-                            Cocycle(const BirthInfo& b, unsigned o):
-                                birth(b), order(o)                                      {}
-
-            ZColumn         cocycle;
-            BirthInfo       birth;
-            unsigned        order;
-
-            bool            operator<(const Cocycle& other) const                       { return order > other.order; }
-            bool            operator==(const Cocycle& other) const                      { return order == other.order; }
-        };
-
-        typedef             boost::optional<BirthInfo>  Death;
-        typedef             std::pair<SimplexIndex,
-                                      Death>            IndexDeathPair;
+        typedef             boost::optional<BirthInfo>                                  Death;
+        typedef             std::pair<SimplexIndex, Death>                              IndexDeathPair;
 
         // return either a SimplexIndex or a Death
         // BI = BoundaryIterator; it should dereference to a SimplexIndex
         template<class BI>
-        IndexDeathPair      add(BI begin, BI end, BirthInfo b, const SimplexData& sd = SimplexData());
+        IndexDeathPair      add(BI begin, BI end, BirthInfo b, bool store = true, const SimplexData& sd = SimplexData());
 
-        void                show_cycles() const;
-
+        void                show_cocycles() const;
+        CocycleIndex        begin()                                                     { return cocycles_.begin(); }
+        CocycleIndex        end()                                                       { return cocycles_.end(); }
 
     private:
-        void                add_cocycle(Cocycle& z1, Cocycle& z2);
+        void                add_cocycle(CocycleCoefficientPair& z1, CocycleCoefficientPair& z2);
 
     private:
         Simplices           simplices_;
         Cocycles            cocycles_;
+        Field               field_;
 };
+        
+// Simplex representation
+template<class BirthInfo_, class SimplexData_, class Field_>
+struct CohomologyPersistence<BirthInfo_, SimplexData_, Field_>::SHead: public SimplexData
+{
+                    SHead(const SHead& other):
+                        SimplexData(other), order(other.order)                  {}  // don't copy row since we can't
+                    SHead(const SimplexData& sd, unsigned o): 
+                        SimplexData(sd), order(o)                               {}
+
+    // intrusive list corresponding to row of s in Z^*, not ordered in any particular order
+    ZRow            row;
+    unsigned        order;
+};
+
+// An entry in a cocycle column; it's also an element in an intrusive list, hence the list_base_hook<>
+typedef             bi::list_base_hook<bi::link_mode<bi::auto_unlink> >         auto_unlink_hook;
+template<class BirthInfo_, class SimplexData_, class Field_>
+struct CohomologyPersistence<BirthInfo_, SimplexData_, Field_>::SNode: public auto_unlink_hook
+{
+                    SNode(const SNode& other):
+                        si(other.si), coefficient(other.coefficient), 
+                        ci(other.ci)                                            {}
+
+                    SNode(SimplexIndex sidx, FieldElement coef, CocycleIndex cidx): 
+                        si(sidx), coefficient(coef), ci(cidx)                   {}
+
+    SimplexIndex    si;
+    FieldElement    coefficient;
+
+    CocycleIndex    ci;                         // TODO: is there no way to get rid of this overhead?
+
+    void            unlink()                    { auto_unlink_hook::unlink(); }
+};
+
+template<class BirthInfo_, class SimplexData_, class Field_>
+struct CohomologyPersistence<BirthInfo_, SimplexData_, Field_>::Cocycle
+{
+                    Cocycle(const BirthInfo& b, unsigned o):
+                        birth(b), order(o)                                      {}
+
+    ZColumn         zcolumn;
+    BirthInfo       birth;
+    unsigned        order;
+
+    bool            operator<(const Cocycle& other) const                       { return order > other.order; }
+    bool            operator==(const Cocycle& other) const                      { return order == other.order; }
+};
+
 
 #include "cohomology-persistence.hpp"
 
--- a/include/topology/cohomology-persistence.hpp	Thu May 14 14:04:43 2009 -0700
+++ b/include/topology/cohomology-persistence.hpp	Thu Jul 09 00:59:32 2009 -0700
@@ -9,52 +9,62 @@
 static rlog::RLogChannel* rlCohomology =                DEF_CHANNEL("topology/cohomology",        rlog::Log_Debug);
 #endif
 
-template<class BirthInfo, class SimplexData>
-class CohomologyPersistence<BirthInfo, SimplexData>::CompareSNode
+template<class BirthInfo, class SimplexData, class Field>
+class CohomologyPersistence<BirthInfo, SimplexData, Field>::CompareSNode
 {
     public:
         bool        operator()(const SNode& s1, const SNode& s2) const                  { return s1.si->order < s2.si->order; }
 };
 
-template<class BirthInfo, class SimplexData>
+template<class BirthInfo, class SimplexData, class Field>
 template<class BI>
-typename CohomologyPersistence<BirthInfo, SimplexData>::IndexDeathPair
-CohomologyPersistence<BirthInfo, SimplexData>::
-add(BI begin, BI end, BirthInfo birth, const SimplexData& sd)
+typename CohomologyPersistence<BirthInfo, SimplexData, Field>::IndexDeathPair
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+add(BI begin, BI end, BirthInfo birth, bool store, const SimplexData& sd)
 {
     // Create simplex representation
     simplices_.push_back(SHead(sd, simplices_.empty() ? 0 : (simplices_.back().order + 1)));
     SimplexIndex    si = boost::prior(simplices_.end());
 
     // Find out if there are cocycles that evaluate to non-zero on the new simplex
-    typedef         std::list<CocycleIndex>                                             Candidates;
+    typedef         std::list<CocycleCoefficientPair>           Candidates;
     Candidates      candidates, candidates_bulk;
     rLog(rlCohomology, "Boundary");
+    bool sign = true;       // TODO: this is very crude, fix later
     for (BI cur = begin; cur != end; ++cur)
     {
-        rLog(rlCohomology, "  %d", (*cur)->order);
+        rLog(rlCohomology, "  %d %d", (*cur)->order, sign);
         for (typename ZRow::const_iterator zcur = (*cur)->row.begin(); zcur != (*cur)->row.end(); ++zcur)
-            candidates_bulk.push_back(zcur->cocycle);
+            candidates_bulk.push_back(std::make_pair(zcur->ci, 
+                                                     (sign ? zcur->coefficient : field_.neg(zcur->coefficient))));
+        sign = !sign;
     }
 
-    candidates_bulk.sort(make_indirect_comparison(std::less<Cocycle>()));
+    candidates_bulk.sort(make_first_comparison(make_indirect_comparison(std::less<Cocycle>())));
     
+#if LOGGING    
     rLog(rlCohomology,  "  Candidates bulk");
     for (typename Candidates::iterator cur  = candidates_bulk.begin(); 
                                        cur != candidates_bulk.end(); ++cur)
-        rLog(rlCohomology, "    %d", (*cur)->order);
+        rLog(rlCohomology, "    %d %d", cur->first->order, cur->second);
+#endif
 
-    // Remove duplicates --- this is really Z_2, we need a more sophisticated
+    // Remove duplicates
     {
         typename Candidates::const_iterator cur = candidates_bulk.begin();
         while (cur != candidates_bulk.end())
         {
             typename Candidates::const_iterator next = cur;
-            unsigned count = 0;
-            while (next != candidates_bulk.end() && *next == *cur) { ++next; ++count; }
+            FieldElement sum = field_.zero();
+            while (next != candidates_bulk.end() && next->first == cur->first) 
+            { 
+                sum = field_.add(sum, next->second);
+                ++next; 
+            }
     
-            if (count % 2)
-                candidates.push_back(*cur);
+            rLog(rlCohomology, "  Bulk candidate %d sum %d", cur->first->order, sum);
+            if (!field_.is_zero(sum))
+                candidates.push_back(CocycleCoefficientPair(cur->first, sum));
     
             cur = next;
         }
@@ -63,16 +73,22 @@
     // Birth
     if (candidates.empty())
     {
-        rLog(rlCohomology,  "Birth");
+        if (!store)
+        {
+            simplices_.pop_back();
+            return std::make_pair(simplices_.begin(), Death());         // TODO: shouldn't return front
+        }
         
         unsigned order = cocycles_.empty() ? 0 : cocycles_.front().order + 1;
         cocycles_.push_front(Cocycle(birth, order));
+        
+        rLog(rlCohomology,  "Birth: %d", cocycles_.front().order);
 
         // set up the cocycle
-        ZColumn& cocycle = cocycles_.front().cocycle;
-        cocycle.push_back(si);
-        cocycle.front().cocycle = cocycles_.begin();
-        si->row.push_back(cocycles_.front().cocycle.front());
+        ZColumn& cocycle = cocycles_.front().zcolumn;
+        cocycle.push_back(SNode(si, field_.id(), cocycles_.begin()));
+        si->row.push_back(cocycles_.front().zcolumn.front());
+        rLog(rlCohomology,  "  Cocyle: %d", si->order);
 
         return std::make_pair(si, Death());
     }
@@ -80,67 +96,97 @@
     // Death
     rLog(rlCohomology,  "Death");
 
-#if 0
+#if LOGGING
     // Debug only, output candidates
     rLog(rlCohomology,  "  Candidates");
-    for (typename Candidates::iterator cur  = candidates.begin(); 
-                                       cur != candidates.end(); ++cur)
-        rLog(rlCohomology, "    %d", (*cur)->order);
+    for (typename Candidates::iterator cur  = candidates.begin(); cur != candidates.end(); ++cur)
+        rLog(rlCohomology, "    %d %d", cur->first->order, cur->second);
 #endif
 
-    Cocycle& z          = *candidates.front();
-    Death d             = z.birth;
+    CocycleCoefficientPair& z   = candidates.front();
+    Death d                     = z.first->birth;
 
     // add z to everything else in candidates
     for (typename Candidates::iterator cur  = boost::next(candidates.begin()); 
                                        cur != candidates.end(); ++cur)
-        add_cocycle(**cur, z);
+        add_cocycle(*cur, z);
 
-    for (typename ZColumn::iterator cur = z.cocycle.begin(); cur != z.cocycle.end(); ++cur)
+    for (typename ZColumn::iterator cur = z.first->zcolumn.begin(); cur != z.first->zcolumn.end(); ++cur)
         cur->unlink();
     
-    cocycles_.erase(candidates.front());
+    cocycles_.erase(candidates.front().first);
 
     return std::make_pair(si, d);
 }
         
-template<class BirthInfo, class SimplexData>
+template<class BirthInfo, class SimplexData, class Field>
 void
-CohomologyPersistence<BirthInfo, SimplexData>::
-show_cycles() const
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+show_cocycles() const
 {
-    std::cout << "Cocycles" << std::endl;
+    std::cout << "Cocycles: " << cocycles_.size() << std::endl;
     for (typename Cocycles::const_iterator cur = cocycles_.begin(); cur != cocycles_.end(); ++cur)
     {
-        std::cout << cur->order << ": ";
-        for (typename ZColumn::const_iterator zcur = cur->cocycle.begin(); zcur != cur->cocycle.end(); ++zcur)
-            std::cout << zcur->si->order << ", ";
+        // std::cout << cur->order << " (" << cur->birth << "): ";
+        for (typename ZColumn::const_iterator zcur = cur->zcolumn.begin(); zcur != cur->zcolumn.end(); ++zcur)
+            std::cout << zcur->coefficient << " * " << zcur->si->order << ", ";
         std::cout << std::endl;
     }
 }
 
-template<class BirthInfo, class SimplexData>
+template<class BirthInfo, class SimplexData, class Field>
 void
-CohomologyPersistence<BirthInfo, SimplexData>::
-add_cocycle(Cocycle& to, Cocycle& from)
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+add_cocycle(CocycleCoefficientPair& to, CocycleCoefficientPair& from)
 {
-    rLog(rlCohomology,  "Adding cocycle %d to %d", from.order, to.order);
+    rLog(rlCohomology,  "Adding cocycle %d to %d", from.first->order, to.first->order);
+
+    ZColumn         nw;
+    FieldElement    multiplier = field_.neg(field_.div(to.second, from.second));
+    CocycleIndex    ci = to.first;
 
-    ZColumn     nw;
+    typename ZColumn::iterator tcur = to.first->zcolumn.begin();
+    typename ZColumn::iterator fcur = from.first->zcolumn.begin();
+    CompareSNode cmp;
+    while (tcur != to.first->zcolumn.end() && fcur != from.first->zcolumn.end())
+    {
+        rLog(rlCohomology, "  %d %d", tcur->si->order, fcur->si->order);
+        if (cmp(*tcur, *fcur))
+        {
+            nw.push_back(*tcur);
+            ++tcur;
+        }
+        else if (cmp(*fcur, *tcur))
+        {
+            nw.push_back(SNode(fcur->si, field_.mul(multiplier, fcur->coefficient), ci));
+            ++fcur;
+        }
+        else        // equality
+        {
+            FieldElement res = field_.mul(multiplier, tcur->coefficient);
+            res = field_.add(fcur->coefficient, res);
+            if (!field_.is_zero(res))
+                nw.push_back(SNode(fcur->si, res, ci));
+            ++tcur; ++fcur;
+        }
+    }
+    for (; tcur != to.first->zcolumn.end(); ++tcur)
+    {
+        rLog(rlCohomology, "  %d", tcur->si->order);
+        nw.push_back(SNode(*tcur));
+    }
+    for (; fcur != from.first->zcolumn.end(); ++fcur)
+    {
+        rLog(rlCohomology, "  %d", fcur->si->order);
+        nw.push_back(SNode(fcur->si, field_.mul(multiplier, fcur->coefficient), ci));
+    }
 
-    std::set_symmetric_difference(to.cocycle.begin(),   to.cocycle.end(),           // this is catered to Z_2
-                                  from.cocycle.begin(), from.cocycle.end(), 
-                                  std::back_insert_iterator<ZColumn>(nw), 
-                                  CompareSNode());
 
-    for (typename ZColumn::iterator cur = to.cocycle.begin(); cur != to.cocycle.end(); ++cur)
+    for (typename ZColumn::iterator cur = to.first->zcolumn.begin(); cur != to.first->zcolumn.end(); ++cur)
         cur->unlink();
 
-    to.cocycle.swap(nw);
+    to.first->zcolumn.swap(nw);
 
-    for (typename ZColumn::iterator cur = to.cocycle.begin(); cur != to.cocycle.end(); ++cur)
-    {
+    for (typename ZColumn::iterator cur = to.first->zcolumn.begin(); cur != to.first->zcolumn.end(); ++cur)
         cur->si->row.push_back(*cur);
-        cur->cocycle = nw[0].cocycle;
-    }
 }
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/field-arithmetic.h	Thu Jul 09 00:59:32 2009 -0700
@@ -0,0 +1,64 @@
+#ifndef __FIELD_ARITHMETIC_H__
+#define __FIELD_ARITHMETIC_H__
+
+#include <vector>
+#include <gmpxx.h>
+
+class ZpField
+{
+    public:
+        typedef     int                                             Element;
+
+                    ZpField(Element p = 2);
+
+        Element     id()  const                                     { return 1; }
+        Element     zero()  const                                   { return 0; }
+
+        Element     neg(Element a) const                            { return p_ - a; }
+        Element     add(Element a, Element b) const                 { return (a+b) % p_; }
+
+        Element     inv(Element a) const                            { return inverses_[a]; }
+        Element     mul(Element a, Element b) const                 { return (a*b) % p_; }
+        Element     div(Element a, Element b) const                 { return mul(a, inv(b)); }
+
+        bool        is_zero(Element a) const                        { return (a % p_) == 0; }
+
+    private:
+        Element                 p_;
+        std::vector<Element>    inverses_;
+};
+
+ZpField::
+ZpField(Element p):
+    p_(p), inverses_(p_)
+{
+    for (Element i = 1; i < p_; ++i)
+        for (Element j = 1; j < p_; ++j)
+            if (mul(i,j) == 1)
+            {
+                inverses_[i] = j;
+                break;
+            }
+}
+
+class QField
+{
+    public:
+        typedef     mpq_class                                       Element;
+
+                    QField()                                        {}
+
+        Element     id()  const                                     { return 1; }
+        Element     zero()  const                                   { return 0; }
+
+        Element     neg(Element a) const                            { return -a; }
+        Element     add(Element a, Element b) const                 { return (a+b); }
+
+        Element     inv(Element a) const                            { return id()/a; }
+        Element     mul(Element a, Element b) const                 { return (a*b); }
+        Element     div(Element a, Element b) const                 { return a/b; }
+
+        bool        is_zero(Element a) const                        { return a == 0; }
+};
+ 
+#endif // __FIELD_ARITHMETIC_H__
--- a/include/utilities/indirect.h	Thu May 14 14:04:43 2009 -0700
+++ b/include/utilities/indirect.h	Thu Jul 09 00:59:32 2009 -0700
@@ -28,6 +28,27 @@
 { return IndirectComparison<Comparison>(cmp); }
 
 
+template<class Comparison_>
+class FirstComparison
+{
+    public:
+        typedef     Comparison_                             Comparison;
+
+                    FirstComparison(Comparison cmp): 
+                        cmp_(cmp)                           {}
+
+        template<class Pair>
+        bool        operator()(const Pair& a, const Pair& b) const
+        { return cmp_(a.first, b.first); }
+
+    private:
+        Comparison  cmp_;
+};
+template<class Comparison>
+FirstComparison<Comparison> make_first_comparison(const Comparison& cmp)
+{ return FirstComparison<Comparison>(cmp); }
+
+
 template<class Comparison>
 struct ThreeOutcomeCompare: public Comparison
 {