sampler.cpp

Go to the documentation of this file.

#include <csignal>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <ctime>
#include <math.h>
#include <qthread.h>
#include <gzstream.h>
#include "coloring.hpp"
#include "random.hpp"
#include "gui.hpp"
#include "point_obs.hpp"
#include "range_obs.hpp"
#include "sonar_obs.hpp"
#include "coloring_mcmc.hpp"
#include "cn94.hpp"
#include "query_point.hpp"
#include "estimation.hpp"
#include "parsecl.hpp"
#include "properties.hpp"
#include "verification.hpp"

using std::cout; 
using std::endl;
using std::istream;
using std::ifstream;
using std::ofstream;
using namespace Arak;
using namespace Arak::Util;

const char* default_props = "\n\
arak.proposal = cn94 \n\
arak.cn94.zeta.bd_not_mr = 0.5 \n\
arak.cn94.zeta.move_not_recolor = 0.5 \n\
arak.cn94.zeta.birth_not_death = 0.5 \n\
arak.cn94.zeta.int_not_boundary = 0.68 \n\
arak.cn94.zeta.bound_not_corner = 0.9 \n\
arak.cn94.zeta.ib_vertex_not_triangle = 0.6 \n\
arak.modified_cn94.zeta.local_not_global_recolor = 0.9 \n\
arak.modified_cn94.zeta.local_not_global_bd_tri_birth = 0.5 \n\
# arak.sampler.coloring_path = /tmp/coloring.out \n\
# arak.sampler.map_coloring_path = /tmp/map-coloring.out \n\
# arak.sampler.estimates_path = /tmp/estimates.out.gz \n\
# arak.sampler.flush_interval = 600.0 \n\
# arak.sampler.samples_path = /tmp/samples.bin \n\
# arak.sampler.record_stride = 1000 \n\
arak.scale = 1.0 \n\
arak.coloring.xmin = 0.0 \n\
arak.coloring.ymin = 0.0 \n\
arak.coloring.width = 1.0 \n\
arak.coloring.height = 1.0 \n\
arak.coloring.rows = 35 \n\
arak.coloring.cols = 35 \n\
arak.mcmc.estimation.stride = 100 \n\
arak.mcmc.estimation.burn_in = 1000 \n\
arak.mcmc.estimation.grid_size = 10000 \n\
";

// Signal handling for interrupts
volatile std::sig_atomic_t interrupted = 0;
extern "C" void sigint(int sig) { interrupted = 1; }

class MCMCThread : public QThread {

 protected:

  ArakMarkovChain& chain;

  const GridColorEstimator& estimates;

  CriticalStatsEstimator* stats;

  QMutex& mutex;

  std::string coloring_path;

  std::string estimates_path;

  std::string map_coloring_path;

  double map_log_likelihood;

  std::ofstream* samples_stream;

  int record_stride;

  double flush_interval;

  void flush() {
    if (!coloring_path.empty()) {
      ofstream clr_out(coloring_path.data());
      mutex.lock();
      clr_out << chain.getState();
      mutex.unlock();
      clr_out.close();
    }
    if (!estimates_path.empty()) {
      ogzstream est_out(estimates_path.data());
      mutex.lock();
      est_out << estimates;
      mutex.unlock();
      est_out.close();
    }
  }

  void processSample() {
    mutex.lock();
    if (stats != NULL)
      stats->update(chain.getState());
    // Record the current sample if needed.
    if ((samples_stream != NULL) &&
  (chain.getNumSteps() % (record_stride + 1) == 0)) {
      chain.getState().writeBinary(*samples_stream);
    }
    double cur_log_likelihood = chain.getLogLikelihood();
    if ((cur_log_likelihood > map_log_likelihood)
  &&
  (!map_coloring_path.empty())
  &&
  // Make sure we're past the burn-in phase to avoid too many
  // initial writes.
  (chain.getNumSamples() > 1)) {
      map_log_likelihood = cur_log_likelihood;
      ofstream clr_out(map_coloring_path.data());
      clr_out << chain.getState();
      clr_out.close();
    }
    mutex.unlock();
  }

  void run() {
    // The number of samples between flushes.  This is determined
    // automatically.
    unsigned long int interval = 0;
    std::time_t last_flush = std::time(NULL);
    while (true) {
      if (interrupted != 0) goto report;
      mutex.lock();
      chain.advance();
      mutex.unlock();
      processSample();
      std::time_t cur_time = std::time(NULL);
      if (std::difftime(cur_time, last_flush) > flush_interval) 
  break;
      else
  interval++;
    }
    // We have determined the interval and are ready for the first flush.
    while (true) {
      // Write the current sample and point color estimates.
      last_flush = std::time(NULL);
      flush();
      std::string timestr(std::asctime(std::localtime(&last_flush)), 23);
      std::cerr << "[" << timestr << "] " 
    << "Preliminary estimates written to " 
    << coloring_path << "/"
    << estimates_path << "." << endl;
      // Continue sampling.
      for (unsigned long int i = 0; i < interval; i++) {
  if (interrupted != 0) goto report;
  mutex.lock();
  chain.advance();
  mutex.unlock();
  processSample();
      }
    }
  report:
    if (stats != NULL) {
      std::cout << std::endl << "VERIFICATION STATISTICS:" << std::endl;
      stats->report(chain.getState().boundary(),
        chain.getDistribution().scale(),
        std::cout);
      std::cout << std::endl;
    }
    if (samples_stream != NULL) samples_stream->close();
  }

