SimplicialComplex.h

Go to the documentation of this file.

#ifndef SimplicialComplex_HEADER
#define SimplicialComplex_HEADER

/***************************************************************************
 * Benoit Hudson (C) 2007
 ***************************************************************************/

#include <boost/array.hpp>
#include <boost/iterator/transform_iterator.hpp>
#include <utilities/IntrusivePtr.h>
#include <utilities/IntrusiveList.h>
#include <utilities/FastHash.h>
#include <utilities/FrontStack.h>
#include <utilities/FrontQueue.h>
#include <utilities/ArrayFunctions.h>
#include <utilities/hash.h>
#include <utilities/option.h>
#include <utilities/SmallList.h>
#include <utilities/Void.h>
#include <utilities/foreach.h>
#include <stack>
#include <ext/slist>

#if 0
#define dprintf(...) printf(__VA_ARGS__)
#else
#define dprintf(...)
#endif

#define FRINGE_USES_INTRUSIVE_LIST

/***************************************************************************
 ** An arbitrary-dimensional simplicial complex, backed internally by a
 ** link structure.  Operations occur in the obvious amount of time (constant
 ** for traversing from one face to the next, or for adding or removing
 ** a simplex to the complex), assuming constant dimension.
 **
 ** We have an external point, represented as NULL.  All real vertices must be
 ** pointers, and they must be non-NULL.
 **
 ** Simplices have data on them, which saves the user from holding their
 ** own hash tables.
 **
 ** Simplices are consistently oriented.  To get finer control of orientation,
 ** use OSimplex.
 **
 ** The template arguments have underscores so that I can typedef them
 ** and reveal them.
 **
 ** Important revisions:
 ** - r301: best simplex-dictionary complex
 ** - r280: first working pointer-based complex; too slow, backed out.
 ** - r276: first working simplex-dictionary complex
 ** - r270: best face-dictionary complex
 ***************************************************************************/
template <unsigned d, class Vertex_, class SimplexData_ = void, class VertexPrinter_ = typename Vertex_::Printer>
class SimplicialComplex {
struct Flipper; // for friendship purposes.
public:
class OSimplex; // Oriented simplex.  See below.
class Simplex;  // Unoriented simplex.  See below.
class Face;     // Unoriented (d-1)-facet.  See below.
struct Star;    // Acts just like a vector<Simplex>.  See below.
struct OptSimplex;  // for friendship

typedef Vertex_ Vertex; 
typedef VertexPrinter_ VertexPrinter; 
typedef SimplexData_ SimplexData; 
typedef typename hudson::PayloadBucket<size_t, SimplexData>::reference data_ref; 
typedef typename hudson::PayloadBucket<size_t, SimplexData>::const_reference data_const_ref; 
/***************************************************************************
 * The internal class that implements a simplex.  Private.
 *
 * We don't publicize this class' existence.  Therefore, everything is public,
 * since the user will never see an ISimplex, and so that I don't have to
 * befriend every class in this file.
 ***************************************************************************/
private:
struct ISimplex : public hudson::SmallMemoryItem<ISimplex> {
    friend struct Flipper;

    boost::array<Vertex*, d+1> verts_; // canonically oriented (not sorted)
    boost::array<ISimplex*, d+1> neigh_;

    // Intrusive-list pointers for searches.  We need two:
    // one for the visited set, one for the fringe.
    ISimplex *nextVisited_;
#ifdef FRINGE_USES_INTRUSIVE_LIST
    ISimplex *nextFringe_;
#endif

    // refcount, then payload goes last; pair them up so that a void
    // user-data item costs zero storage.
    hudson::PayloadBucket<unsigned, SimplexData> refcount_and_user_;

    unsigned& refcount() { return refcount_and_user_.data; }

    // Status of the simplex -- born (in a star); in the complex; or removed
    // (but not yet deallocated).  Stored using the nextVisited_ pointer, since
    // only an inserted simplex can be visited.
    enum status {
      BORN = -2,
      INSERTED = hudson::IntrusiveListDefaultOut,
      DEAD = -3
    };
    BOOST_STATIC_ASSERT(hudson::IntrusiveListDefaultOut != BORN);
    BOOST_STATIC_ASSERT(hudson::IntrusiveListDefaultOut != DEAD);

    void init() {
#ifndef NDEBUG
      BOOST_FOREACH(const Vertex *v, verts_) { assert(v); }
#endif
      neigh_.assign(NULL);

      refcount() = 0;
      // This gets clobbered by setStatus() below.
      // nextVisited_ = (ISimplex*)hudson::IntrusiveListDefaultOut;
#ifdef FRINGE_USES_INTRUSIVE_LIST
      nextFringe_ = (ISimplex*)hudson::IntrusiveListDefaultOut;
#endif
      setStatus(BORN);
    }

