cn94.cpp

Go to the documentation of this file.

#include <limits>
#include <iostream>
#include <iterator>
#include <vector>
#include "cn94.hpp"

using namespace Arak;
using namespace Arak::Geometry;

const char* CN94Proposal::moveNames[numMoveTypes] = {
  "interior triangle birth",
  "interior triangle death",
  "recolor quadrilateral",
  "move interior vertex",
  "move vertex along boundary",
  "move vertex past corner",
  "interior vertex birth",
  "interior vertex death",
  "boundary triangle birth",
  "boundary triangle death",
  "corner cut birth",
  "corner cut death",
  "invalid move"
};

Point CN94Proposal::sampleInCastingBox(const Point& p,
               Arak::Util::Random& random) const {
  double x = CGAL::to_double(p.x());
  double y = CGAL::to_double(p.y());
  return Point(random.uniform(x - castingBoxRadius,
            x + castingBoxRadius),
         random.uniform(y - castingBoxRadius,
            y + castingBoxRadius));
}

bool CN94Proposal::inCastingBox(const Point& p, const Point& q) const {
  return ((fabs(CGAL::to_double(p.x() - q.x())) <= castingBoxRadius) &&
    (fabs(CGAL::to_double(p.y() - q.y())) <= castingBoxRadius));
}


void CN94Proposal::InteriorTriangleBirth::reset(const Point& u, const Point& v, const Point& w) {
  this->u = u;
  this->v = v;
  this->w = w;
  vh = Coloring::VertexHandle();
}

void CN94Proposal::InteriorTriangleBirth::execute(Coloring& c) {
  vh = c.newInteriorTriangle(u, v, w, Coloring::DO_WITHOUT_TEST_TAG());
  assert(vh.valid());
}
  
void CN94Proposal::InteriorTriangleBirth::undo(Coloring& c) {
  assert(c.deleteIntTriangle(vh, Coloring::DO_WITHOUT_TEST_TAG()));
}
  

void CN94Proposal::InteriorTriangleDeath::reset(Coloring::VertexHandle vertex,
            bool isShort) {
  vh = vertex;
  this->isShort = isShort;
  u = vertex->point();
  const Coloring::VertexHandle intVertex = vertex->getNextVertex();
  const Coloring::VertexHandle bdVertex = intVertex->getNextVertex();
  v = bdVertex->point();
  w = intVertex->point();
}

void CN94Proposal::InteriorTriangleDeath::execute(Coloring& c) {
  assert(vh.valid());
  assert(c.deleteIntTriangle(vh, Coloring::DO_WITHOUT_TEST_TAG()));
}
  
void CN94Proposal::InteriorTriangleDeath::undo(Coloring& c) {
  assert(!vh.valid());
  vh = c.newInteriorTriangle(u, v, w, Coloring::DO_WITHOUT_TEST_TAG());
  assert(vh.valid());
}
  

void CN94Proposal::Recolor::reset(Coloring::IntEdgeHandle edge1, 
          Coloring::IntEdgeHandle edge2,
          bool joinSourceToSource,
          bool convex) {
  this->edge1 = edge1;
  this->edge2 = edge2;
  this->joinSourceToSource = joinSourceToSource;
  this->convex = convex;
}

void CN94Proposal::Recolor::execute(Coloring& c) {
  Coloring::VertexHandle a = edge1->getPrevVertex();
  Coloring::VertexHandle b = edge1->getNextVertex();
  assert(c.recolorQuadrilateral(edge1, edge2, 
        joinSourceToSource, Coloring::DO_WITHOUT_TEST_TAG()));
  // Update joinSourceToSource for a potential undo.
  joinSourceToSource = 
    (((a == edge1->getPrevVertex()) && (b == edge2->getPrevVertex())) ||
     ((b == edge1->getPrevVertex()) && (a == edge2->getPrevVertex())) ||
     ((a == edge1->getNextVertex()) && (b == edge2->getNextVertex())) ||
     ((b == edge1->getNextVertex()) && (a == edge2->getNextVertex())));
}

void CN94Proposal::Recolor::undo(Coloring& c) {
  execute(c);
}


void CN94Proposal::MoveBdVertexPastCorner::reset(Coloring::BdEdgeHandle edge, 
             const Point& point) { 
  this->edge = edge;
  this->oldLoc = (edge->getPrevVertex()->type() == Coloring::Vertex::CORNER) ?
    edge->getNextVertex()->point() : edge->getPrevVertex()->point();
  this->newLoc = point;
}

void CN94Proposal::MoveBdVertexPastCorner::execute(Coloring& c) {  
  assert(c.moveBdVertexPastCorner(edge, newLoc, Coloring::DO_WITHOUT_TEST_TAG()));
}

void CN94Proposal::MoveBdVertexPastCorner::undo(Coloring& c) { 
  assert(edge.valid());
  assert(c.moveBdVertexPastCorner(edge, oldLoc, Coloring::DO_WITHOUT_TEST_TAG()));
}


void CN94Proposal::InteriorVertexBirth::reset(Coloring::IntEdgeHandle edge, 
                const Point& point) {
  this->prevVertex = edge->getPrevVertex();
  this->nextVertex = edge->getNextVertex();
  this->point = point;
}

void CN94Proposal::InteriorVertexBirth::execute(Coloring& c) {
  Coloring::IntEdgeHandle edge = prevVertex->getNextIntEdge();
  assert(c.splitEdge(edge, point, Coloring::DO_WITHOUT_TEST_TAG()));
}

void CN94Proposal::InteriorVertexBirth::undo(Coloring& c) {
  Coloring::VertexHandle newVertex = prevVertex->getNextVertex();
  assert(c.deleteVertex(newVertex, Coloring::DO_WITHOUT_TEST_TAG()));
}


void CN94Proposal::InteriorVertexDeath::reset(Coloring::VertexHandle vertex) {
  this->prevVertex = vertex->getPrevVertex();
  this->point = vertex->point();
  Coloring::VertexHandle nextVertex = vertex->getNextVertex();
  newEdgeLength = sqrt(CGAL::to_double(squared_distance(prevVertex->point(), 
              nextVertex->point())));
}

void CN94Proposal::InteriorVertexDeath::execute(Coloring& c) {
  assert(prevVertex.valid());
  Coloring::VertexHandle vertex = prevVertex->getNextVertex();
  assert(c.deleteVertex(vertex, Coloring::DO_WITHOUT_TEST_TAG()));
}

void CN94Proposal::InteriorVertexDeath::undo(Coloring& c) {
  Coloring::IntEdgeHandle edge = prevVertex->getNextIntEdge();
  assert(c.splitEdge(edge, point, Coloring::DO_WITHOUT_TEST_TAG()));
}


void CN94Proposal::MoveInteriorVertex::reset(Coloring::VertexHandle vertex, 
               const Point& point) {
  this->vertex = vertex;
  oldLoc = vertex->point();
  newLoc = point;
}

void CN94Proposal::MoveInteriorVertex::execute(Coloring& c) {  
  assert(c.moveIntVertex(vertex, newLoc, Coloring::DO_WITHOUT_TEST_TAG()));
}

void CN94Proposal::MoveInteriorVertex::undo(Coloring& c) { 
  assert(c.moveIntVertex(vertex, oldLoc, Coloring::DO_WITHOUT_TEST_TAG()));
}


void CN94Proposal::MoveVertexAlongBdEdge::reset(Coloring::VertexHandle vertex, const Point& point) { 
  this->vertex = vertex;
  this->oldLoc = vertex->point();
  this->newLoc = point;
}

void CN94Proposal::MoveVertexAlongBdEdge::execute(Coloring& c) {  
  assert(c.moveVertexAlongBd(vertex, newLoc, Coloring::DO_WITHOUT_TEST_TAG()));
}

void CN94Proposal::MoveVertexAlongBdEdge::undo(Coloring& c) { 
  assert(c.moveVertexAlongBd(vertex, oldLoc, Coloring::DO_WITHOUT_TEST_TAG()));
}


void CN94Proposal::BoundaryTriangleBirth::reset(Coloring::BdEdgeHandle e,
            const Point& u, 
            const Point& v, 
            const Point& w) {
  this->edge = e;
  this->u = u;
  this->v = v;
  this->w = w;
  vh = Coloring::VertexHandle();
}