 public:

  MCMCThread(ArakMarkovChain& chain,
       const GridColorEstimator& estimates,
       QMutex& mutex,
       const Arak::Util::PropertyMap& props) : 
    chain(chain), 
    estimates(estimates),
    mutex(mutex),
    map_log_likelihood(-std::numeric_limits<double>::infinity()) 
  { 
    // Parse the properties.
    if (hasp(props, "arak.sampler.coloring_path"))
      coloring_path = getp(props, "arak.sampler.coloring_path");
    else
      coloring_path = "";
    if (hasp(props, "arak.sampler.map_coloring_path"))
      map_coloring_path = getp(props, "arak.sampler.map_coloring_path");
    else 
      map_coloring_path = "";
    if (hasp(props, "arak.sampler.estimates_path"))
      estimates_path = getp(props, "arak.sampler.estimates_path");
    else
      estimates_path = "";
    if (hasp(props, "arak.sampler.flush_interval"))
      assert(parse(getp(props, "arak.sampler.flush_interval"), 
       flush_interval));
    else
      flush_interval = 600.0; // 10 minutes
    if (hasp(props, "arak.sampler.verify")) {
      int k;
      assert(parse(getp(props, "arak.sampler.verify"), k));
      stats = new CriticalStatsEstimator(k, chain.getState());
    } else
      stats = NULL;
    if (hasp(props, "arak.sampler.samples_path")) {
      const std::string& samples_path = 
  getp(props, "arak.sampler.samples_path");
      samples_stream = new std::ofstream(samples_path.data(), 
           std::ios::binary |
           std::ios::out);
      if (hasp(props, "arak.sampler.record_stride"))
  assert(parse(getp(props, "arak.sampler.record_stride"), 
         record_stride));
    } else
      samples_stream = NULL;
  }

  ~MCMCThread() {
    if (stats != NULL) delete stats;
    if (samples_stream != NULL) delete samples_stream;
  }
};