    /* Initialize using the exactly d+1 vertices in the given vector. */
    ISimplex(const vector<Vertex*>& vs) {
      assert(vs.size() == d + 1);
      for(size_t i = 0; i < d + 1; ++i) {
        verts_[i] = vs[i];
      }
      init();
    }

    /* Initialize using the exactly d+1 vertices at the given indices of
       the vertex vector.  I assert the */
    ISimplex(const vector<Vertex*>& vs,
             const std::vector<size_t>& ids,
             ssize_t zeroindex) {
#ifndef NDEBUG
      size_t n = vs.size(); // size() doesn't always get inlined.
#endif
      assert(ids.size() == d + 1);
      for(unsigned i = 0; i < d + 1; ++i) {
        unsigned ii = ids[i] - zeroindex;
        assert(ii < n);
        verts_[i] = vs[ii];
      }
      init();
    }

    /* Initialize given the array, which is known at compile time to have
       the proper size. */
    ISimplex(const boost::array<Vertex*, d+1>& vs) : verts_(vs) {
      init();
    }

    ~ISimplex() { dprintf("deleting %s\n", toString().c_str()); }

    void grab() { refcount()++; }
    void drop() {
      refcount()--;
      if (refcount() == 0) {
        delete this;
      }
    }
    static void add_ref(ISimplex *s) { s->grab(); }
    static void release(ISimplex *s) { s->drop(); }

    unsigned indexOf(const ISimplex *neigh) const {
      for (unsigned i = 0; i <= d; ++i) {
        if (neigh_[i] == neigh) {
          return i;
        }
      }
      abort();
    }
    Vertex *newApex(const ISimplex *origin) const {
      return verts_[indexOf(origin)];
    }

    size_t hash() const {
      return hudson::hashArray(verts_, hudson::hash_by_cast<Vertex*>());
    }

    std::string toString() const {
      string s = "<";
      VertexPrinter print;
      BOOST_FOREACH(const Vertex *v, verts_) {
        if (v) {
          s += " ";
          s += print(v);
        } else {
          s += " (null)";
        }
      }
      s += " >";
      return s;
    }

    Vertex *operator[] (unsigned i) const {
      return get(i);
    }
    Vertex *get(unsigned i) const {
      return verts_[i];
    }

    ISimplex *getN(unsigned i) const {
      return neigh_[i];
    }

    bool has(const Vertex *q) const {
      BOOST_FOREACH(const Vertex *v, verts_) {
        if (v == q) { return true; }
      }
      return false;
    }
    template <class BinaryPredicate>
    bool has(const Vertex *q, const BinaryPredicate& pred) const {
      BOOST_FOREACH(const Vertex *v, verts_) {
        if (pred(v, q)) { return true; }
      }
      return false;
    }
    unsigned indexOf(const Vertex *v) const {
      for(unsigned i = 0; i <= d; ++i) {
        if (get(i) == v) { return i; }
      }
      abort();
    }

    const boost::array<Vertex*, d+1>& toArray() const {
      return verts_;
    }
    typedef typename boost::array<Vertex*, d+1>::const_iterator vertex_iterator;
    vertex_iterator begin_vertices() const {
      return verts_.begin();
    }
   vertex_iterator end_vertices() const {
      return verts_.end();
    }

    status getStatus() const {
      switch((ssize_t)nextVisited_) {
        case -2: return BORN;
        case -3: return DEAD;
        default: return INSERTED;
      }
    }
    void setStatus(status s) {
      nextVisited_ = (ISimplex*)s;
    }
};

private:
/****************************************
 * Handle the IntrusiveList for marking
 * simplices visited.
 ****************************************/
struct ILAdaptor_visited {
  static ISimplex *getNext(ISimplex *is) { return is->nextVisited_; }
  static void setNext(ISimplex *is, ISimplex *next) { is->nextVisited_ = next; }
};

#ifdef FRINGE_USES_INTRUSIVE_LIST
/****************************************
 * Handle the IntrusiveList for putting
 * simplices on the search fringe.
 ****************************************/
struct ILAdaptor_fringe {
  static ISimplex *getNext(ISimplex *is) { return is->nextFringe_; }
  static void setNext(ISimplex *is, ISimplex *next) { is->nextFringe_ = next; }
};
#endif

/***************************************************************************
 * A class that represents the ith face of a simplex (across from the ith
 * vertex).
 *
 * Not ref-counted, so we don't leak this to the caller; it's used internally
 * during the modification routines.
 ***************************************************************************/
private:
struct IFace {
  ISimplex *simplex_;
  unsigned face_;
  IFace() : simplex_(NULL) { }
  IFace(ISimplex *s, unsigned i): simplex_(s), face_(i) { }
  IFace(const OSimplex& os): simplex_(os.get()), face_(simplex_->indexOf(os[d])) { }
  IFace(const OSimplex& os, unsigned i): simplex_(os.get()), face_(simplex_->indexOf(os[i])) { }