void CN94Proposal::BoundaryTriangleBirth::execute(Coloring& c) {
  assert(edge.valid());
  assert(c.newBoundaryTriangle(edge, u, v, w, Coloring::TEST_IF_VALID_TAG()));
  vh = c.newBoundaryTriangle(edge, u, v, w, Coloring::DO_WITHOUT_TEST_TAG());
  assert(vh.valid());
  edge = Coloring::BdEdgeHandle();
}
  
void CN94Proposal::BoundaryTriangleBirth::undo(Coloring& c) {
  assert(vh.valid());
  assert(c.deleteBdTriangle(vh, Coloring::DO_WITHOUT_TEST_TAG()));
}
  

void CN94Proposal::BoundaryTriangleDeath::reset(Coloring::VertexHandle vertex) {
  assert(vertex->type() == Coloring::Vertex::BOUNDARY);
  Coloring::VertexHandle u = vertex;
  assert(u->hasNextIntEdge());
  Coloring::VertexHandle w = u->getNextVertex();
  assert(w->type() == Coloring::Vertex::INTERIOR);
  Coloring::VertexHandle v = w->getNextVertex();
  assert(v->type() == Coloring::Vertex::BOUNDARY);
  vh = vertex;
  wp = *w;
  up = *u;
  vp = *v;
  // Compute the previous and next boundary vertices.
  if (u->getNextBdEdge()->getNextVertex() == v) {
    prevBoundaryVertex = u->getPrevBdEdge()->getPrevVertex();
    nextBoundaryVertex = v->getNextBdEdge()->getNextVertex();
  } else if (u->getPrevBdEdge()->getPrevVertex() == v) {
    prevBoundaryVertex = v->getPrevBdEdge()->getPrevVertex();
    nextBoundaryVertex = u->getNextBdEdge()->getNextVertex();
  } else
    assert(false);
}

void CN94Proposal::BoundaryTriangleDeath::execute(Coloring& c) {
  assert(vh.valid());
  assert(c.deleteBdTriangle(vh, Coloring::DO_WITHOUT_TEST_TAG()));
}
  
void CN94Proposal::BoundaryTriangleDeath::undo(Coloring& c) {
  assert(!vh.valid());
  Coloring::BdEdgeHandle e = prevBoundaryVertex->getNextBdEdge();
  vh = c.newBoundaryTriangle(e, up, vp, wp, Coloring::DO_WITHOUT_TEST_TAG());
  assert(vh.valid());
}
  

void CN94Proposal::CornerCutBirth::reset(Coloring::VertexHandle corner,
           const Point& u, 
           const Point& v) {
  assert(corner->type() == Coloring::Vertex::CORNER);
#ifdef EXACT_GEOM_CXN
  assert(corner->getPrevBdEdge()->has_on(u));
  assert(corner->getNextBdEdge()->has_on(v));
#endif
  this->corner = corner;
  this->u = u;
  this->v = v;
  vh = Coloring::VertexHandle();
}

void CN94Proposal::CornerCutBirth::execute(Coloring& c) {
  assert(corner.valid());
  vh = c.newCornerTriangle(corner, u, v, Coloring::DO_WITHOUT_TEST_TAG());
  assert(vh.valid());
}
  
void CN94Proposal::CornerCutBirth::undo(Coloring& c) {
  assert(vh.valid());
  assert(c.deleteCornerTriangle(vh, Coloring::DO_WITHOUT_TEST_TAG()));
}
  

void CN94Proposal::CornerCutDeath::reset(Coloring::VertexHandle corner) {
  assert(corner->type() == Coloring::Vertex::CORNER);
  Coloring::VertexHandle prevBdVertex = 
    corner->getPrevBdEdge()->getPrevVertex();
  Coloring::VertexHandle nextBdVertex = 
    corner->getNextBdEdge()->getNextVertex();
  assert(prevBdVertex->type() == Coloring::Vertex::BOUNDARY);
  assert(nextBdVertex->type() == Coloring::Vertex::BOUNDARY);
  if (prevBdVertex->hasNextIntEdge()) {
    vh = prevBdVertex;
    assert(vh->getNextVertex() == nextBdVertex);
  } else if (nextBdVertex->hasNextIntEdge()) {
    vh = nextBdVertex;
    assert(vh->getNextVertex() == prevBdVertex);
  } else assert(false);
  this->corner = corner;
  this->u = *prevBdVertex;
  this->v = *nextBdVertex;
}

void CN94Proposal::CornerCutDeath::execute(Coloring& c) {
  assert(vh.valid());
  assert(c.deleteCornerTriangle(vh, Coloring::DO_WITHOUT_TEST_TAG()));
}
  
void CN94Proposal::CornerCutDeath::undo(Coloring& c) {
  assert(!vh.valid());
  vh = c.newCornerTriangle(corner, u, v, Coloring::DO_WITHOUT_TEST_TAG());
  assert(vh.valid());
}
  

CN94Proposal::CN94Proposal(const Arak::Util::PropertyMap& props) :
  itd(*this), 
  ivd(*this),
  curMoveType(INVALID_MOVE),
  castingBoxRadius(1.0)
{ 
  using namespace Arak::Util;
  if (hasp(props, "arak.cn94.casting_box_side")) {
    assert(parse(getp(props, "arak.cn94.casting_box_side"), castingBoxRadius));
    assert(castingBoxRadius > 0.0);
    castingBoxRadius /= 2.0;
  } else if (hasp(props, "arak.scale")) {
    double scale;
    assert(parse(getp(props, "arak.scale"), scale));
    // Smart initialization (used in Nicholl's implementation)
    castingBoxRadius = 1.0 / (4.0 * scale * sqrt(M_PI));
  }
  // std::cerr << "Using casting box radius: " << castingBoxRadius << std::endl;
  // Initialize the parameters of the proposal.
  assert(parse(getp(props, "arak.cn94.zeta.bd_not_mr"), ZETA_BD_NOT_MR));
  assert(parse(getp(props, "arak.cn94.zeta.move_not_recolor"), 
         ZETA_MOVE_NOT_RECOLOR));
  assert(parse(getp(props, "arak.cn94.zeta.birth_not_death"), 
         ZETA_BIRTH_NOT_DEATH));
  assert(parse(getp(props, "arak.cn94.zeta.int_not_boundary"), 
         ZETA_INT_NOT_BOUNDARY));
  assert(parse(getp(props, "arak.cn94.zeta.bound_not_corner"), 
         ZETA_BOUND_NOT_CORNER));
  assert(parse(getp(props, "arak.cn94.zeta.ib_vertex_not_triangle"), 
         ZETA_IB_VERTEX_NOT_TRIANGLE));
  // Initialize the statistics.
  for (int i = 0; i < numMoveTypes; i++) {
    numProposed[i] = 0;
    numAccepted[i] = 0;
  }
}

CN94Proposal::~CN94Proposal() { }

ColoringMove& CN94Proposal::move() {
  switch (curMoveType) {
  case INT_TRIANGLE_BIRTH:         return itb;
  case INT_TRIANGLE_DEATH:         return itd;
  case RECOLOR:                    return rec;
  case MOVE_INT_VERTEX:            return miv;
  case MOVE_BD_VERTEX_ALONG_BD:    return mbb;
  case MOVE_BD_VERTEX_PAST_CORNER: return mbc;
  case INT_VERTEX_BIRTH:           return ivb;
  case INT_VERTEX_DEATH:           return ivd;
  case BD_TRIANGLE_BIRTH:          return btb;
  case BD_TRIANGLE_DEATH:          return btd;
  case CORNER_CUT_BIRTH:           return ccb;
  case CORNER_CUT_DEATH:           return ccd;
  case INVALID_MOVE:
    return (ColoringMove&)(NullColoringMove::instance);
  default:
    assert(false);
  }
}