int sampler(int argc, 
      char** argv,
      bool visualize,
      bool toolbar,
      double refreshRateHz,
      const Arak::Util::PropertyMap& props) {
  // Form the initial coloring.
  Coloring c(props);
  if (hasp(props, "arak.sampler.init_state_path")) {
    const std::string& path = getp(props, "arak.sampler.init_state_path");
    ifstream in(path.data());
    in >> c;
    if (!in.good()) {
      std::cerr << "Cannot read initial coloring from " << path << std::endl;
      exit(1);
    }
  }
  // Form the prior process.
  ArakProcess* prior = NULL;
  using namespace Arak::Util;
  if (hasp(props, "arak.prior")) {
    const std::string& prior_type = getp(props, "arak.prior");
    if (prior_type == "standard") 
      prior = new ArakPrior(c, props);
    else
      std::cerr << "Unrecognized prior type: " 
    << prior_type << std::endl;
  } else 
    prior = new ArakPrior(c, props);
  // Form the posterior process.
  ArakProcess* process = NULL;
  using namespace Arak::Util;
  if (hasp(props, "arak.obs_type")) {
    const std::string& obs_type = getp(props, "arak.obs_type");
    if (obs_type == "gaussian") 
      process = new ArakPosteriorGaussianObs(*prior, props);
    else if (obs_type == "bernoulli") 
      process = new ArakPosteriorBernoulliObs(*prior, props);
    else if (obs_type == "range") 
      process = new ArakPosteriorRangeObs(*prior, props);
    else if (obs_type == "sonar") 
      process = new ArakPosteriorSonarObs(*prior, props);
    else
      std::cerr << "Unrecognized observation type: " 
    << obs_type << std::endl;
  } else 
    process = prior;
  
  // Form the proposal distribution.
  CN94Proposal* proposal = NULL;
  const std::string& proposal_type = getp(props, "arak.proposal");
  if (proposal_type == "cn94") 
    proposal = new CN94Proposal(props);
  else if (proposal_type == "modified_cn94") 
    proposal = new ModifiedCN94Proposal(props);
  else {
    std::cerr << "Unrecognized proposal distribution: "
        << proposal_type << std::endl;
    exit(1);
  }

  // Choose the random number generator.
  Arak::Util::Random& random = Arak::Util::default_random;
  std::cout << "PRNG: " << random << std::endl;

  // Form the Markov chain.
  ArakMarkovChain* chain;
  if (hasp(props, "arak.sampler.chain_type")) {
    const std::string& chain_type = getp(props, "arak.sampler.chain_type");
    if (chain_type == "standard") 
      chain = new ArakMarkovChain(*process, *proposal, c, 
          props, default_random);
    else if (chain_type == "annealed") 
      chain = new AnnealedArakMarkovChain(*process, *proposal, c, 
            props, default_random);
    else if (chain_type == "hill-climbing") {
      chain = new HillClimbingArakMarkovChain(*process, *proposal, c, 
                props, default_random);
      std::cerr << "Warning: hill-climbing in an Arak process is usually a bad"
    << std::endl
    << "         idea; try annealing instead."
    << std::endl;
    } else {
      std::cerr << "Unrecognized Markov chain type: "
    << chain_type << std::endl;
      exit(1);
    }
  } else
    chain = new ArakMarkovChain(*process, *proposal, c, 
        props, default_random);

  // Form the point set color estimator.
  GridColorEstimator estimator(*chain, props);
  
  // Create the sampling thread and the GUI application.
  QMutex mutex;
  MCMCThread thread(*chain, estimator, mutex, props);

  // Start the inference thread and transfer control to the GUI.
  std::time_t start = std::time(NULL);
  thread.start();
  int result = 0;

  // If visualization was requested, start the GUI.
  if (visualize) {
    QApplication *app = InitGui(argc, argv, 
        480, 480, 
        c, *process, estimator,
        mutex,
        toolbar, 
        refreshRateHz);
    result = app->exec();
    interrupted = 1;
  } else {
    if (std::signal(SIGINT, sigint) == SIG_ERR)
      std::cerr << "Warning: cannot set signal handler." << endl;
    else
      cout << "Signal INT (press Ctrl-C) to stop sampling..." << endl;
  }

  // Wait for the sampler to stop.
  thread.wait();
  std::time_t end = std::time(NULL);

  // Report the results.
  cout << "Number of samples: " << chain->getNumSteps() 
       << " (" << double(chain->getNumSteps()) / difftime(end, start)
       << " samples per second)" << endl;
  cout << "Acceptance ratio: " << chain->acceptance() 
       << " (" << double(chain->getNumMoves()) / difftime(end, start)
       << " moves per second)" << endl;
  proposal->writeStatistics(cout);

  // Write the final coloring and point color estimates.
  if (hasp(props, "arak.sampler.coloring_path")) {
    std::string coloring_path = getp(props, "arak.sampler.coloring_path");
    ofstream clr_out(coloring_path.data());
    clr_out << chain->getState();
    clr_out.close();
    cout << "Final coloring written to " << coloring_path << "." << endl;
  }
  if (hasp(props, "arak.sampler.estimates_path")) {
    std::string estimates_path = getp(props, "arak.sampler.estimates_path");
    ogzstream est_out(estimates_path.data());
    est_out << estimator;
    est_out.close();
    cout << "Final estimates written to " << estimates_path << "." << endl;
  }

  // Deallocate the chain, prior, process, and proposal.
  delete chain;
  if (prior != process)
    delete process;
  delete prior;
  delete proposal;
  
  // Return the result of the GUI.
  return result;
}

int main(int argc, char** argv) {
  using namespace Arak::Util;
  // Parse the command line.
  CommandLine cl;
  CommandLine::Option visualize("-v", "--visualize", 
        "visualizes the state of the sampler");
  cl.add(visualize);
  CommandLine::Parameter<double> 
    refresh("-r", "--refresh-rate", 
      "the rate (in Hz) that the visualization is updated", 1.0);
  cl.add(refresh);
  CommandLine::Option toolbar("-t", "--toolbar", 
            "uses the toolbar in the visualization");
  cl.add(toolbar);
  CommandLine::MultiParameter<std::string>
    propFiles("-p", "--prop-file", "specifies a property file");
  cl.add(propFiles);
  CommandLine::MultiParameter<std::string>
    propDefs("-D", "--define-prop", "defines a property value");
  cl.add(propDefs);

  if (!cl.parse(argc, argv, std::cerr)) {
    cl.printUsage(std::cerr);
    exit(1);
  }

  // Parse the default properties.
  Arak::Util::PropertyMap properties;
  {
    std::string defaultProps(default_props);
    std::istringstream in(defaultProps);
    in >> properties;
  }

  // Parse the properties files.
  const std::vector<std::string>& propFilePaths = propFiles.values();
  for (std::vector<std::string>::const_iterator it = propFilePaths.begin();
       it != propFilePaths.end(); it++) {
    std::string path = *it;
    std::ifstream in(path.data());
    in >> properties;
  }
  // Parse the overriding property definitions.
  const std::vector<std::string>& propDefinitions = propDefs.values();
  for (std::vector<std::string>::const_iterator it = propDefinitions.begin();
       it != propDefinitions.end(); it++) {
    std::string def = *it;
    std::istringstream in(def);
    in >> properties;
  }

  std::cout << "Properties:" << std::endl << properties << std::endl;

  // Invoke the program.
  return sampler(argc, argv, 
     visualize.supplied(),
     toolbar.supplied(),
     refresh.value(),
     properties);
}

---