  ISimplex *getN() const { return simplex_->getN(face_); }
  IFace getNFace() const {
    ISimplex *n = getN();
    if (n) {
      return IFace(n, n->indexOf(simplex_));
    } else {
      return IFace();
    }
  }
  operator bool() const { return simplex_ != NULL; }
  bool operator== (const ISimplex *other) const {
    return simplex_ == other;
  }
  bool operator!= (const ISimplex *other) const {
    return simplex_ != other;
  }
  string toString() const {
    if (simplex_) {
      char buffer[4000];
      sprintf(buffer, " face %u", face_);
      return simplex_->toString() + buffer;
    } else {
      return "(nil)";
    }
  }

  void relink(const IFace& other) {
    simplex_->neigh_[face_] = other.simplex_;
  }
  void relink(ISimplex *other) {
    simplex_->neigh_[face_] = other;
  }
};


/***************************************************************************
 * A reference-counted pointer to the internal structures.
 *
 * As long as the simplicial complex hasn't been deleted,
 * it is safe to hold the Simplex, even after a removeSimplex().
 *
 * Simplices are canonically oriented according to the orientation
 * defined at construction time (or after a flipOrientation() call).
 ***************************************************************************/
public:
class Simplex : public hudson::SmallMemoryItem<Simplex> {
  private:
    friend class SimplicialComplex;
    friend class OSimplex;
    friend struct Flipper;
    friend class Star;
    friend struct OptSimplex;

    hudson::IntrusivePtr<ISimplex> s_;
    ISimplex *get() const { return s_.get(); }
    void reset(ISimplex *is) { s_ = is; }
  public:
    Simplex(const Simplex& other): s_(other.s_) { }

    Simplex(ISimplex *s): s_(s) { }

    static const unsigned dimension = d; 
    static unsigned dim() { return d; } 
    size_t hash() const { return s_->hash(); }

    std::string toString() const { return s_->toString(); }

    typedef typename ISimplex::vertex_iterator vertex_iterator;
    typedef vertex_iterator const_iterator;
    typedef vertex_iterator iterator;
    vertex_iterator begin() const { return s_->begin_vertices(); }
    vertex_iterator end() const { return s_->end_vertices(); }
    Vertex *operator[] (unsigned i) const { return s_->get(i); } 
    Vertex *get(unsigned i) const { return s_->get(i); } 
    unsigned indexOf(const Vertex *v) const { return s_->indexOf(v); } 
    bool has(const Vertex *v) const { return s_->has(v); } 
    template <class BinaryPredicate>
    bool has(const Vertex *v, const BinaryPredicate& pred) const {
      return s_->has(v,pred);
    }

    bool operator== (const Simplex& other) const { return s_ == other.s_; }
    bool operator!= (const Simplex& other) const { return s_ != other.s_; }
};

/***************************************************************************
 * An oriented handle on an Simplex.
 *
 * While simplices are canonically oriented, OSimplex allows
 * holding the simplex in some different orientation.  Use
 * the switch, getFace, and getNeighbour operations to navigate.
 *
 * Note that this class is somewhat heavyweight: in a performance-critical
 * setting, it may be better to avoid OSimplex.
 *
 * Equality and hash functions are specifically not defined on this class:
 * users typically want equality to be simplex equality regardless of
 * orientation and would be shocked to have it be based on the particular
 * orientation.  I force them to make their wishes explicit by not making
 * operator== available.
 ***************************************************************************/
class OSimplex {
  private:
    friend class SimplicialComplex;
    friend class Face;
    friend struct Star;
    friend struct Flipper;

    boost::array<Vertex*, d+1> t_;
    Simplex s_;

    OSimplex(ISimplex *s) : t_(s->toArray()), s_(s) { }

    OSimplex(ISimplex *s, unsigned i) : t_(s->toArray()), s_(s) {
      t_.multiSwapInPlace(i); // this is d-i swaps
      if (d > 1 && (((d-i) % 2) == 1)) {
        t_.swapInPlace(0);
      }
    }

