point_obs.cpp

Go to the documentation of this file.

#include <fstream>
#include <CGAL/IO/Qt_widget.h>
#include "point_obs.hpp"
#include "query_point.hpp"

using namespace Arak;
using namespace Arak::Geometry;

ArakPosteriorGaussianObs::ArakPosteriorGaussianObs(const ArakProcess& prior,
               const Arak::Util::PropertyMap& props) 
  : ArakProcess(prior.getColoring()), prior(prior), lp(0.0) {
  using namespace Arak::Util;
  assert(parse(getp(props, "arak.gaussian_obs.white_mean"), whiteMean));
  assert(parse(getp(props, "arak.gaussian_obs.white_var"), whiteVar));
  assert(parse(getp(props, "arak.gaussian_obs.black_mean"), blackMean));
  assert(parse(getp(props, "arak.gaussian_obs.black_var"), blackVar));
  const std::string& datafile_path = 
    getp(props, "arak.gaussian_obs.data_file");
  std::ifstream in(datafile_path.data());
  assert(in.good());
  int n;
  in >> n;
  std::vector<Geometry::Point> pts;
  points.reserve(n);
  obs.reserve(n);
  potentials.reserve(n);
  for (int i = 0; i < n; i++) {
    Geometry::Point p;
    double o;
    in >> p;
    in >> o;
    points.push_back(p);
    obs.push_back(o);
    double blackPotential = 
      (o - blackMean) * (o - blackMean) / 2.0 * blackVar;
    double whitePotential = 
      (o - whiteMean) * (o - whiteMean) / 2.0 * whiteVar;
    potentials.push_back(new QueryPointPotential(lp,
             blackPotential,
             whitePotential));
  }
  assert(!in.bad());
  QueryPointIndex index(points, potentials);
  assert(lp == 0.0);
  c.addQueryPoints(index);
}

double ArakPosteriorGaussianObs::logMeasure() const {
  return prior.logMeasure();
}

double ArakPosteriorGaussianObs::potential() const {
  return prior.potential() + lp;
}

ArakPosteriorGaussianObs::~ArakPosteriorGaussianObs() { 
  for (int i = 0; i < int(potentials.size()); i++)
    delete potentials[i];
}

void ArakPosteriorGaussianObs::visualize(CGAL::Qt_widget& widget) const {
  // Compute the mapping from observations to colors.
  double white = std::min<double>(blackMean - 3.0 * sqrt(blackVar),
          whiteMean - 3.0 * sqrt(whiteVar));
  double black = std::max<double>(blackMean + 3.0 * sqrt(blackVar),
          whiteMean + 3.0 * sqrt(whiteVar));
  for (int i = 0; i < int(points.size()); i++) {
    CGAL::Color c;
    if (obs[i] <= white)
      c = CGAL::WHITE;
    else if (obs[i] >= black) 
      c = CGAL::BLACK;
    else {
      double w = (obs[i] - white) / (black - white);
      int v = 255 - int(w * 255.0);
      c = CGAL::Color(v, v, v);
    }
    widget << c << points[i];
  }
}

ArakPosteriorBernoulliObs::ArakPosteriorBernoulliObs(const ArakProcess& prior,
                 const Arak::Util::PropertyMap& props) 
  : ArakProcess(prior.getColoring()), prior(prior), lp(0.0) {
  using namespace Arak::Util;
  assert(parse(getp(props, "arak.bernoulli_obs.white_prob"), whiteProb));
  assert(parse(getp(props, "arak.bernoulli_obs.black_prob"), blackProb));
  const std::string& datafile_path = 
    getp(props, "arak.bernoulli_obs.data_file");
  std::ifstream in(datafile_path.data());
  assert(in.good());
  int n;
  in >> n;
  points.reserve(n);
  obs.reserve(n);
  potentials.reserve(n);
  for (int i = 0; i < n; i++) {
    Geometry::Point p;
    int o;
    in >> p;
    in >> o;
    points.push_back(p);
    obs.push_back(o == 0);
    double blackPotential = (o == 0) ? -ln(1.0 - blackProb) : -ln(blackProb);
    double whitePotential = (o == 0) ? -ln(whiteProb) : -ln(1.0 - whiteProb);
    potentials.push_back(new QueryPointPotential(lp,
             blackPotential,
             whitePotential));
  }
  assert(!in.bad());
  QueryPointIndex index(points, potentials);
  c.addQueryPoints(index);
}

double ArakPosteriorBernoulliObs::logMeasure() const {
  return prior.logMeasure();
}

double ArakPosteriorBernoulliObs::potential() const {
  return prior.potential() + lp;
}

ArakPosteriorBernoulliObs::~ArakPosteriorBernoulliObs() { 
  for (int i = 0; i < int(potentials.size()); i++)
    delete potentials[i];
}

void ArakPosteriorBernoulliObs::visualize(CGAL::Qt_widget& widget) const {
  for (int i = 0; i < int(points.size()); i++) {
    CGAL::Color c;
    if (obs[i])
      c = CGAL::BLACK;
    else 
      c = CGAL::RED;
    widget << c << points[i];
  }
}

---