SimplicialHelpers.h

Go to the documentation of this file.

#ifndef SimplicialHelpers_HEADER
#define SimplicialHelpers_HEADER

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

#include "SimplicialComplex.h"
#include <iterator>
#include <vector>
#include <utilities/FastHash.h>
#include <utilities/hash.h>
#include <utilities/count_iterator.h>


/***************************************************************************
 * \file
 * Helpers to use with SimplicialComplex.
 *
 * The idiom is to have a function foo(), which uses a namespace foo_.
 * This is because we can't have templated functions defining structures
 * in their own namespace, for some reason.
 ***************************************************************************/
namespace Simplicial {


template <class Vertex> class VertexSet : public hashers::hash_set<Vertex*, hudson::hash_by_cast<Vertex*> > { };

/***************************************************************************
 * \name Debugging: Checking the shape of a complex
 *
 * The following functions and structures allow checking that the complex
 * has the expected shape.  I can only think of this being useful for
 * debugging and testing.
 ***************************************************************************/
template <class Complex, class BinaryPredicate = std::equal_to<const typename Complex::Vertex*> >
struct Element {
  typedef typename Complex::Vertex Vertex;
  typedef typename Complex::Simplex Simplex;

  private:
  vector<Vertex*> verts_;

  public:
  template <class iterator>
  Element(const iterator& begin, const iterator& end) : verts_(begin, end) { }

  Element(const vector<Vertex*>& v):verts_(v) { }

  template <class eiterator>
  Element(const vector<Vertex*>& v, eiterator begin, eiterator end, int offset) {
    while(begin != end) {
      verts_.push_back(v[*begin - offset]);
      ++begin;
    }
  }

  bool matches(const Simplex& s) const {
    BinaryPredicate pred;
    for (size_t i = 0; i < verts_.size(); ++i) {
      if (!s.has(verts_[i], pred)) {
        return false;
      }
    }
    return true;
  }
};

namespace checkComplex_ {
  template <class Complex>
    struct Checker {
      typedef Simplicial::Element<Complex> Element;
      typedef typename Complex::Simplex Simplex;

      const std::vector<Element>& correct_;
      size_t nseen_;

      Checker(const std::vector<Element>& c): correct_(c), nseen_(0) { }
      static bool canEnter(const Simplex&) { return true; }
      bool choose(const Simplex& s) {
        nseen_++;
        for(size_t i = 0; i < correct_.size(); ++i) {
          if (correct_[i].matches(s)) return false;
        }
        return true;
      }
    };
}
template <class Complex>
void
checkComplex(const Complex& cplx, const vector<Element<Complex> >& correct) {
  checkComplex_::Checker<Complex> checker(correct);
  typename Complex::OptSimplex bad
    = cplx.findSimplex(checker, cplx.getHandle());
  if (bad.isSome()) {
    printf("first bad simplex: %s\n", (*bad).toString().c_str());
    assert(0 && "bad simplex");
  }
  assert(checker.nseen_ == correct.size());
  static unsigned testnum = 0;
  printf("check #%u passed\n", testnum++);
}

template <class Complex>
void
checkComplex(const Complex& c,
    const vector<typename Complex::Vertex*>& vertices,
    const vector<vector<size_t> >& elements,
    ssize_t offset = 0) {
  typedef Element<Complex> Element;
  vector<Element> check;
  for(size_t i = 0; i < elements.size(); ++i) {
    check.push_back(Element(vertices, elements[i].begin(), elements[i].end(), offset));
  }
  checkComplex(c, check);
}
/***************************************************************************
 * \name Debugging: printing.
 *
 * Debugging support for printing a complex.  This uses the Simplex::toString()
 * call, which is not guaranteed to be pretty.
 ***************************************************************************/
namespace printComplex_ {
  template <class Complex>
  struct printer {
    vector<string> strs;
    static bool canEnter(const typename Complex::Simplex&) { return true; }
    bool choose(const typename Complex::Simplex& os) {
      strs.push_back(os.toString());
      return false;
    }
  };
}
template <class Complex>
void printComplex(const Complex& cplx) {
  printComplex_::printer<Complex> p;
  cplx.dfsBySimplex(p, cplx.getHandle());
  std::sort(p.strs.begin(), p.strs.end());
  printf("%u elements:\n", (unsigned)p.strs.size());
  for(size_t i = 0; i < p.strs.size(); ++i) {
    printf("  %s\n", p.strs[i].c_str());
  }
}
/***************************************************************************
 * \name Collecting the star of a vertex.
 *
 * Collect all the simplices around a named vertex (its star).
 ***************************************************************************/
namespace collectStar_ {
  template <class Complex>
  struct finder {
    typedef typename Complex::Simplex Simplex;
    typedef typename Complex::Vertex Vertex;
    const Vertex *v;
    finder(const Vertex *v) : v(v) { }
    static bool canEnter(const Simplex&) { return true; }
    bool choose(const Simplex& os) const { return os.has(v); }
  };