    OSimplex(ISimplex *s, const boost::array<Vertex*, d+1>& verts)
      : t_(verts), s_(s) { }

    // Set the neighbour through the face that we represent.
    void setN(ISimplex *is) {
      s_->setN(s_->indexOf(t_[d]), is);
    }

    // we can't change what the OSimplex points to by getting its ISimplex,
    // so it's const.
    ISimplex *get() const { return s_.get(); }

    void reset(ISimplex *is) { s_.reset(is); }
  public:

    static const unsigned dim = d;

    bool isInternal() const { return t_[d] != NULL; }

    bool has(const Vertex *v) const { return t_.has(const_cast<Vertex*>(v)); }
    template <class BinaryPredicate>
    bool has(const Vertex *v, const BinaryPredicate& pred ) const {
      return t_.has(const_cast<Vertex*>(v), pred);
    }

#if 0
    // dangerous; don't use these.
    bool operator==(OSimplex s) const { return t == s.t; }
    uint32_t hash() const { return t.hash(); }
#endif

    Vertex *getUnchecked(unsigned i) const { return t_[i]; }
    Vertex *get(unsigned i) const {
      Vertex *v = getUnchecked(i);
      assert(v);
      return v;
    }
    Vertex *operator[](unsigned i) const { return get(i); }
    const boost::array<Vertex*, d+1>& toArray() const { return t_; }

    const Simplex& toSimplex() const { return s_; }
    Simplex& toSimplex() { return s_; }

    operator const Simplex& () const { return s_; }
    operator Simplex& () { return s_; }

    std::string toString() const {
      string s = "<";
      VertexPrinter print;
      BOOST_FOREACH(const Vertex *v, t_) {
        if (v) {
          s += " ";
          s += print(v);
        } else {
          s += " (null)";
        }
      }
      s += " >";
      return s;
    }
};


/************************************************************
 * A face is a (d-1) simplex.
 *
 * Faces are canonically ordered (the vertices are sorted by pointer value).
 * Thus, a face can be derived from either simplex that contains it.
 *
 * Note that this makes faces somewhat heavyweight.
 ************************************************************/
public:
class Face {
  private:
    boost::array<Vertex*, d> t;
    void init(const boost::array<Vertex*, d+1>& s, unsigned i) {
      for(unsigned j = 0; j < i; ++j) {
        t[j] = s[j];
      }
      for(unsigned j = i+1; j <= d; ++j) {
        t[j-1] = s[j];
      }
      hudson::sort(t);
    }

    // Friend for the ISimplex constructor.
    friend class SimplicialComplex;

    Face(const ISimplex *s, unsigned i) { init(s->toArray(), i); }

  public:
    Face(const OSimplex& os): t(os.toArray()) { t.sort(); }

    Face(const Simplex& s, unsigned i) { init(s.toArray(), i); }

    bool operator==(const Face& f) const { return t == f.t; }
    bool operator!=(const Face& f) const { return !(t == f.t); }

    size_t hash() const {
      return hudson::hashArray(t, hudson::hash_by_cast<Vertex*>());
    }

    std::string toString() const {
      string s = "<";
      VertexPrinter print;
      BOOST_FOREACH(const Vertex *v, t) {
        if (v) {
          s += " ";
          s += print(v);
        } else {
          s += " (null)";
        }
      }
      s += " >";
      return s;
    }
};

/************************************************************
 * A face2 is a (d-2) simplex.
 *
 * This is used in computing a cavity, where we know one of
 * the vertices (the apex of the star), to reduce the number
 * of things we need to sort and hash.
 ************************************************************/
public:
class Face2 {
  private:
    boost::array<Vertex*, d - 1> t;
    void init(const boost::array<Vertex*, d+1>& s, unsigned i1, unsigned i2) {
      assert(i1 != i2);
      assert(i1 <= d);
      assert(i2 <= d);
      if (i1 > i2) {
        // i1 should be < i2; fix that
        unsigned tmp = i1;
        i1 = i2;
        i2 = tmp;
      }
      for(unsigned j = 0; j < i1; ++j) {
        t[j] = s[j];
      }
      for(unsigned j = i1+1; j < i2; ++j) {
        t[j-1] = s[j];
      }
      for(unsigned j = i2+1; j <= d; ++j) {
        t[j-2] = s[j];
      }
      hudson::sort(t);
    }

    friend class SimplicialComplex;

    Face2(const ISimplex *s, unsigned i1, unsigned i2) {
      init(s->toArray(), i1, i2);
    }

