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) 2016 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
 * Copyright (c) 2017 Caltech, Institute for Quantum Information and Matter, 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 <climits>

#include <string>
#include <limits>
#include <random>
#include <sstream>
#include <stdexcept>

#include <boost/serialization/serialization.hpp>
#include <boost/serialization/vector.hpp>

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


namespace Tomographer {
namespace MHRWTasks {




template<typename MHWalkerParams_, typename IterCountIntType_ = int,
         typename RngSeedType_ = std::mt19937::result_type>
struct TOMOGRAPHER_EXPORT CDataBase
{
  typedef IterCountIntType_ IterCountIntType;
  typedef MHWalkerParams_ MHWalkerParams;

  typedef RngSeedType_ RngSeedType;

  typedef MHRWParams<MHWalkerParams, IterCountIntType> MHRWParamsType;


  template<typename MHRWParamsType>
  CDataBase(MHRWParamsType&& p, RngSeedType base_seed_)
    : mhrw_params(std::forward<MHRWParamsType>(p)),
      base_seed(base_seed_),
      task_seeds()
  {
  }

  template<typename MHRWParamsType>
  CDataBase(MHRWParamsType&& p, std::vector<RngSeedType> task_seeds_)
    : mhrw_params(std::forward<MHRWParamsType>(p)),
      base_seed(0),
      task_seeds(std::move(task_seeds_))
  {
  }

  CDataBase() : mhrw_params(), base_seed(), task_seeds() { }


  MHRWParamsType mhrw_params;

  RngSeedType base_seed;

  std::vector<RngSeedType> task_seeds;


  template<typename TaskNoCountIntType>
  inline RngSeedType getTaskInput(TaskNoCountIntType k) const
  {
    if (task_seeds.size() > 0) {
      // specific task seeds given (e.g. from random device)
      if (k < 0 || (std::size_t)k >= task_seeds.size()) {
        throw std::out_of_range(streamstr("getTaskInput(): k="<<k<<" out of range; seed list size is = "
                                          << task_seeds.size()));
      }
      return task_seeds[(std::size_t)k];
    } else {
      // derive seed from base seed
      
      // empirically it's noticeably better to feed the RNG sequential numbers rather than
      // try to shuffle bits around (!!)
      return base_seed + (RngSeedType)k;
    }
  }


  inline void printBasicCDataMHRWInfo(std::ostream & str) const
  {
    str << "\tstep            : " << std::setprecision(4) << mhrw_params.mhwalker_params << "\n"
        << "\t# iter. / sweep = " << mhrw_params.n_sweep << "\n"
        << "\t# therm. sweeps = " << mhrw_params.n_therm << "\n"
        << "\t# run sweeps    = " << mhrw_params.n_run << "\n";
  }

  inline std::string getBasicCDataMHRWInfo() const
  {
    std::ostringstream ss;
    printBasicCDataMHRWInfo(ss);
    return ss.str();
  }

private:
  friend boost::serialization::access;
  template<typename Archive>
  void serialize(Archive & a, unsigned int /* version */)
  {
    a & mhrw_params;
    a & base_seed;
    a & task_seeds;
  }
};







template<typename MHRWStatsResultsType_, typename IterCountIntType, typename MHWalkerParams>
struct TOMOGRAPHER_EXPORT MHRandomWalkTaskResult
  : public virtual Tools::NeedOwnOperatorNew<MHRWStatsResultsType_>::ProviderType
{
  typedef MHRWStatsResultsType_ MHRWStatsResultsType;
    
  typedef MHRWParams<MHWalkerParams, IterCountIntType> MHRWParamsType;
  
  template<typename MHRWStatsResultsTypeInit,
           typename MHRWParamsTypeInit>
  MHRandomWalkTaskResult(MHRWStatsResultsTypeInit && stats_results_,
                         MHRWParamsTypeInit && mhrw_params_,
                         double acceptance_ratio_)
    : stats_results(std::forward<MHRWStatsResultsTypeInit>(stats_results_)),
      mhrw_params(std::forward<MHRWParamsTypeInit>(mhrw_params_)),
      acceptance_ratio(acceptance_ratio_)
  {
  }

  template<typename MHRWStatsResultsTypeInit, typename MHRandomWalkType>
  MHRandomWalkTaskResult(MHRWStatsResultsTypeInit && stats_results_,
                         const MHRandomWalkType & mhrandomwalk)
    : stats_results(std::forward<MHRWStatsResultsTypeInit>(stats_results_)),
      mhrw_params(mhrandomwalk.mhrwParams()),
      acceptance_ratio(mhrandomwalk.hasAcceptanceRatio() ?
                       mhrandomwalk.acceptanceRatio() :
                       std::numeric_limits<double>::quiet_NaN())
  {
  }


  TOMOGRAPHER_ENABLED_IF(std::is_default_constructible<MHRWStatsResultsType>::value)
  MHRandomWalkTaskResult() : stats_results(), mhrw_params(), acceptance_ratio() { }

    
  MHRWStatsResultsType stats_results;
    
  MHRWParamsType mhrw_params;

  double acceptance_ratio;


  MHRandomWalkTaskResult(MHRandomWalkTaskResult && ) = default;
  MHRandomWalkTaskResult(const MHRandomWalkTaskResult & ) = default;

private:
  friend boost::serialization::access;
  template<typename Archive,
           typename MHRWStatsResultsType2 = MHRWStatsResultsType>
  void serialize(Archive & a, unsigned int /* version */)
  {
    MHRWStatsResultsType2 & stats_results_ref = stats_results;
    a & stats_results_ref;
    a & mhrw_params;
    a & acceptance_ratio;
  }
};




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