void CN94Proposal::sample(const Coloring& state, 
        Arak::Util::Random& random,
        bool reversible) {
  curMoveType = UNSPECIFIED_MOVE;
  if (random.uniform() < ZETA_BD_NOT_MR) {
    // Birth/death group
    if (random.uniform() < ZETA_BIRTH_NOT_DEATH) {
      // Birth group
      if (random.uniform() < ZETA_INT_NOT_BOUNDARY) {
  // Interior birth group
  if (random.uniform() < ZETA_IB_VERTEX_NOT_TRIANGLE) {
    // Vertex birth move
    sampleInteriorVertexBirth(state, random);
  } else {
    // Triangle birth move
    sampleInteriorTriangleBirth(state, random);
  }
      } else {
  // Boundary birth group
  if (random.uniform() < ZETA_BOUND_NOT_CORNER) {
    // Boundary triangle birth move
    sampleBoundaryTriangleBirth(state, random);
  } else {
    // Corner cut birth move
    sampleCornerCutBirth(state, random);
  }
      }
    } else {
      // Death group
      if (random.uniform() < ZETA_INT_NOT_BOUNDARY) {
  // Interior death group
  sampleInteriorDeath(state, random, reversible);
      } else {
  // Boundary death group
  if (random.uniform() < ZETA_BOUND_NOT_CORNER) {
    // Boundary edge death move
    sampleBoundaryTriangleDeath(state, random, reversible);
  } else {
    // Corner cut death move
    sampleCornerCutDeath(state, random);
  }
      }
    }
  } else {
    // Move/recolor group
    if (random.uniform() < ZETA_MOVE_NOT_RECOLOR) {
      // Move group
      if (random.uniform() < ZETA_INT_NOT_BOUNDARY) {
  // Relocate interior vertex move
  sampleMoveInteriorVertex(state, random);
      } else {
  // Relocate boundary vertex move
  if (random.uniform() < ZETA_BOUND_NOT_CORNER) {
    // Relocate vertex along boundary move
    sampleMoveVertexAlongBd(state, random);
  } else {
    // Relocate boundary vertex past corner move
    sampleMoveBdVertexPastCorner(state, random);
  }
      }
    } else {
      // Recolor move
      sampleRecolor(state, random);
    }
  }
  assert(curMoveType != UNSPECIFIED_MOVE);
  // Update the proposal statistics.
  numProposed[curMoveType]++;
  // std::cerr << "Proposed " << moveNames[int(curMoveType)] << std::endl;
}