  public:
    bool operator==(const Face2& f) const { return t == f.t; }
    bool operator!=(const Face2& f) const { return !(t == f.t); }

    size_t hash() const {
      return hudson::hashArray(t, hudson::hash_by_cast<Vertex*>());
    }

    std::string toString() const {
      string s = "<";
      VertexPrinter print;
      BOOST_FOREACH(const Vertex *v, t) {
        if (v) {
          s += " ";
          s += print(v);
        } else {
          s += " (null)";
        }
      }
      s += " >";
      return s;
    }
};

/***************************************************************************
 * The underlying data structure is a ref-counted link structure.
 * It is rooted at the handle_.
 ***************************************************************************/
private:
Simplex handle_;

/***************************************************************************
 * \name Internal modification routines.
 *
 * To add simplices, we need to determine the neighbourhood information.
 * We do that using a map from each face to its one or two neighbours;
 * this is identical to the face-based simplicial complex, but we only hold
 * this temporarily.  Of course, that also means we need to build it up
 * every time we need it.
 *
 * Remember to set the handle_ after adding simplices (by now, the handle is
 * probably invalid because we probably just deleted some simplices).
 *
 * The hashmap we use is my own, and we point out to it that the destructors
 * need not be called.  This is also why I use dumb pointers rather than
 * smart ones in the pair; I know the simplices are owned by someone else,
 * so I don't need to refcount this usage.
 *
 * To remove a simplex, we merely unlink it from its neighbours, mark
 * it dead, and decrement the reference count.
 ***************************************************************************/
typedef std::pair<IFace,IFace> FaceNeigh;
typedef hudson::FastHashMap<Face, FaceNeigh, hudson::hash_by_call<Face>, false> FaceHashMap;
typedef hudson::FastHashMap<Face2, FaceNeigh, hudson::hash_by_call<Face2>, false> Face2HashMap;
private:

/* When the second face of the neighbour appears, immediately link the two.
 * Note that the null face won't appear; however:
 * (1) neighbours of the cavity always have an new neighbour in the star
 * (2) the star simplices default to having a null neighbour
 */
template <class Map, class Face>
static void addFaceToFaceMap(Map& facemap, const Face& f, const IFace& s) {
  dprintf("  adding face %s => %s\n", f.toString().c_str(), s.toString().c_str());
  FaceNeigh& p = facemap[f];
  if (!p.first) {
    p.first = s;
    assert(!p.second);
    assert(p.second != s.simplex_);
  } else {
#ifndef NDEBUG
    // assert(p.second == NULL); // but do it with a printout
    if (p.second) {
      struct { string operator()(const IFace& s) const { return s.toString(); } } printer;
      string V = printer(s);
      string v1 = printer(p.first);
      string v2 = printer(p.second);
      printf("adding %s => %s\n", f.toString().c_str(), V.c_str());
      printf("  face %s conflicts: %s and %s.\n",
          f.toString().c_str(), v1.c_str(), v2.c_str());
      assert(0 && "there is already a simplex here");
    }
#endif
    assert(p.first != s.simplex_);
    p.second = s;

    /* Now, link. */
    p.first.relink(p.second);
    p.second.relink(p.first);
  }
}

static void addSimplexToFaceMap(FaceHashMap& facemap, ISimplex *s) {
  for(unsigned i = 0; i <= d; i++) {
    IFace iF(s, i);
    Face f (iF.simplex_, iF.face_);
    addFaceToFaceMap(facemap, f, iF);
  }
}



/***************************************************************************
 * Put the simplex into the complex.
 * We assume it's been linked already.
 ***************************************************************************/
private:
void addSimplex(ISimplex *s) {
  dprintf("  committing simplex %s\n", s->toString().c_str());
  assert(s->getStatus() == ISimplex::BORN);

  s->setStatus(ISimplex::INSERTED);
  s->refcount()++;
}

/***************************************************************************
 * Remove a single simplex.
 * Warning: this may disconnect the complex, and it may invalidate the
 * handle_.  Therefore, it is private.
 ***************************************************************************/
private:
void removeSimplex(const Simplex& s) {
  dprintf("  removing simplex %s\n", s.toString().c_str());
  assert(isMember(s));

  ISimplex *is = s.get();

  /* Unlink from the neighbours. */
  BOOST_FOREACH(ISimplex *neigh, is->neigh_) {
    if (neigh == NULL) { continue; }
    unsigned back = neigh->indexOf(is);
    neigh->neigh_[back] = NULL;
  }
  // Update the refcount.
  removeSimplexNoUnlink(is);
}

void removeSimplexNoUnlink(ISimplex *is) {
  assert(is->getStatus() == ISimplex::INSERTED);

  // Someone's holding on to the Simplex, so there's at least one ref
  // there, plus we're just now removing the complex's link, for two.
  assert(is->refcount() >= 2);
  is->refcount()--;

  is->setStatus(ISimplex::DEAD);
}

/***************************************************************************
 * \name Initialization
 *
 * There are two constructors: one creates a single simplex, which can then
 * be turned into any other manifold simplicial complex using replaceCavity
 * operations.
 *
 * The other constructor creates an entire complex.  This complex is not
 * guaranteed to be manifold or orientable.  Use at your own risk (or benefit
 * -- you can't create a torus with the first constructor).
 *
 * Removing a simplex requires making sure that
 ***************************************************************************/

/****************************************
 * Constructor: create a single simplex.
 ****************************************/
public:
SimplicialComplex(const vector<Vertex*>& verts) : handle_(new ISimplex(verts), verts) {
  assert(verts.size() == d + 1); // and no more
  dprintf("creating the simplex\n");
  addSimplex(handle_.get());
  dprintf("done\n");
}

/****************************************
 * Constructor: Construct a set of elements.
 *
 * Each element is made up of a tuple of vertices, but is named by a tuple of
 * indices.  If zeroindex is non-zero, then the indices need not be zero-based.
 * The orientations are as per the elements vector.  We do not check for
 * consistency (this allows creating an unorientable surface, which may or may
 * not be a problem).
 ****************************************/
public:
SimplicialComplex(const std::vector<Vertex*>& verts,
                  const std::vector<std::vector<size_t> >& elements, ssize_t zeroindex = 0)
: handle_(new ISimplex(verts, elements[0], zeroindex))
{
  assert(elements.size() > 0);
  dprintf("creating a complex\n");

  FaceHashMap facemap;
  vector<ISimplex*> simplices;
  simplices.reserve(elements.size());
  facemap.reserve(elements.size() * (d+1));

  /* Element 0 is already in handle_, so skip it in the loop. */
  simplices.push_back(handle_.get());
  addSimplexToFaceMap(facemap, handle_.get());
  size_t n = elements.size();
  for(size_t i = 1; i < n; ++i) {
    ISimplex *s = new ISimplex(verts, elements[i], zeroindex);
    simplices.push_back(s);
    addSimplexToFaceMap(facemap, s);
  }

  /* Link up the simplices. */
  BOOST_FOREACH(ISimplex *is, simplices) {
    addSimplex(is);
  }
  dprintf("done\n");
}

/***************************************************************************
 * \name Dimension
 *
 * Return the topological dimension of the complex.
 ***************************************************************************/
public:
static unsigned dim() { return d; }
static const unsigned dimension = d;
/***************************************************************************
 * \name Entering the complex.
 *
 * The complex maintains a handle; this is necessary after initialization,
 * and is useful besides that time.  On insertion of new vertices, the handle
 * is kept up to date; the user will only need to call setHandle() when removing
 * simplices (namely, if we remove the handle).
 *
 * Simplices in the complex are canonically oriented.  If the orientation is
 * the "wrong" one -- if you wish to view the complex from the opposite
 * orientation -- call flipOrientation().  This makes in-memory changes, and
 * takes linear time.  Prior handles into the complex are also flipped.
 ***************************************************************************/
/***********************************
 * Check whether the given simplex appears in the complex.
 ***********************************/
public:
bool isMember(const Simplex& s) const {
  return s.get()->getStatus() == ISimplex::INSERTED;
}

public:
const Simplex& getHandle() const {
  assert(isMember(handle_));
  return handle_;
}

OSimplex orientSimplex(const Simplex& s) const {
  return OSimplex(s.get());
}

public:
void setHandle(const Simplex& s) {
  assert(isMember(s));
  handle_ = s;
}

private:
// Helper for flipOrientation
struct Flipper {
  static bool canEnter(const Simplex&) { return true; }
  static bool choose(const Simplex& s) {
    /* Note that this changes the complex in place, but
       we're not in the middle of reading s's neighbours,
       so it turns out to be safe. */
    hudson::swapInPlace(s.get()->verts_, 0);
    hudson::swapInPlace(s.get()->neigh_, 0);
    return false;
  }
};
/****************************************
 * Flip the orientation of the complex.
 * This operation is linear-time.
 ****************************************/
public:
void flipOrientation() {
  Flipper flip;
  dfsBySimplex(flip, getHandle());
}
/***************************************************************************
 * \name Inserting a vertex
 *
 * Vertex insertion happens as follows:
 * -# Compute a <i>cavity</i>: a connected set of simplices.
 * -# Create a vertex.
 * -# Call computeStar to compute the replacement of the cavity.  Inspect the
 *    star if you want; you may safely discard it.
 * -# Call replaceCavity, which commits the star to the structure.
 *
 * Steps 3-4 can be combined into one by calling the other form of insertVertex.
 *
 * The star is guaranteed to be oriented consistently with the prior orientation.
 ***************************************************************************/
/***************************************************************************
 * Represents the star of a vertex.
 *
 * The star acts like a const vector<Simplex>: it has size(), operator[], and
 * begin/end.
 *
 * Internally, we hold the apex, the simplices of the star, and their
 * (possibly-null) neighbours outside the cavity.  This is all required for
 * initially constructing the star, and is also required for linking the star
 * in during replaceCavity.
 ***************************************************************************/
public:
struct Star {
  typedef hudson::SmallList<std::pair<Simplex,IFace> > StarData;
  typedef typename std::pair<Simplex,IFace> pair;
  typedef typename StarData::const_iterator pair_it;