  template <class Complex, class output_iterator>
  struct collector {
    typedef typename Complex::Simplex Simplex;
    typedef typename Complex::Vertex Vertex;
    const Vertex *v;
    output_iterator out;
    collector(const Vertex *v, const output_iterator& out) : v(v), out(out) { }
    bool canEnter(const Simplex& os) const { return os.has(v); }
    bool choose(const Simplex& os) { *out = os; ++out; return false; }
  };
}

template <class Complex, class output_iterator>
output_iterator collectStar(const Complex& cplx,
                 const typename Complex::Vertex *v,
                 const typename Complex::Simplex& seed,
                 const output_iterator& out)
{
  assert (seed.has(v));
  collectStar_::collector<Complex, output_iterator> c(v, out);
  cplx.dfsBySimplex(c, seed);
  return c.out;
}

template <class Complex, class output_iterator>
output_iterator slowlyCollectStar(const Complex& cplx,
                                  const typename Complex::Vertex *v,
                                  output_iterator out)
{
  typedef typename Complex::Simplex Simplex;
  typedef typename Complex::OptSimplex OptSimplex;

  collectStar_::finder<Complex> f(v);
  OptSimplex start = cplx.findSimplex(f, cplx.getHandle());
  if (start.isSome()) {
    out = collectStar(cplx, v, *start, out);
  }
  return out;
}
/***************************************************************************
 * \name Collect all the simplices in the mesh.
 ***************************************************************************/

namespace collectSimplices_ {
  template <class Complex, class output_iterator>
  struct C {
    typedef typename Complex::Simplex Simplex;
    output_iterator out;
    C(output_iterator oo) : out(oo) { }
    static bool canEnter(const Simplex&) { return true; }
    bool choose(const Simplex& s) { *out = s; ++out; return false; }
  };
}

template <class Complex, class output_iterator>
output_iterator collectSimplices(const Complex& cplx, output_iterator out) {
  collectSimplices_::C<Complex, output_iterator> c(out);
  cplx.dfsBySimplex(c, cplx.getHandle());
  return c.out;
}

template <class Complex>
void collectSimplices(const Complex& cplx, std::vector<typename Complex::Simplex>& out) {
  collectSimplices(cplx, std::back_inserter(out));
}

template <class Complex>
size_t countSimplices(const Complex& cplx) {
  return collectSimplices(cplx, hudson::count_iterator<typename Complex::Simplex>()).count;
}


/***************************************************************************
 * \name Collect all the vertices in the mesh.
 *
 * We use a set for this task.  At the moment, we also iterate over the set,
 * which has the unfortunate effect of essentially randomizing the order of
 * vertices.
 ***************************************************************************/
namespace collectVertices_ {
  template <class Simplex, class Vertex>
  void addVertices(const Simplex& os, VertexSet<Vertex>& set) {
    for (size_t i = 0; i <= os.dimension; ++i) {
      set.insert(os[i]);
    }
  }

  template <class Complex>
  struct C {
    typedef typename Complex::Simplex Simplex;
    typedef typename Complex::Vertex Vertex;
    VertexSet<Vertex> set;
    C() { }
    C(const Complex& cplx) {
      cplx.dfsBySimplex(*this, cplx.getHandle());
    }
    static bool canEnter(const Simplex&) { return true; }
    bool choose(const Simplex& os) { addVertices(os, set); return false; }
  };
}

template <class iterator/*of Simplex */, class Vertex>
void collectVertices(iterator begin, iterator end,
    std::vector<Vertex*>& out) {
  VertexSet<Vertex> set;
  while (begin != end) {
    collectVertices_::addVertices(*begin, set); ++begin;
  }
  std::copy(set.begin(), set.end(), std::back_inserter(out));
}

template <class Complex, class output_iterator>
output_iterator collectVertices(const Complex& cplx,
    output_iterator out) {
  collectVertices_::C<Complex> c(cplx);
  return std::copy(c.set.begin(), c.set.end(), out);
}

template <class Complex>
void collectVertices(const Complex& cplx,
    std::vector<typename Complex::Vertex*>& out) {
  collectVertices(cplx, std::back_inserter(out));
}

template <class Complex>
size_t countVertices(const Complex& cplx) {
  collectVertices_::C<Complex> c(cplx);
  return c.set.size();
}


/***************************************************************************
 * \name Serialization code.
 *
 * This is the question of turning our complex into something we can store
 * on disk.  We translate it here into an array of vertices and an array
 * of indices (namely, something from which we could create a new complex).
 *
 * The implementation is quite fast given the output format of a vector
 * of vectors of integers: I only perform one pass over the complex.
 * However, a vector of vectors is agonizingly slow due to the memory
 * allocation cost.
 ***************************************************************************/
template <class Complex>
struct Serialize_ {
  typedef typename Complex::Simplex Simplex;
  typedef typename Complex::Vertex Vertex;
  typedef hashers::hash_map<Vertex*, unsigned, hudson::hash_by_cast<Vertex*> > VMap;

  vector<Vertex*>& verts;
  vector<vector<size_t> >& elements;
  VMap vmap;
  Serialize_(vector<Vertex*>& v, vector<vector<size_t> >& e)
    : verts(v), elements(e) { }

  static bool canEnter(const Simplex&) { return true; }

  bool choose(const Simplex& os) {
    elements.push_back(vector<size_t>());
    vector<size_t>& elt = elements.back();
    elt.resize(Complex::dimension+1);
    for (size_t j = 0; j <= Complex::dimension; ++j) {
      // Look up the vertex in the table, or insert it if it isn't there.
      size_t id;
      size_t n = vmap.size();
      pair<typename VMap::iterator, bool> p = vmap.insert(make_pair(os[j], n));
      if(p.second) {
        // Succeeded => this is a new vertex with id = n.
        verts.push_back(os[j]);
        id = n;
      } else {
        // Failed => we can read out the old ID.
        id = (*(p.first)).second;
      }
      // Now we know the ID, so we can write to the element array.
      elt[j] = id;
    }
    return false;
  }
};

template <class Complex>
void serialize(const Complex& cplx,
    std::vector<typename Complex::Vertex*>& verts,
    std::vector<std::vector<size_t> >& elements)
{
  Serialize_<Complex> s(verts, elements);
  cplx.dfsBySimplex(s, cplx.getHandle());
}

} // namespace Simplicial
#endif

---