void CN94Proposal::sampleInteriorTriangleBirth(const Coloring& state,
                 Arak::Util::Random& random) {
  // Sample the vertices of the triangle.
  Point u = state.randomPoint(random);
  Point v = sampleInCastingBox(u, random);
  Point w = sampleInCastingBox(u, random);
  if (state.interior(v) &&
      state.interior(w) &&
      state.newInteriorTriangle(u, v, w, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = INT_TRIANGLE_BIRTH;
    itb.reset(u, v, w);
  } else {
    // This sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

void CN94Proposal::proposeBdTriangleDeath(const Coloring& state,
            Coloring::VertexHandle prevVertex,
            Coloring::VertexHandle vertex,
            Coloring::VertexHandle nextVertex,
            bool reversible) {
  if (reversible) {
    // Check to see if w satisfies the distance contraint.
    Segment edge(prevVertex->point(), nextVertex->point());
    Point p = edge.supporting_line().projection(vertex->point());
    if (!(CGAL::collinear_are_ordered_along_line(edge.source(),
             p,
             edge.target())
    &&
    (CGAL::squared_distance(vertex->point(), p) <= 
     edge.squared_length()))) {
      // This boundary triangle death could not be reversed by a
      // boundary triangle birth.
      curMoveType = INVALID_MOVE;
      return;
    }
  }
  if (state.deleteBdTriangle(prevVertex, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = BD_TRIANGLE_DEATH;
    btd.reset(prevVertex);
  } else {
    // The sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

void CN94Proposal::proposeIntVertexDeath(const Coloring& state,
           Coloring::VertexHandle prevVertex,
           Coloring::VertexHandle vertex,
           Coloring::VertexHandle nextVertex,
           bool reversible) {
  if (reversible) {
    // Check that this move is reversible.
    Geometry::Segment edge(prevVertex->point(), nextVertex->point());
    Geometry::Point projection = 
      edge.supporting_line().projection(vertex->point());
    // If the vertex is farther than one half of the edge's length
    // from the edge that will be created, then it could not have
    // been created by a birth.  Also, if the vertex's projection
    // does not lie on the edge that will be created, then it could
    // not have been created by a birth.
    double sd = 
      CGAL::to_double(CGAL::squared_distance(vertex->point(), 
               projection));
    double squared_limit = CGAL::to_double(edge.squared_length()) / 4.0;
    if (!((sd <= squared_limit) &&
    CGAL::collinear_are_ordered_along_line(edge.source(),
             projection,
             edge.target()))) {
      // This move could not be reversed by a interior vertex birth.
      curMoveType = INVALID_MOVE;
      return;
    }
  }
  if (state.deleteVertex(vertex, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = INT_VERTEX_DEATH;
    ivd.reset(vertex);
  } else {
    // The sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

void CN94Proposal::proposeIntTriangleDeath(const Coloring& state,
             Coloring::VertexHandle prevVertex,
             Coloring::VertexHandle vertex,
             Coloring::VertexHandle nextVertex,
             bool reversible) {
  // Check if the move is possible.
  if (!state.deleteIntTriangle(prevVertex, Coloring::TEST_IF_VALID_TAG())) {
    // The move is not valid.
    curMoveType = INVALID_MOVE;
    return;
  }
  if (!reversible) {
    curMoveType = INT_TRIANGLE_DEATH;
    itd.reset(vertex, false);
    return;
  }
  // Check if this move could be reversed.
  const Geometry::Point& u = vertex->point();
  const Geometry::Point& v = nextVertex->point();
  const Geometry::Point& w = prevVertex->point();
  int uvClose = inCastingBox(u, v) ? 1 : 0;
  int uwClose = inCastingBox(u, w) ? 1 : 0;
  int vwClose = inCastingBox(v, w) ? 1 : 0;
  int numClose = uvClose + uwClose + vwClose;
  if (numClose < 2) {
    // At least two of the edges are too long.  This triangle could
    // not have been generated by an interior triangle birth.
    curMoveType = INVALID_MOVE;
  } else if (numClose == 2) {
    // This is a "long" triangle that could have been generated by
    // an interior triangle birth.
    curMoveType = INT_TRIANGLE_DEATH;
    itd.reset(vertex, false);
  } else /* numClose == 3 */ {
    // This is a "short" triangle that could have been generated by
    // an interior triangle birth.
    curMoveType = INT_TRIANGLE_DEATH;
    itd.reset(vertex, true);
  }
}

void CN94Proposal::sampleInteriorDeath(const Coloring& state,
               Arak::Util::Random& random,
               bool reversible) {
  // If there are no vertices, the sampling failed.
  if (state.numVertices(Coloring::Vertex::INTERIOR) == 0) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Sample a vertex uniformly from the interior of the coloring.
  Coloring::VertexHandle vertex = 
    state.randomVertex(Coloring::Vertex::INTERIOR, random);
  // Examine the neighborhood of the vertex to determine the move
  // type.
  Coloring::VertexHandle nextVertex = vertex->getNextVertex();
  Coloring::VertexHandle prevVertex = vertex->getPrevVertex();
  if ((nextVertex->type() == Coloring::Vertex::BOUNDARY) &&
      (prevVertex->type() == Coloring::Vertex::BOUNDARY)) {
    // Both adjacent vertices are on the boundary.
    if ((nextVertex->getNextBdEdge()->getNextVertex() == prevVertex) ||
  (nextVertex->getPrevBdEdge()->getPrevVertex() == prevVertex)) {
      // The adjacent vertices are connected by a single boundary
      // edge, so this is a structure that could have been built by a
      // boundary triangle birth.  Propose a boundary triangle death
      // move.
      proposeBdTriangleDeath(state, prevVertex, vertex, nextVertex,
           reversible);
    } else if (nextVertex->getNextCornerVertex() != 
         prevVertex->getNextCornerVertex()) {
      // The adjacent vertices are on distinct boundary faces.  Try to
      // propose an interior vertex death move.  We must check that
      // this vertex could have been created by an interior vertex
      // birth move.
      proposeIntVertexDeath(state, prevVertex, vertex, nextVertex,
          reversible);
    } else {
      // The adjacent vertices are on the same boundary face but are
      // not connected by a single boundary edge.  In this case, the
      // sampling failed.
      curMoveType = INVALID_MOVE;
    }
  } else if ((nextVertex->type() == Coloring::Vertex::INTERIOR) &&
       (prevVertex->type() == Coloring::Vertex::INTERIOR) &&
       (nextVertex->getNextVertex() == prevVertex)) {
    // This vertex is part of an interior triangle.  Propose an
    // interior triangle death move.
    proposeIntTriangleDeath(state, prevVertex, vertex, nextVertex,
          reversible);
  } else {
    // Otherwise, propose an interior vertex death move.
    proposeIntVertexDeath(state, prevVertex, vertex, nextVertex,
        reversible);
  }
}

void CN94Proposal::sampleRecolor(const Coloring& state,
         Arak::Util::Random& random) {
  // If there are no edges, sample an invalid move.
  if (state.numInteriorEdges() < 2) {
    // This sampling failed.
    curMoveType = INVALID_MOVE;
    return;
  }
  Coloring::IntEdgeHandle edge1 = 
    state.getInteriorEdge(random.uniform(state.numInteriorEdges()));
  Coloring::IntEdgeHandle edge2 = 
    state.getInteriorEdge(random.uniform(state.numInteriorEdges()));
  while (edge1 == edge2)
    edge2 = state.getInteriorEdge(random.uniform(state.numInteriorEdges()));
  bool joinSourceToSourcePossible = 
    !CGAL::do_intersect(Segment(*(edge1->getPrevVertex()), 
        *(edge2->getPrevVertex())),
      Segment(*(edge1->getNextVertex()), 
        *(edge2->getNextVertex())));
  bool joinSourceToTargetPossible = 
    !CGAL::do_intersect(Segment(*(edge1->getPrevVertex()), 
        *(edge2->getNextVertex())),
      Segment(*(edge1->getNextVertex()), 
        *(edge2->getPrevVertex())));
  assert(joinSourceToSourcePossible || joinSourceToTargetPossible);
  bool convex = joinSourceToSourcePossible xor joinSourceToTargetPossible;
  bool joinSourceToSource = 
    (joinSourceToSourcePossible && !joinSourceToTargetPossible) ||
    (joinSourceToSourcePossible && joinSourceToTargetPossible &&
     random.bernoulli());
  if (state.recolorQuadrilateral(edge1, edge2, joinSourceToSource,
         Coloring::TEST_IF_VALID_TAG())) {
    // Propose the recoloring move.
    curMoveType = RECOLOR;
    rec.reset(edge1, edge2, joinSourceToSource, convex);
  } else {
    curMoveType = INVALID_MOVE;
  }
}

void CN94Proposal::sampleMoveBdVertexPastCorner(const Coloring& state, 
            Arak::Util::Random& random) {
  // Sample a corner vertex.
  Coloring::VertexHandle corner = state.randomVertex(Coloring::Vertex::CORNER);
  // Sample an incident edge.
  bool useNext = random.bernoulli();
  Coloring::BdEdgeHandle edge = (useNext ? 
         corner->getNextBdEdge() : 
         corner->getPrevBdEdge());
  Coloring::VertexHandle bdVertex = (useNext ? 
             edge->getNextVertex() : 
             edge->getPrevVertex());
  // If the next vertex is not a boundary vertex (i.e., it is a
  // corner), then this sampling failed.
  if (bdVertex->type() != Coloring::Vertex::BOUNDARY) {
    curMoveType = INVALID_MOVE;
    return;
  }
  Coloring::VertexHandle adjVertex = (bdVertex->hasNextIntEdge() ? 
              bdVertex->getNextVertex() : 
              bdVertex->getPrevVertex());
  // If the adjacent vertex is a boundary vertex on the target face,
  // this sampling failed.
  if ((adjVertex->type() == Coloring::Vertex::BOUNDARY) &&
      ((useNext && (adjVertex->getNextCornerVertex() == corner)) ||
       (!useNext && (adjVertex->getPrevCornerVertex() == corner)))) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Sample a new location for the vertex on the other edge incident
  // to the corner.
  Coloring::BdEdgeHandle otherEdge = (useNext ? 
              corner->getPrevBdEdge() :
              corner->getNextBdEdge()); 
  Point p = otherEdge->randomPoint(random);
  // Check that this move is valid.
  if (state.moveBdVertexPastCorner(edge, p, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = MOVE_BD_VERTEX_PAST_CORNER;
    mbc.reset(edge, p);
  } else {
    // This sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

void CN94Proposal::sampleInteriorVertexBirth(const Coloring& state,
               Arak::Util::Random& random) {
  // If there are no edges in the interior, sample an invalid move.
  if (state.numInteriorEdges() == 0) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Sample an edge uniformly from the interior.
  const Coloring::IntEdgeHandle e = 
    state.getInteriorEdge(random.uniform(state.numInteriorEdges()));
  // Sample a point uniformly on the edge.
  Point q = e->randomPoint(random);
  // Compute the normal vector to the edge.
  Vector v = e->to_vector().perpendicular(CGAL::CLOCKWISE);
  // Normalize the vector.
  v = v / sqrt(CGAL::to_double(v * v));
  // Choose the distance from e.
  double halfLength = e->length() / 2.0;
  Point p = q + (v * random.uniform(-halfLength, halfLength));
  // Check that the move is valid.
  if (state.interior(p) &&
      state.splitEdge(e, p, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = INT_VERTEX_BIRTH;
    ivb.reset(e, p);
  } else {
    // This sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

void CN94Proposal::sampleMoveInteriorVertex(const Coloring& state, 
              Arak::Util::Random& random) {
  // If there are no vertices in the interior, sample an invalid move.
  if (state.numVertices(Coloring::Vertex::INTERIOR) == 0) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Sample a vertex uniformly from the interior.
  Coloring::VertexHandle v = 
    state.randomVertex(Coloring::Vertex::INTERIOR, random);
  Point p = sampleInCastingBox(v->point(), random);
  // Check to see if the move is valid.
  if (state.interior(p) &&
      state.moveIntVertex(v, p, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = MOVE_INT_VERTEX;
    miv.reset(v, p);
  } else {
    // This sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

void CN94Proposal::sampleMoveVertexAlongBd(const Coloring& state, 
             Arak::Util::Random& random) {
  // If there are no vertices on the boundary, sample an invalid move.
  if (state.numVertices(Coloring::Vertex::BOUNDARY) == 0) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Sample a boundary vertex uniformly.
  Coloring::VertexHandle v = 
    state.randomVertex(Coloring::Vertex::BOUNDARY, random);
  Point a = *(v->getPrevBdEdge()->getPrevVertex());
  Point b = *(v->getNextBdEdge()->getNextVertex());
  Geometry::Vector va(CGAL::ORIGIN, a);
  Geometry::Vector vb(CGAL::ORIGIN, b);
  Geometry::Kernel::FT x(random.uniform());
  Geometry::Vector vp = (va * x) + vb - (vb * x);
  Geometry::Point p(vp.x(), vp.y());
#ifdef EXACT_GEOM_CXN
  assert(v->getPrevBdEdge()->has_on(p) || v->getNextBdEdge()->has_on(p));
#endif
  // Check that the move is valid.
  if ((a != p) &&
      (b != p) &&
      state.moveVertexAlongBd(v, p, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = MOVE_BD_VERTEX_ALONG_BD;
    mbb.reset(v, p);
  } else {
    // The sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

void CN94Proposal::sampleBoundaryTriangleBirth(const Coloring& state,
                 Arak::Util::Random& random) {
  // Sample the boundary edge.
  const Coloring::BdEdgeHandle e = 
    state.getBoundaryEdge(random.uniform(state.numBoundaryEdges()));
  // Sample the boundary vertices of the triangle.
  Point u = e->randomPoint(random);
  Point v = e->randomPoint(random);
  // Reject the move if the sampled points are vertices of the
  // boundary edge; this can happen when the boundary edge is short.
  if ((e->source() == u) ||
      (e->target() == u) ||
      (e->source() == v) ||
      (e->target() == v)) {
    curMoveType = INVALID_MOVE;
    return;
  }
  /* 
   * Sample the interior vertex.
   */
  // Sample a point uniformly on the new boundary edge.
  Point q = random.uniform(u, v);
  // Compute the normal vector to the edge.
  Vector n = e->to_vector().perpendicular(CGAL::CLOCKWISE);
  // Normalize the vector.
  n = n / sqrt(CGAL::to_double(n * n));
  // Choose the distance from e.
  double distance = 
    random.uniform(0.0, sqrt(CGAL::to_double(CGAL::squared_distance(u, v))));
  Point w = q + (n * distance);
  // If we chose the wrong perpendicular direction, flip it.
  if (!state.interior(w))
    w = q - n * distance;
  // Check that the move is valid.
  if ((u != v) && (v != w) && (w != u) &&
      state.interior(w) &&
      state.newBoundaryTriangle(e, u, v, w, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = BD_TRIANGLE_BIRTH;
    btb.reset(e, u, v, w);
  } else {
    // The sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

void CN94Proposal::sampleBoundaryTriangleDeath(const Coloring& state,
                 Arak::Util::Random& random,
                 bool reversible) {
  // If there are no edges on the boundary, sample an invalid move.
  if (state.numBoundaryEdges() == 0) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Sample a boundary edge.
  const Coloring::BdEdgeHandle e = 
    state.getBoundaryEdge(random.uniform(state.numBoundaryEdges()));
  // Check to see if this edge is part of a boundary triangle.
  const Coloring::VertexHandle u = e->getPrevVertex();
  const Coloring::VertexHandle v = e->getNextVertex();
  // Make sure both vertices are boundary vertices.
  if ((u->type() == Coloring::Vertex::CORNER) || 
      (v->type() == Coloring::Vertex::CORNER)) {
    curMoveType = INVALID_MOVE; return;
  }
  // Make sure they are connected to a common central vertex.
  if ((u->hasNextIntEdge() == v->hasNextIntEdge()) ||
      (u->hasNextIntEdge() && 
       (u->getNextVertex() != v->getPrevVertex())) ||
      (v->hasNextIntEdge() && 
       (v->getNextVertex() != u->getPrevVertex()))) {
    curMoveType = INVALID_MOVE; return;
  }
  // Propose the boundary triangle death move.
  if (u->hasNextIntEdge())
    proposeBdTriangleDeath(state, u, u->getNextVertex(), v, reversible);
  else
    proposeBdTriangleDeath(state, v, v->getNextVertex(), u, reversible);
}

void CN94Proposal::sampleCornerCutBirth(const Coloring& state,
          Arak::Util::Random& random) {
  Coloring::VertexHandle corner = 
    state.randomVertex(Coloring::Vertex::CORNER, random);
  Coloring::BdEdgeHandle prevBdEdge = corner->getPrevBdEdge();
  Coloring::BdEdgeHandle nextBdEdge = corner->getNextBdEdge();
  Point u = prevBdEdge->randomPoint(random);
  Point v = nextBdEdge->randomPoint(random);
  // Check that (u, v) is a compatible edge.  Also ensure that the
  // sampled points are not vertices; this can happen when the
  // boundary edges are short.
  if ((u == prevBdEdge->source()) ||
      (u == prevBdEdge->target()) ||
      (v == nextBdEdge->source()) ||
      (v == nextBdEdge->target()) ||
      !state.newCornerTriangle(corner, u, v, Coloring::TEST_IF_VALID_TAG())) {
    // This sampling failed.
    curMoveType = INVALID_MOVE;
  } else {
    curMoveType = CORNER_CUT_BIRTH;
    ccb.reset(corner, u, v);
  }
}

void CN94Proposal::sampleCornerCutDeath(const Coloring& state,
          Arak::Util::Random& random) {
  Coloring::VertexHandle corner = 
    state.randomVertex(Coloring::Vertex::CORNER, random);
  Coloring::VertexHandle prevBdVertex = 
    corner->getPrevBdEdge()->getPrevVertex();
  Coloring::VertexHandle nextBdVertex = 
    corner->getNextBdEdge()->getNextVertex();
  if ((prevBdVertex->type() == Coloring::Vertex::BOUNDARY) &&
      (nextBdVertex->type() == Coloring::Vertex::BOUNDARY) &&
      ((prevBdVertex->hasNextIntEdge() && 
  (prevBdVertex->getNextVertex() == nextBdVertex) &&
  (state.deleteCornerTriangle(prevBdVertex, 
            Coloring::TEST_IF_VALID_TAG())))
       ||
       (nextBdVertex->hasNextIntEdge() && 
  (nextBdVertex->getNextVertex() == prevBdVertex) &&
  (state.deleteCornerTriangle(nextBdVertex, 
            Coloring::TEST_IF_VALID_TAG()))))) {
    ccd.reset(corner);
    curMoveType = CORNER_CUT_DEATH;
  } else {
    // This sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

double CN94Proposal::ll(const Coloring& c) {
  switch (curMoveType) {
  case INT_TRIANGLE_BIRTH: {
    double a = castingArea();
    double A = c.area();
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        ZETA_INT_NOT_BOUNDARY *
        (1.0 - ZETA_IB_VERTEX_NOT_TRIANGLE) * 
        (inCastingBox(itb.v, itb.w) ? 6.0 : 2.0) / (A * a * a));
  }
  case INT_TRIANGLE_DEATH: 
    // The extra factor of three is necessary because there are three
    // ways to propose removal of an interior triangle (one per
    // vertex).
    return ln(ZETA_BD_NOT_MR * 
        (1.0 - ZETA_BIRTH_NOT_DEATH) *
        ZETA_INT_NOT_BOUNDARY *
        3.0 / 
        double(c.numVertices(Coloring::Vertex::INTERIOR)));
  case RECOLOR: {
    int n = c.numInteriorEdges();
    return ln(((1.0 - ZETA_BD_NOT_MR) * 
         (1.0 - ZETA_MOVE_NOT_RECOLOR) *
         2.0 * (rec.isConvex() ? 1.0 : 0.5)) /
        double(n * (n - 1)));
  }
  case MOVE_INT_VERTEX: {
    return ln((1.0 - ZETA_BD_NOT_MR) * 
        ZETA_MOVE_NOT_RECOLOR *
        ZETA_INT_NOT_BOUNDARY /
        (double(c.numVertices(Coloring::Vertex::INTERIOR)) *
         castingArea()));
  }
  case MOVE_BD_VERTEX_ALONG_BD: {
    return ln((1.0 - ZETA_BD_NOT_MR) * 
        ZETA_MOVE_NOT_RECOLOR *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        ZETA_BOUND_NOT_CORNER /
        (double(c.numVertices(Coloring::Vertex::BOUNDARY)) *
         (mbb.vertex->getNextBdEdge()->length() + 
    mbb.vertex->getPrevBdEdge()->length())));
  }
  case MOVE_BD_VERTEX_PAST_CORNER: {
    // Find the other edge incident to the corner.
    Coloring::BdEdgeHandle ep = 
      (mbc.edge->getPrevVertex()->type() == Coloring::Vertex::CORNER) ?
      mbc.edge->getPrevVertex()->getPrevBdEdge() :
      mbc.edge->getNextVertex()->getNextBdEdge();
    return ln((1.0 - ZETA_BD_NOT_MR) * 
        ZETA_MOVE_NOT_RECOLOR *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        (1.0 - ZETA_BOUND_NOT_CORNER) /
        (2.0 * double(c.numVertices(Coloring::Vertex::CORNER)) 
         * ep->length()));
  }
  case INT_VERTEX_BIRTH: {
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        ZETA_INT_NOT_BOUNDARY *
        ZETA_IB_VERTEX_NOT_TRIANGLE /
        (double(c.numInteriorEdges()) 
         * CGAL::to_double(squared_distance(*(ivb.prevVertex), 
              *(ivb.nextVertex)))));
  }
  case INT_VERTEX_DEATH: {
    return ln(ZETA_BD_NOT_MR * 
        (1.0 - ZETA_BIRTH_NOT_DEATH) *
        ZETA_INT_NOT_BOUNDARY /
        double(c.numVertices(Coloring::Vertex::INTERIOR)));
  }
  case BD_TRIANGLE_BIRTH: {
    double a = CGAL::to_double(btb.edge->segment().squared_length());
    double b = CGAL::to_double(squared_distance(btb.u, btb.v));
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        ZETA_BOUND_NOT_CORNER *
        2.0 / 
        (double(c.numBoundaryEdges()) * a * b));
  }
  case BD_TRIANGLE_DEATH:  {
    return ln(// Probability proposed by interior vertex death
        (ZETA_BD_NOT_MR * 
         (1.0 - ZETA_BIRTH_NOT_DEATH) *
         ZETA_INT_NOT_BOUNDARY /
         double(c.numVertices(Coloring::Vertex::INTERIOR)))
        +
        // Probability proposed by boundary edge death
        (ZETA_BD_NOT_MR * 
         (1.0 - ZETA_BIRTH_NOT_DEATH) *
         (1.0 - ZETA_INT_NOT_BOUNDARY) *
         ZETA_BOUND_NOT_CORNER / 
         double(c.numBoundaryEdges())));
  }
  case CORNER_CUT_BIRTH: {
    double lp = ccb.corner->getPrevBdEdge()->length();
    double np = ccb.corner->getNextBdEdge()->length();
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        (1.0 - ZETA_BOUND_NOT_CORNER) /
        (double(c.numVertices(Coloring::Vertex::CORNER)) * lp * np));    
  }
  case CORNER_CUT_DEATH:
    return ln(ZETA_BD_NOT_MR * 
        (1.0 - ZETA_BIRTH_NOT_DEATH) *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        (1.0 - ZETA_BOUND_NOT_CORNER) /
        double(c.numVertices(Coloring::Vertex::CORNER)));
  case INVALID_MOVE:
    // We want to force the sampler to reject the current move.  To do
    // this, we say its proposal probability is one and its reverse
    // proposal probability is zero.  TODO: this should be done more
    // cleanly.
    return 0.0;
  default:
    assert(false);
  }
}
  
double CN94Proposal::rll(const Coloring& c) {
  switch (curMoveType) {
  case INT_TRIANGLE_BIRTH:
    // Reverse of: INT_TRIANGLE_DEATH
    return ln(ZETA_BD_NOT_MR * 
        (1.0 - ZETA_BIRTH_NOT_DEATH) *
        ZETA_INT_NOT_BOUNDARY *
        3.0 / 
        double(3 + c.numVertices(Coloring::Vertex::INTERIOR)));
  case INT_TRIANGLE_DEATH: {
    // Reverse of: INT_TRIANGLE_BIRTH
    double a = castingArea();
    double A = c.area();
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        ZETA_INT_NOT_BOUNDARY *
        (1.0 - ZETA_IB_VERTEX_NOT_TRIANGLE) * 
        (itd.isShort ? 6.0 : 2.0) / 
        (A * a * a));
  } 
  case RECOLOR: {
    // Reverse of: RECOLOR
    int n = c.numInteriorEdges();
    return ln((1.0 - ZETA_BD_NOT_MR) * 
        (1.0 - ZETA_MOVE_NOT_RECOLOR) *
        2.0 * (rec.isConvex() ? 1.0 : 0.5) /
        double(n * (n - 1)));
  }
  case MOVE_INT_VERTEX: {
    // Reverse of: MOVE_INT_VERTEX
    return ln((1.0 - ZETA_BD_NOT_MR) * 
        ZETA_MOVE_NOT_RECOLOR *
        ZETA_INT_NOT_BOUNDARY /
        (double(c.numVertices(Coloring::Vertex::INTERIOR)) *
         castingArea()));
  }
  case MOVE_BD_VERTEX_ALONG_BD: {
    // Reverse of: MOVE_BD_VERTEX_ALONG_BD
    return ln((1.0 - ZETA_BD_NOT_MR) * 
        ZETA_MOVE_NOT_RECOLOR *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        ZETA_BOUND_NOT_CORNER /
        (double(c.numVertices(Coloring::Vertex::BOUNDARY)) *
         (mbb.vertex->getNextBdEdge()->length() + 
    mbb.vertex->getPrevBdEdge()->length())));
  }
  case MOVE_BD_VERTEX_PAST_CORNER: {
    // Reverse of: MOVE_BD_VERTEX_PAST_CORNER
    Coloring::BdEdgeHandle e = mbc.edge;
    Coloring::BdEdgeHandle f = 
      ((e->getPrevVertex()->type() == Coloring::Vertex::CORNER) ?
       e->getNextVertex()->getNextBdEdge() :
       e->getPrevVertex()->getPrevBdEdge());
    return ln((1.0 - ZETA_BD_NOT_MR) * 
        ZETA_MOVE_NOT_RECOLOR *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        (1.0 - ZETA_BOUND_NOT_CORNER) /
        (2.0 * double(c.numVertices(Coloring::Vertex::CORNER)) 
         * (e->length() + f->length())));
  }
  case INT_VERTEX_BIRTH: {
    // Reverse of: INT_VERTEX_DEATH
    return ln(ZETA_BD_NOT_MR * 
        (1.0 - ZETA_BIRTH_NOT_DEATH) *
        ZETA_INT_NOT_BOUNDARY /
        double(1 + c.numVertices(Coloring::Vertex::INTERIOR)));
  }
  case INT_VERTEX_DEATH: {
    // Reverse of: INT_VERTEX_BIRTH
    double l = ivd.newEdgeLength;
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        ZETA_INT_NOT_BOUNDARY *
        ZETA_IB_VERTEX_NOT_TRIANGLE /
        (double(c.numInteriorEdges() - 1) * l * l));
  }
  case BD_TRIANGLE_BIRTH: {
    // Reverse of: BD_TRIANGLE_DEATH
    return ln(// Probability reverse proposed by interior vertex death
        (ZETA_BD_NOT_MR * 
         (1.0 - ZETA_BIRTH_NOT_DEATH) *
         ZETA_INT_NOT_BOUNDARY /
         double(c.numVertices(Coloring::Vertex::INTERIOR) + 1))
        +
        // Probability reverse proposed by boundary edge death
        (ZETA_BD_NOT_MR * 
         (1.0 - ZETA_BIRTH_NOT_DEATH) *
         (1.0 - ZETA_INT_NOT_BOUNDARY) *
         ZETA_BOUND_NOT_CORNER / 
         double(c.numBoundaryEdges() + 2)));
  }
  case BD_TRIANGLE_DEATH: {
    // Reverse of: BD_TRIANGLE_BIRTH
    double a = CGAL::to_double(squared_distance(*(btd.prevBoundaryVertex),
            *(btd.nextBoundaryVertex)));
    double b = CGAL::to_double(squared_distance(btd.up, btd.vp));
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        ZETA_BOUND_NOT_CORNER *
        2.0 / 
        (double(c.numBoundaryEdges() - 2) * a * b));
  }
  case CORNER_CUT_BIRTH: {
    // Reverse of: CORNER_CUT_DEATH
    return ln(ZETA_BD_NOT_MR * 
        (1.0 - ZETA_BIRTH_NOT_DEATH) *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        (1.0 - ZETA_BOUND_NOT_CORNER) /
        double(c.numVertices(Coloring::Vertex::CORNER)));
  }
  case CORNER_CUT_DEATH: {
    // Reverse of: CORNER_CUT_BIRTH
    Coloring::BdEdgeHandle prevBdEdge = ccb.corner->getPrevBdEdge();
    Coloring::BdEdgeHandle nextBdEdge = ccb.corner->getNextBdEdge();
    double lp = prevBdEdge->length() + 
      prevBdEdge->getPrevVertex()->getPrevBdEdge()->length();
    double np = nextBdEdge->length() + 
      nextBdEdge->getNextVertex()->getNextBdEdge()->length();
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        (1.0 - ZETA_BOUND_NOT_CORNER) /
        (double(c.numVertices(Coloring::Vertex::CORNER)) * lp * np));
  }
  case INVALID_MOVE:
    // We want to force the sampler to reject the current move.  To do
    // this, we say its proposal probability is one and its reverse
    // proposal probability is zero.  TODO: this should be done more
    // cleanly.
    return -std::numeric_limits<double>::infinity();
  default:
    assert(false);
  }
}

ModifiedCN94Proposal::ModifiedCN94Proposal(const Arak::Util::PropertyMap& props)
  : CN94Proposal(props) { 
  // Initialize the parameters of the proposal.
  using namespace Arak::Util;
  assert(parse(getp(props, 
        "arak.modified_cn94.zeta.local_not_global_recolor"), 
         ZETA_LOCAL_NOT_GLOBAL_RECOLOR));
  assert(parse(getp(props, 
        "arak.modified_cn94.zeta.local_not_global_bd_tri_birth"), 
         ZETA_LOCAL_NOT_GLOBAL_BD_TRI_BIRTH));
  assert(parse(getp(props, 
        "arak.modified_cn94.zeta.slide_not_move_int_vertex"), 
         ZETA_SLIDE_NOT_MOVE_INT_VERTEX));
}

void ModifiedCN94Proposal::sampleRecolor(const Coloring& state,
           Arak::Util::Random& random) {
  // With some probability use the global recolor move.
  if (random.uniform() >= ZETA_LOCAL_NOT_GLOBAL_RECOLOR) {
    CN94Proposal::sampleRecolor(state, random);
    return;
  }
  // If there are no edges, sample an invalid move.
  if (state.numInteriorEdges() < 2) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Sample an interior edge uniformly from the coloring.
  Coloring::IntEdgeHandle edge1 = 
    state.getInteriorEdge(random.uniform(state.numInteriorEdges()));
  // Sample a cell from the interior edge index that is crossed by
  // this edge.
  const Coloring::InteriorEdge::CellEntryList& cellEntryList = 
    edge1->cellEntries();
  int numEntries = cellEntryList.size();
  Coloring::InteriorEdge::CellEntryList::const_iterator eit = 
    cellEntryList.begin();
  std::advance(eit, random.uniform(numEntries));
  const Coloring::InteriorEdge::CellEntry entry = *eit;
  const Coloring::InteriorEdgeIndex::Cell& cell = entry.getCell();
  // Get the list of edges that cross this cell.
  const Coloring::InteriorEdgeIndex::Cell::ItemList& intEdgeList = 
    cell.getItemList();
  int numEdges = intEdgeList.size();
  // If there is only one edge crossing this cell, the move is invalid.
  if (numEdges == 1) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Sample an edge from the list that is distinct from edge1.
  int i = random.uniform(numEdges - 1);
  typedef Coloring::InteriorEdgeIndex::Cell::ItemList::const_iterator Iterator;
  Coloring::IntEdgeHandle edge2;
  for (Iterator it = intEdgeList.begin(); it != intEdgeList.end(); ++it) {
    edge2 = *it;
    if (edge2 == edge1) {
      ++it; 
      continue;
    }
    if (i-- == 0) break;
  }
  // Sample the way in which we will recolor.
  bool joinSourceToSourcePossible = 
    !CGAL::do_intersect(Segment(*(edge1->getPrevVertex()), 
        *(edge2->getPrevVertex())),
      Segment(*(edge1->getNextVertex()), 
        *(edge2->getNextVertex())));
  bool joinSourceToTargetPossible = 
    !CGAL::do_intersect(Segment(*(edge1->getPrevVertex()), 
        *(edge2->getNextVertex())),
      Segment(*(edge1->getNextVertex()), 
        *(edge2->getPrevVertex())));
  assert(joinSourceToSourcePossible || joinSourceToTargetPossible);
  bool convex = joinSourceToSourcePossible xor joinSourceToTargetPossible;
  bool joinSourceToSource = 
    (joinSourceToSourcePossible && !joinSourceToTargetPossible) ||
    (joinSourceToSourcePossible && joinSourceToTargetPossible &&
     random.bernoulli());
  // Propose the recoloring move.
  if (state.recolorQuadrilateral(edge1, edge2, joinSourceToSource,
         Coloring::TEST_IF_VALID_TAG())) {
    // Propose the recoloring move.
    curMoveType = RECOLOR;
    rec.reset(edge1, edge2, joinSourceToSource, convex);
  } else {
    curMoveType = INVALID_MOVE;
  }
}

void ModifiedCN94Proposal::sampleBoundaryTriangleBirth(const Coloring& state,
                   Arak::Util::Random& random) {
  // With some probability use the global recolor move.
  if (random.uniform() >= ZETA_LOCAL_NOT_GLOBAL_BD_TRI_BIRTH) {
    CN94Proposal::sampleBoundaryTriangleBirth(state, random);
    return;
  }
  // Sample the boundary edge.
  const Coloring::BdEdgeHandle e = 
    state.getBoundaryEdge(random.uniform(state.numBoundaryEdges()));
  // Sample the first boundary vertex of the triangle.
  Point u = e->randomPoint(random);
  // Compute a vector aligned with the edge of length castingBoxRadius.
  Vector a = e->to_vector();
  a = a * (castingBoxRadius / sqrt(CGAL::to_double(a * a)));
  // Sample the second vertex of the triangle.
  Point v = random.uniform(u + a, u - a);
  // If this point is not on the edge, the move is rejected.
  if (!CGAL::collinear_are_ordered_along_line(e->source(), v, e->target())) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Reject the move if the sampled points are vertices of the
  // boundary edge; this can happen when the boundary edge is short.
  if ((e->source() == u) ||
      (e->target() == u) ||
      (e->source() == v) ||
      (e->target() == v)) {
    curMoveType = INVALID_MOVE;
    return;
  }
  /* 
   * Sample the interior vertex.
   */
  // Sample a point uniformly on the new boundary edge.
  Point q = random.uniform(u, v);
  // Compute the normal vector to the edge.
  Vector n = e->to_vector().perpendicular(CGAL::CLOCKWISE);
  // Normalize the vector.
  n = n / sqrt(CGAL::to_double(n * n));
  // Choose the distance from e.
  double distance = 
    random.uniform(0.0, sqrt(CGAL::to_double(CGAL::squared_distance(u, v))));
  Point w = q + (n * distance);
  // If we chose the wrong perpendicular direction, flip it.
  if (!state.interior(w))
    w = q - n * distance;
  // Check that the move is valid.
  if ((u != v) && (v != w) && (w != u) &&
      state.interior(w) &&
      state.newBoundaryTriangle(e, u, v, w, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = BD_TRIANGLE_BIRTH;
    btb.reset(e, u, v, w);
  } else {
    // The sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

inline double localRecolorSampleProb(const Coloring& state,
             Coloring::IntEdgeHandle edge1,
             Coloring::IntEdgeHandle edge2) {
  int n = state.numInteriorEdges();
  int k1 = edge1->cellEntries().size();
  int k2 = edge2->cellEntries().size();
  double p = 0.0;
  // We must iterate over all grid cells that are crossed by edge1 and
  // edge2.  To do this, we iterate over edge1's cells.
  typedef Coloring::InteriorEdgeIndex::Cell Cell;
  typedef Coloring::InteriorEdge::CellEntry CellEntry;
  typedef Coloring::InteriorEdge::CellEntryList CellEntryList;
  typedef CellEntryList::const_iterator CellEntryIterator;
  const CellEntryList& cel1 = edge1->cellEntries();
  for (CellEntryIterator it1 = cel1.begin(); it1 != cel1.end(); it1++) {
    const Cell& c1 = it1->getCell();
    const Cell::ItemList& edgesInCell = c1.getItemList();
    if (std::find(edgesInCell.begin(), 
      edgesInCell.end(), edge2) != edgesInCell.end()) {
      // This cell is crossed by edge1 and edge2.
      int size = edgesInCell.size();
      // Add in probability we first sampled edge1, then sampled this
      // cell from the list of cells containing edge1, then sampled
      // edge2 from the other edges in this cell.
      p += 1.0 / double(n * k1 * (size - 1));
      // Add in probability we first sampled edge2, then sampled this
      // cell from the list of cells containing edge2, then sampled
      // edge1 from the other edges in this cell.
      p += 1.0 / double(n * k2 * (size - 1));
    }
  }    
  return p;
}

void ModifiedCN94Proposal::sampleSlideInteriorVertex(const Coloring& state,
                 Arak::Util::Random& random) {
  // If there are no vertices in the interior, sample an invalid move.
  if (state.numVertices(Coloring::Vertex::INTERIOR) == 0) {
    curMoveType = INVALID_MOVE;
    return;
  }
  // Sample a vertex uniformly from the interior.
  Coloring::VertexHandle u = 
    state.randomVertex(Coloring::Vertex::INTERIOR, random);
  // Sample an edge incident to u.
  bool next = random.bernoulli();
  Coloring::IntEdgeHandle e = next ? u->getNextIntEdge() : u->getPrevIntEdge();
  Coloring::VertexHandle v = next ? e->getNextVertex() : e->getPrevVertex();
  // Sample the new location.
  Geometry::Vector e_vec = (v->point() - u->point());
  Point p = random.uniform(u->point() + (e_vec / 2.0),
         u->point() - e_vec);
  // Check to see if the move is valid.
  if (state.interior(p) &&
      state.moveIntVertex(u, p, Coloring::TEST_IF_VALID_TAG())) {
    curMoveType = MOVE_INT_VERTEX; // slide_not_move flag indicates slide
    slide_not_move = true;
    siv.reset(u, next, p);
  } else {
    // This sampling failed.
    curMoveType = INVALID_MOVE;
  }
}

void ModifiedCN94Proposal::sampleMoveInteriorVertex(const Coloring& state,
                Arak::Util::Random& random) {
  // With some probability sample a move interior vertex move.
  if (random.uniform() >= ZETA_SLIDE_NOT_MOVE_INT_VERTEX) {
    CN94Proposal::sampleMoveInteriorVertex(state, random);
    slide_not_move = false;
    return;
  }
  // With the remaining probability, sample a slide vertex move.
  sampleSlideInteriorVertex(state, random);
}

double ModifiedCN94Proposal::ll(const Coloring& state) {
  if (curMoveType == RECOLOR) {
    int n = state.numInteriorEdges();
    double p = localRecolorSampleProb(state, rec.edge1, rec.edge2);
    return ln((1.0 - ZETA_BD_NOT_MR) * 
        (1.0 - ZETA_MOVE_NOT_RECOLOR) *
        (// Global probability
         ((1.0 - ZETA_LOCAL_NOT_GLOBAL_RECOLOR) *
    2.0 * (rec.isConvex() ? 1.0 : 0.5) /
    double(n * (n - 1)))
         + 
         // Local probability
         (ZETA_LOCAL_NOT_GLOBAL_RECOLOR * p)));
  } else if (curMoveType == BD_TRIANGLE_BIRTH) {
    int n = state.numBoundaryEdges();
    double a = CGAL::to_double(btb.edge->segment().squared_length());
    double b = CGAL::to_double(squared_distance(btb.u, btb.v));
    bool couldHaveBeenLocal = (b <= castingBoxRadius * castingBoxRadius);
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        ZETA_BOUND_NOT_CORNER *
        (
         // Global probability
         ((1.0 - ZETA_LOCAL_NOT_GLOBAL_BD_TRI_BIRTH) *
    2.0 / 
    (double(n) * a * b))
         +
         // Local probability
         (couldHaveBeenLocal ? 
    (ZETA_LOCAL_NOT_GLOBAL_BD_TRI_BIRTH *
     2.0 / 
     (double(n) * sqrt(a) * 2.0 * castingBoxRadius * b)) 
    :
    0.0)));
  } else if (curMoveType == MOVE_INT_VERTEX) {
    if (slide_not_move) {
      // Slide interior vertex.
      Coloring::IntEdgeHandle e = 
  siv.next ? siv.vertex->getNextIntEdge() : siv.vertex->getPrevIntEdge();
      return ln((1.0 - ZETA_BD_NOT_MR) * 
    ZETA_MOVE_NOT_RECOLOR *
    ZETA_INT_NOT_BOUNDARY *
    ZETA_SLIDE_NOT_MOVE_INT_VERTEX /
    (double(state.numVertices(Coloring::Vertex::INTERIOR)) /
     (3.0 * e->length())));
      // The factor of 3 comes from 2 (edges) * (1.5 * length)
    } else {
      // Move interior vertex.
      return ln((1.0 - ZETA_BD_NOT_MR) * 
    ZETA_MOVE_NOT_RECOLOR *
    ZETA_INT_NOT_BOUNDARY *
    (1.0 - ZETA_SLIDE_NOT_MOVE_INT_VERTEX) /
    (double(state.numVertices(Coloring::Vertex::INTERIOR)) *
     castingArea()));
    }
  } else
    return CN94Proposal::ll(state);
}
  
double ModifiedCN94Proposal::rll(const Coloring& state) {
  if (curMoveType == RECOLOR) {
    // Reverse of: RECOLOR
    int n = state.numInteriorEdges();
    double p = localRecolorSampleProb(state, rec.edge1, rec.edge2);
    return ln((1.0 - ZETA_BD_NOT_MR) * 
        (1.0 - ZETA_MOVE_NOT_RECOLOR) *
        (// Global probability
         ((1.0 - ZETA_LOCAL_NOT_GLOBAL_RECOLOR) *
    2.0 * (rec.isConvex() ? 1.0 : 0.5) /
    double(n * (n - 1)))
         + 
         // Local probability
         (ZETA_LOCAL_NOT_GLOBAL_RECOLOR * p)));
  } else if (curMoveType == BD_TRIANGLE_DEATH) {
    // Reverse of: BD_TRIANGLE_BIRTH
    int n = state.numBoundaryEdges() - 2;
    double a = CGAL::to_double(squared_distance(*(btd.prevBoundaryVertex),
            *(btd.nextBoundaryVertex)));
    double b = CGAL::to_double(squared_distance(btd.up, btd.vp));
    bool couldHaveBeenLocal = (b <= castingBoxRadius * castingBoxRadius);
    return ln(ZETA_BD_NOT_MR * 
        ZETA_BIRTH_NOT_DEATH *
        (1.0 - ZETA_INT_NOT_BOUNDARY) *
        ZETA_BOUND_NOT_CORNER *
        (
         // Global probability
         ((1.0 - ZETA_LOCAL_NOT_GLOBAL_BD_TRI_BIRTH) *
    2.0 / 
    (double(n) * a * b))
         +
         // Local probability
         (couldHaveBeenLocal ? 
    (ZETA_LOCAL_NOT_GLOBAL_BD_TRI_BIRTH *
     2.0 / 
     (double(n) * 
      sqrt(a) * 2.0 * castingBoxRadius * b)) 
    :
    0.0)));
  } else if (curMoveType == MOVE_INT_VERTEX) {
    if (slide_not_move) {
      // Reverse of: SLIDE_INT_VERTEX (= MOVE_INT_VERTEX + flag)
      Coloring::VertexHandle nbr = 
  siv.next ? siv.vertex->getNextVertex() : siv.vertex->getPrevVertex();
      double new_length = 
  sqrt(CGAL::to_double(CGAL::squared_distance(siv.newLoc, 
                nbr->point())));
      return ln((1.0 - ZETA_BD_NOT_MR) * 
    ZETA_MOVE_NOT_RECOLOR *
    ZETA_INT_NOT_BOUNDARY *
    ZETA_SLIDE_NOT_MOVE_INT_VERTEX /
    (double(state.numVertices(Coloring::Vertex::INTERIOR)) /
     (3.0 * new_length)));
    } else {
      // Reverse of: MOVE_INT_VERTEX
      return ln((1.0 - ZETA_BD_NOT_MR) * 
    ZETA_MOVE_NOT_RECOLOR *
    ZETA_INT_NOT_BOUNDARY *
    (1.0 - ZETA_SLIDE_NOT_MOVE_INT_VERTEX) /
    (double(state.numVertices(Coloring::Vertex::INTERIOR)) *
     castingArea()));
    }
  } else
    return CN94Proposal::rll(state);
}

---