  Star() : apex_(NULL), size_(0) { }

  void clear() {
    star_.clear();
    apex_ = NULL;
    size_ = 0;
  }

  template <class sdata_iterator>
  class iterator_base : public boost::iterator_facade
                                < iterator_base<sdata_iterator>
                                , const Simplex
                                , boost::forward_traversal_tag
                                >
  {
    sdata_iterator it_;
    friend struct Star;
    explicit iterator_base (const sdata_iterator& it): it_(it) { }
    public:
    void increment() { ++it_; }
    bool equal(const iterator_base& other) const { return it_ == other.it_; }
    const Simplex& dereference() const { return (*it_).first; }
  };
  typedef iterator_base<typename StarData::iterator> iterator;
  typedef iterator_base<typename StarData::const_iterator> const_iterator;

  iterator begin() { return iterator(star_.begin()); }
  iterator end() { return iterator(star_.end()); }

  const_iterator begin() const { return const_iterator(star_.begin()); }
  const_iterator end() const { return const_iterator(star_.end()); }

  struct neigh_iterator : public boost::iterator_facade
                                < neigh_iterator
                                , Simplex
                                , boost::forward_traversal_tag
                                , Simplex /* not a ref: we construct it here */
                                >
  {
    pair_it it_;
    friend struct Star;
    explicit neigh_iterator (const pair_it& it): it_(it) { }
    public:
    void increment() { ++it_; }
    bool equal(const iterator& other) const { return it_ == other.it_; }
    const Simplex& dereference() const {
      const IFace& iface = (*it_).second;
      return Simplex(iface.simplex_);
    }
  };
  neigh_iterator begin_neighbours() const { return star_.begin(); }
  neigh_iterator end_neighbours() const { return star_.end(); }

