mhrwtasks.h

Go to the documentation of this file.

/* This file is part of the Tomographer project, which is distributed under the
 * terms of the MIT license.
 *
 * The MIT License (MIT)
 *
 * Copyright (c) 2015 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef TOMOGRAPHER_MHWALKER_TASKS_H
#define TOMOGRAPHER_MHWALKER_TASKS_H

#include <string>
#include <limits>
#include <random>

#include <tomographer/tools/fmt.h>
#include <tomographer/mhrw.h>
#include <tomographer/multiproc.h> // StatusReport Base


namespace Tomographer {


namespace MHRWTasks
{

  template<typename CountIntType_ = unsigned int, typename RealType_ = double>
  struct CDataBase
  {
    CDataBase(int base_seed_ = 0)
      : base_seed(base_seed_)
    {
    }

    typedef CountIntType_ CountIntType;
    typedef RealType_ RealType;

    CountIntType n_sweep;
    CountIntType n_therm;
    CountIntType n_run;
    RealType step_size;

    int base_seed;
  };

  template<typename MHRWStatsCollectorResultType_, typename CountIntType>
  struct MHRandomWalkTaskResult
  {
    typedef MHRWStatsCollectorResultType_ MHRWStatsCollectorResultType;
    
    template<typename MHRWStatsCollectorResultTypeInit, typename MHRandomWalkType>
    MHRandomWalkTaskResult(MHRWStatsCollectorResultTypeInit && stats_collector_result_,
                           const MHRandomWalkType & mhrandomwalk)
      : stats_collector_result(std::forward<MHRWStatsCollectorResultTypeInit>(stats_collector_result_)),
        n_sweep(mhrandomwalk.n_sweep()),
        n_therm(mhrandomwalk.n_therm()),
        n_run(mhrandomwalk.n_run()),
        step_size(mhrandomwalk.step_size()),
        acceptance_ratio(mhrandomwalk.acceptance_ratio())
    {
    }
    
    const MHRWStatsCollectorResultType stats_collector_result;
    
    const CountIntType n_sweep;
    const CountIntType n_therm;
    const CountIntType n_run;
    const double step_size;
    const double acceptance_ratio;
  };

  template<typename MHRandomWalkTaskCData,
           typename Rng = std::mt19937>
  struct MHRandomWalkTask
  {
    typedef typename MHRandomWalkTaskCData::CountIntType CountIntType;

    typedef MHRandomWalkTaskResult<typename MHRandomWalkTaskCData::MHRWStatsCollectorResultType,
                                   CountIntType> Result;

    struct StatusReport : public MultiProc::StatusReport
    {
      StatusReport(double fdone = 0.0, const std::string & msg = std::string(),
                   CountIntType kstep_ = 0, CountIntType n_sweep_ = 0,
                   CountIntType n_therm_ = 0, CountIntType n_run_ = 0,
                   double acceptance_ratio_ = 0.0)
        : MultiProc::StatusReport(fdone, msg),
          kstep(kstep_),
          n_sweep(n_sweep_),
          n_therm(n_therm_),
          n_run(n_run_),
          acceptance_ratio(acceptance_ratio_),
          n_total_iters(n_sweep*(n_therm+n_run))
      {
      }
      CountIntType kstep;
      CountIntType n_sweep;
      CountIntType n_therm;
      CountIntType n_run;
      double acceptance_ratio;
      CountIntType n_total_iters;
    };
    typedef StatusReport StatusReportType;

  private:
    typename Rng::result_type _seed;

    Result * result;

  public:

    static inline int get_input(int k, const MHRandomWalkTaskCData * pcdata)
    {
      return pcdata->base_seed + k;
    }

    template<typename LoggerType>
    MHRandomWalkTask(int inputseed, const MHRandomWalkTaskCData * /*pcdata*/, LoggerType & logger)
      : _seed(inputseed), result(NULL)
    {
      logger.longdebug("MHRandomWalkTask", "() inputseed=%d", inputseed);
    }

    ~MHRandomWalkTask()
    {
      if (result != NULL) {
        delete result;
      }
    }


    template<typename LoggerType, typename TaskManagerIface>
    inline void run(const MHRandomWalkTaskCData * pcdata, LoggerType & logger,
                    TaskManagerIface * tmgriface)
    {
      Rng rng(_seed); // seeded random number generator

      // the user's stats collector
      auto stats = pcdata->createStatsCollector(logger);
      typedef decltype(stats) MHRWStatsCollectorType;

      // our own "stats collector" which checks if we need to send a status report back
      typedef StatusReportCheck<TaskManagerIface, MHRWStatsCollectorType> OurStatusReportCheck;
      OurStatusReportCheck statreportcheck(this, &stats, tmgriface);
      
      typedef MultipleMHRWStatsCollectors<MHRWStatsCollectorType, OurStatusReportCheck>
        OurStatsCollectors;
      OurStatsCollectors ourstatscollectors(stats, statreportcheck);

      auto mhwalker = pcdata->createMHWalker(rng, logger);
      typedef decltype(mhwalker) MHWalkerType;

      MHRandomWalk<Rng, MHWalkerType, OurStatsCollectors, LoggerType, CountIntType> rwalk(
          // MH random walk parameters
          pcdata->n_sweep,
          pcdata->n_therm,
          pcdata->n_run,
          pcdata->step_size,
          // the MHWalker
          mhwalker,
          // our stats collectors
          ourstatscollectors,
          // a random number generator
          rng,
          // and a logger
          logger
          );
      
      rwalk.run();

      result = new Result(stats.getResult(), rwalk);
    }

    inline const Result & getResult() const
    {
      return *result;
    }

  private:
    template<typename TaskManagerIface, typename MHRWStatsCollectorType>
    struct StatusReportCheck
    {
      StatusReportCheck(MHRandomWalkTask * mhrwtask_, MHRWStatsCollectorType * stats_, TaskManagerIface *tmgriface_)
        : mhrwtask(mhrwtask_), stats(stats_), tmgriface(tmgriface_)
      { }

      MHRandomWalkTask *mhrwtask;
      MHRWStatsCollectorType *stats;
      TaskManagerIface *tmgriface;

      inline void init() { }
      inline void thermalizing_done() { }
      inline void done() { }

      template<typename PointType, typename FnValueType, typename MHRandomWalk>
      inline void raw_move(
          CountIntType k, bool is_thermalizing, bool, bool, double, const PointType &, FnValueType,
          const PointType &, FnValueType, MHRandomWalk & rw
          )
      {
        // see if we should provide a status report
        //        fprintf(stderr, "StatusReportCheck::raw_move(): testing for status report requested!\n");
        if (tmgriface->status_report_requested()) {
          // prepare & provide status report
          CountIntType totiters = rw.n_sweep()*(rw.n_therm()+rw.n_run());
          double fdone = (double)k/totiters;
          double accept_ratio = std::numeric_limits<double>::quiet_NaN();
          bool warn_accept_ratio = false;
          if (rw.has_acceptance_ratio()) {
            accept_ratio = rw.acceptance_ratio();
            warn_accept_ratio = (accept_ratio > 0.35 || accept_ratio < 0.2);
          }
          std::string msg = Tools::fmts(
              "%s %lu/(%lu=%lu*(%lu+%lu)) : %5.2f%% done  [%saccept ratio=%.2f%s]",
              ( ! is_thermalizing
                ? "iteration"
                : "[therm.] "),
              (unsigned long)k, (unsigned long)totiters, (unsigned long)rw.n_sweep(),
              (unsigned long)rw.n_therm(), (unsigned long)rw.n_run(),
              fdone*100.0,
              warn_accept_ratio ? "!!** " : "",
              accept_ratio,
              warn_accept_ratio ? " **!!" : ""
              );
          typedef MHRWStatsCollectorStatus<MHRWStatsCollectorType> MHRWStatsCollectorStatusType;
          if (MHRWStatsCollectorStatusType::CanProvideStatus) {
            std::string nlindent = "\n    ";
            msg += nlindent;
            std::string s = MHRWStatsCollectorStatusType::getStatus(stats);
            for (std::size_t j = 0; j < s.size(); ++j) {
              if (s[j] == '\n') {
                msg += nlindent;
              } else {
                msg += s[j];
              }
            }
          }
          tmgriface->submit_status_report(StatusReport(fdone, msg, k, rw.n_sweep(), rw.n_therm(),
                                                       rw.n_run(), accept_ratio));
        }
        //        fprintf(stderr, "StatusReportCheck::raw_move(): done\n");
      }

      template<typename PointType, typename FnValueType, typename MHRandomWalk>
      inline void process_sample(CountIntType, CountIntType, const PointType &, FnValueType, MHRandomWalk &)
      {
      }
    };
  };



} // MHWalkerTasks



} // namespace Tomographer



#endif