  typedef MHRWParams<MHWalkerParams, IterCountIntType> MHRWParamsType;

  typedef MHRandomWalkTaskResult<typename MHRandomWalkTaskCData::MHRWStatsResultsType,
                                 IterCountIntType, MHWalkerParams> ResultType;

  typedef MHRWStatusReport<MHRWParamsType> StatusReportType;

private:
  typename Rng::result_type _seed;

  ResultType * result;


public:

  template<typename LoggerType>
  MHRandomWalkTask(typename Rng::result_type inputseed, const MHRandomWalkTaskCData * /*pcdata*/,
                   LoggerType & logger)
    : _seed(inputseed), result(NULL)
  {
    logger.longdebug("MHRandomWalkTask", [&](std::ostream & stream) {
        stream << "inputseed=" << inputseed;
      });
  }

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

private:

  template<typename LoggerType, typename TaskMgrIface>
  struct RunInnerCallable
  {
    const MHRandomWalkTaskCData * pcdata;
    Rng & rng;
    LoggerType & baselogger;
    TaskMgrIface * tmgriface;

    ResultType ** ppresult;

    RunInnerCallable(const MHRandomWalkTaskCData * pcdata_, Rng & rng_, LoggerType & logger_,
                     TaskMgrIface * tmgriface_, ResultType ** ppresult_)
      : pcdata(pcdata_), rng(rng_), baselogger(logger_), tmgriface(tmgriface_), ppresult(ppresult_)
    {
    }
    ~RunInnerCallable() {
    }

    template<typename MHWalkerType, typename MHRWStatsCollectorType, typename MHRWControllerType>
    inline void operator()(MHWalkerType & mhwalker, MHRWStatsCollectorType & stats,
                           MHRWControllerType & controller)
    {
      // here we actually run the stuff

      auto logger = Logger::makeLocalLogger("Tomographer::MHRWTasks::MHRandomWalkTask::run()/run", baselogger);

      // our own "stats collector" which checks if we need to send a status report back
      typedef PredStatusReportMHRWStatsCollector<MHRWParamsType> OurStatusReportCheck;
      OurStatusReportCheck statreportcheck(
          // predicate
          [&]() { return this->tmgriface->statusReportRequested(); },
          // send-status-function
          [&](StatusReportType report) { this->tmgriface->submitStatusReport(std::move(report)); }
          );

      typedef MultipleMHRWStatsCollectors<MHRWStatsCollectorType, OurStatusReportCheck>
        OurStatsCollectors;
      OurStatsCollectors ourstatscollectors(stats, statreportcheck);

      logger.longdebug("About to creat MHRandomWalk instance");

      typedef MHRandomWalk<Rng,MHWalkerType,OurStatsCollectors,MHRWControllerType,LoggerType,IterCountIntType>
        MHRandomWalkType;

      MHRandomWalkType rwalk(
          // MH random walk parameters
          pcdata->mhrw_params,
          // the MHWalker
          mhwalker,
          // our stats collectors
          ourstatscollectors,
          // the random walk controller
          controller,
          // a random number generator
          rng,
          // and a logger
          baselogger
          );
      
      logger.longdebug("MHRandomWalk object created, running...");

      rwalk.run();

      logger.longdebug("MHRandomWalk run finished.");

      *ppresult = new ResultType(stats.stealResult(), rwalk);
    }

    template<typename MHWalkerType, typename MHRWStatsCollectorType>
    inline void operator()(MHWalkerType & mhwalker, MHRWStatsCollectorType & stats)
    {
      MHRWNoController ctrl;
      operator()<MHWalkerType, MHRWStatsCollectorType, MHRWNoController>(mhwalker, stats, ctrl) ;
    }
  };


public:
  template<typename LoggerType, typename TaskManagerIface>
  inline void run(const MHRandomWalkTaskCData * pcdata, LoggerType & baselogger, TaskManagerIface * tmgriface)
  {
    auto logger = Logger::makeLocalLogger("Tomographer::MHRWTasks::MHRandomWalkTask::run()", baselogger) ;

    logger.debug([&](std::ostream & stream) {
        stream << "Starting random walk, using rng seed = " << _seed;
      });
    
    Rng rng(_seed); // seeded random number generator

    RunInnerCallable<LoggerType, TaskManagerIface> run_object(pcdata, rng, baselogger, tmgriface, &result);

    pcdata->setupRandomWalkAndRun(rng, baselogger, run_object);
  }

  inline ResultType getResult() const
  {
    return *result;
  }

  inline ResultType stealResult()
  {
    return std::move(*result);
  }
};



} // namespace MHRWTasks



} // namespace Tomographer



#endif