  size_t size() const { return size_; }

  iterator erase(const iterator& it) {
    size_--;
    return iterator(star_.erase(it.it_));
  }

  private:
  friend class SimplicialComplex;

  /* pairs of <simplex, neighbour> */
  StarData star_;
  Vertex *apex_;
  size_t size_;

  pair_it beginPairs() const { return star_.begin(); }
  pair_it endPairs() const { return star_.end(); }
  static ISimplex *getSimplex(const pair_it& it) {
    return (*it).first.get();
  }
  static IFace& getNeighbour(const pair_it& it) {
    return (*it).second;
  }

  void setApex(Vertex *apex) {
    assert(apex_ == NULL);
    apex_ = apex;
  }

  void addBoundaryFace(IFace face) {
    // Build the new star simplex
    boost::array<Vertex*, d+1> starverts = face.simplex_->toArray();
    starverts[face.face_] = apex_;
    ISimplex *is = new ISimplex(starverts);

    // Remember its neighbour, if any.
    IFace neigh = face.getNFace();
    dprintf("  neighbour %s\n", neigh.toString().c_str());

    // Take ownership of is, so that it doesn't go away
    // when it gets used.
    Simplex s(is);
    star_.push_front(make_pair(s, neigh));
    size_++;
  }
};

struct Cavity {
  typedef hudson::SmallList<Simplex> storage_type;
  typedef typename storage_type::const_iterator iterator;
  typedef typename storage_type::const_iterator const_iterator;

  Cavity(): size_(0) { }

