estimation.hpp

Go to the documentation of this file.

#ifndef _ESTIMATION_HPP
#define _ESTIMATION_HPP

#include <iostream>
#include <iomanip>
#include "query_point.hpp"
#include "coloring_mcmc.hpp"

#define POINT_SET_ESTIMATE_HEADER "\
#######################################################################\n\
#                                                                      \n\
# This file represents estimates of the colors of points under an      \n\
# Arak process.  The first line below consists of the number of samples\n\
# used to compute the estimates (needed to merge estimates) and then   \n\
# the number of point estimates.  Each line thereafter represents an   \n\
# point color estimate: the first two numbers are the X and Y          \n\
# coordinates of the point, and the third number is the estimated      \n\
# probability that the point is colored black.                         \n\
#                                                                      \n\
#######################################################################\n"

#define GRID_ESTIMATE_HEADER "\
#######################################################################\n\
#                                                                      \n\
# This file represents estimates of the colors of a grid of points     \n\
# under an Arak process.  The first line below consists of the number  \n\
# of samples used to compute the estimates (needed to merge estimates).\n\
# The second line has the number of rows and columns in the grid.  The \n\
# third line gives the boundary of the grid.  Then the estimates are   \n\
# written out as a matrix (one line per row) in reverse-row order (so  \n\
# the first element corresponds to the upper left corner of the grid.  \n\
# Each estimate is the probability the cell's center point is black.   \n\
#                                                                      \n\
#######################################################################\n"

namespace Arak {

  class PointColorEstimator : public QueryPointListener {

    friend class PointSetColorEstimator;
    friend class GridColorEstimator;

  protected:

    const ArakMarkovChain& chain;

    const Geometry::Point point;

    unsigned long int updated;

    Color color;

    unsigned long int numBlack;

    unsigned long int numWhite;

    void update(Color color, unsigned long int time) {
      assert(updated <= time);
      // This is the number of samples that the point stayed its
      // previous color.
      unsigned long int duration = time - updated;
      if (this->color == BLACK)
  numBlack += duration;
      else if (this->color == WHITE)
  numWhite += duration;
      this->color = color;
      updated = time;
    }

  public:
    
    PointColorEstimator(const ArakMarkovChain& chain,
      const Geometry::Point& point,
      Color color) : 
      QueryPointListener(), chain(chain), point(point), 
      updated(chain.getNumSamples()), color(color), 
      numBlack(0), numWhite(0) { }

    virtual ~PointColorEstimator() { }

    virtual void recolor(Color color) {
      update(color, chain.getNumSamples());
    }

    double estimate() const { 
      unsigned long int time = chain.getNumSamples();
      unsigned long int curNumBlack = 
  numBlack + ((color == BLACK) ? time - updated : 0);
      unsigned long int curNumWhite = 
  numWhite + ((color == BLACK) ? 0 : time - updated);
      if ((curNumBlack == 0) && (curNumWhite == 0))
  return 0.5;
      else
  return (double(curNumBlack) / 
    (double(curNumBlack) + double(curNumWhite)));
    }

  };

  class PointSetColorEstimator : public QueryPointListenerFactory {

  protected:

    const ArakMarkovChain& chain;

    std::list<PointColorEstimator*> estimators;

    QueryPointIndex* index;

  public:

    PointSetColorEstimator(ArakMarkovChain& chain,
         const Arak::Util::PropertyMap& props);

    virtual ~PointSetColorEstimator() {
      while (!estimators.empty()) {
  PointColorEstimator* estimator = estimators.front();
  estimators.pop_front();
  delete estimator;
      }
      // TODO: stop listening to coloring
      delete index;
    }

    virtual PointColorEstimator* create(const QueryPoint& p) {
      PointColorEstimator* estimator = 
  new PointColorEstimator(chain, p, p.color());
      estimators.push_front(estimator);
      return estimator;
    }

    virtual void visualize(CGAL::Qt_widget& widget) const;

    template<typename charT, typename traits>
    void write(std::basic_ostream<charT,traits>& out) const {
      // Write some comment information about the format of the file.
      out << POINT_SET_ESTIMATE_HEADER;
      // Write out the number of samples and the number of points.
      out << chain.getNumSamples() << " " 
    << estimators.size() << std::endl;
      // Now write out each point color estimate.
      out << std::scientific << std::setw(16) << std::setprecision(12);
      typedef std::list<PointColorEstimator*>::const_iterator Iterator;
      for (Iterator it = estimators.begin(); it != estimators.end(); it++) {
  const PointColorEstimator* estimator = *it;
  out << estimator->point << " " << estimator->estimate() << std::endl;
      }
    }
  }; // End of class: PointSetColorEstimator

  inline CGAL::Qt_widget& operator<<(CGAL::Qt_widget& widget, 
             PointSetColorEstimator& p) {
    p.visualize(widget);
    return widget;
  }

  template<typename charT, typename traits>
  std::basic_ostream<charT,traits>& 
  operator<<(std::basic_ostream<charT,traits>& out,
       const PointSetColorEstimator& p) {
    p.write(out);
    return out;
  }

  class GridColorEstimator : public Coloring::Listener {

  protected:

    const ArakMarkovChain& chain;

    std::list<PointColorEstimator*> estimators;

    typedef Grid<PointColorEstimator*> EstimatorGrid;

    EstimatorGrid* grid;

  public:

    GridColorEstimator(ArakMarkovChain& chain, 
           const Arak::Util::PropertyMap& props);

    virtual ~GridColorEstimator() {
      while (!estimators.empty()) {
  PointColorEstimator* estimator = estimators.front();
  estimators.pop_front();
  delete estimator;
      }
      delete grid;
    }

    virtual void recolored(const Geometry::Point& a,
         const Geometry::Point& b,
         const Geometry::Point& c);
    
    virtual void recolored(const Geometry::Point& a,
         const Geometry::Point& b,
         const Geometry::Point& c,
         const Geometry::Point& d) {
      recolored(a, b, c);
      recolored(a, d, c);
    }
    
    virtual void visualize(CGAL::Qt_widget& widget) const;

    template<typename charT, typename traits>
    void write(std::basic_ostream<charT,traits>& out) const {
      // Write some comment information about the format of the file.
      out << GRID_ESTIMATE_HEADER;
      // Write out the number of samples and the number of points.
      out << chain.getNumSamples() << std::endl;
      // Write out the number of rows and columns.
      out << grid->numRows() << " " << grid->numCols() << std::endl;
      // Write out the boundary of the grid.
      out << grid->boundary() << std::endl;
      // Now write out each point color estimate.
      out << std::scientific << std::setw(16) << std::setprecision(12);
      for (int i = grid->numRows() - 1; i >= 0; i--) {
  for (int j = 0; j < grid->numCols(); j++) {
    const EstimatorGrid::Cell& cell = grid->getCell(i, j);
    const PointColorEstimator* pce = *(cell.getItemList().begin());
    out << pce->estimate() << " ";
  }
  out << std::endl;
      }
    }
  }; // End of class: GridColorEstimator

  inline CGAL::Qt_widget& operator<<(CGAL::Qt_widget& widget, 
             GridColorEstimator& p) {
    p.visualize(widget);
    return widget;
  }

  template<typename charT, typename traits>
  std::basic_ostream<charT,traits>& 
  operator<<(std::basic_ostream<charT,traits>& out,
       const GridColorEstimator& p) {
    p.write(out);
    return out;
  }

}

#endif

---