  iterator begin() { return cavity_.begin(); }
  iterator end() { return cavity_.end(); }
  const_iterator begin() const { return cavity_.begin(); }
  const_iterator end() const { return cavity_.end(); }

  bool empty() const { return cavity_.empty(); }
  size_t size() const { return size_; }

  void add(const Simplex& s) { cavity_.push_front(s); size_++; }
  void clear() { cavity_.clear(); size_ = 0; }

  private:
  storage_type cavity_;
  size_t size_;
};

void computeStar(const Cavity& cavity, Vertex *apex, Star& output)
{
  assert(!cavity.empty());
  output.setApex(apex);

  VisitedSet vset(this);
  BOOST_FOREACH(const Simplex& s, cavity) {
    assert(isMember(s));
    vset.insert(s.get());
  }

  BOOST_FOREACH(const Simplex& s, cavity) {
    ISimplex *is = s.get();

    for(unsigned j = 0; j <= d; ++j) {
      ISimplex *neigh = is->getN(j);
      if (!neigh || !vset.has(neigh)) {
        output.addBoundaryFace(IFace(is, j));
      }
    }
  }
}

/***************************************************************************
 * Given a cavity, and a replacement for it, perform the replacement.
 *
 * Algorithm:
 * -# remove the old links to the cavity
 * -# compute all the links between star simplices
 * -# add the in-star links
 * -# add the star-neighbour links
 ***************************************************************************/
public:
void replaceCavity(const Cavity& cavity, Star& star)
{
  dprintf("replaceCavity\n");

  // Remove the old.
  BOOST_FOREACH(const Simplex& s, cavity) {
    removeSimplexNoUnlink(s.get());
  }

  // Compute the face map and link faces.
  //
  // There are exactly d+1 faces for each star simplex.
  // Exactly one face points outside, so d point to other
  // star simplices.
  //
  // I can handle the exterior face specially, and then all the interior faces
  // share the apex.  So I use Face2, saving the sorting, hashing, and comparison
  // of one vertex from each of the cavity faces.
  Face2HashMap map (star.size() * d);

  // I merge this loop with the one that adds the simplices
  // in to the mesh.
  BOUNDED_FOREACH(typename Star::pair p, star.beginPairs(), star.endPairs()) {
    ISimplex *isimplex = p.first.get();
    IFace& exteriorNeighbour = p.second;
    unsigned apexIndex = isimplex->indexOf(star.apex_);

    /* Link the new simplex to its neighbour on the exterior of the cavity, if
     * there is a neighbour. */
    if (exteriorNeighbour) {
      exteriorNeighbour.relink(isimplex);
      IFace apexFace (isimplex, apexIndex);
      apexFace.relink(exteriorNeighbour);
    }

    /* Iterate over all faces except the face across from the
     * apex, linking through the given cavity-interior face.
     * The apex is always part of the cavity face, so we simply
     * omit it from the search structure.
     *
     * Break the loop in two and keep indices sorted to
     * let the optimizer do its work. */
    for (unsigned i = 0 ; i < apexIndex; ++i) {
      addFaceToFaceMap(map, Face2(isimplex, i, apexIndex), IFace(isimplex, i));
    }
    for (unsigned i = apexIndex + 1 ; i <= d; ++i) {
      addFaceToFaceMap(map, Face2(isimplex, apexIndex, i), IFace(isimplex, i));
    }

    /* Add the simplex to the mesh.  Note that it may not yet
       be fully linked in until the loop completes. */
    addSimplex(isimplex);
  }

  // Handle may be invalid; set it to something we know works.
  handle_ = *star.begin();
  dprintf("done\n");
}

/***************************************************************************
 * Given a star-shaped cavity, and a vertex, compute the star and install
 * it into the complex.
 ***************************************************************************/
public:
void replaceCavity(const Cavity& cavity, Vertex *v) {
  Star star;
  computeStar(cavity, v, star);
  replaceCavity(cavity, star);
}

/***************************************************************************
 * Remove a single simplex.
 * Checks that the removed simplex is not the handle.
 * To remove a simplex that is the handle, call setHandle() first.
 ***************************************************************************/
public:
void checkedRemoveSimplex(const Simplex& s) {
  assert(s != getHandle());
  removeSimplex(s);
}
/***************************************************************************
 * \name Data access.
 *
 * We can handle void data.  However, it is an error to call these functions
 * when the data is void.
 ***************************************************************************/
data_ref getDataRW(const Simplex& s) {
  return s.get()->refcount_and_user_.payload_ref();
}

data_const_ref getData(const Simplex& s) const {
  return s.get()->refcount_and_user_.payload_cref();