mhrw.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_MHRW_H
#define _TOMOGRAPHER_MHRW_H

#include <cstddef>

#include <limits>
#include <random>

#include <tomographer2/tools/loggers.h>
#include <tomographer2/tools/fmt.h>
#include <tomographer2/tools/cxxutil.h>
#include <tomographer2/mhrw_bin_err.h>
#include <tomographer2/tools/needownoperatornew.h>



namespace Tomographer {








enum {
  MHUseFnValue = 1,
  MHUseFnLogValue,
  MHUseFnRelativeValue
};


namespace tomo_internal {

  template<typename MHWalker, int UseFnSyntaxType>
  struct MHRandomWalk_helper_decide_jump
  {
    typedef typename MHWalker::PointType PointType;
    typedef typename MHWalker::FnValueType FnValueType;

    static inline FnValueType get_ptval(MHWalker & mhwalker, const PointType & curpt)
    {
      (void)mhwalker; (void)curpt;
      tomographer_assert(0 && "UNKNOWN UseFnSyntaxType: Not implemented");
    }
    static inline double get_a_value(MHWalker & /*mhwalker*/, const PointType & /*newpt*/, FnValueType /*newptval*/,
                                     const PointType & /*curpt*/, FnValueType /*curptval*/)
    {
      tomographer_assert(0 && "UNKNOWN UseFnSyntaxType: Not implemented");
      return 0;
    }
  };

  template<typename MHWalker>
  struct MHRandomWalk_helper_decide_jump<MHWalker, MHUseFnValue>
  {
    typedef typename MHWalker::PointType PointType;
    typedef typename MHWalker::FnValueType FnValueType;

    static inline FnValueType get_ptval(MHWalker & mhwalker, const PointType & curpt)
    {
      return mhwalker.fnVal(curpt);
    }
    static inline double get_a_value(MHWalker & /*mhwalker*/, const PointType & /*newpt*/, double newptval,
                                     const PointType & /*curpt*/, double curptval)
    {
      return (newptval / curptval);
    }
  };

  template<typename MHWalker>
  struct MHRandomWalk_helper_decide_jump<MHWalker, MHUseFnLogValue>
  {
    typedef typename MHWalker::PointType PointType;
    typedef typename MHWalker::FnValueType FnValueType;

    static inline FnValueType get_ptval(MHWalker & mhwalker, const PointType & curpt)
    {
      return mhwalker.fnLogVal(curpt);
    }
    static inline double get_a_value(MHWalker & /*mhwalker*/, const PointType & /*newpt*/, FnValueType newptval,
                                     const PointType & /*curpt*/, FnValueType curptval)
    {
      return (newptval > curptval) ? 1.0 : exp(newptval - curptval);
    }
  };

  template<typename MHWalker>
  struct MHRandomWalk_helper_decide_jump<MHWalker, MHUseFnRelativeValue>
  {
    typedef typename MHWalker::PointType PointType;
    typedef int FnValueType; // dummy FnValueType

    static inline FnValueType get_ptval(MHWalker & /*mhwalker*/, const PointType & /*curpt*/)
    {
      return 0;
    }
    static inline double get_a_value(MHWalker & mhwalker, const PointType & newpt, FnValueType /*newptval*/,
                                     const PointType & curpt, FnValueType /*curptval*/)
    {
      return mhwalker.fnRelVal(newpt, curpt);
    }
  };
};



template<typename CountIntType_, typename StepRealType_>
struct MHRWParams
{
  typedef CountIntType_ CountIntType;
  typedef StepRealType_ StepRealType;

  explicit MHRWParams()
    : step_size(0), n_sweep(0), n_therm(0), n_run(0)
  {
  }
  MHRWParams(StepRealType step_size_, CountIntType n_sweep_, CountIntType n_therm_, CountIntType n_run_)
    : step_size(step_size_), n_sweep(n_sweep_), n_therm(n_therm_), n_run(n_run_)
  {
  }

  StepRealType step_size;

  CountIntType n_sweep;

  CountIntType n_therm;

  CountIntType n_run;
};



template<typename CountIntType, typename StepRealType>
std::ostream & operator<<(std::ostream & str, const MHRWParams<CountIntType,StepRealType> & p)
{
  str << "MHRWParams(step_size=" << p.step_size << ",n_sweep=" << p.n_sweep
      << ",n_therm=" << p.n_therm << ",n_run=" << p.n_run << ")";
  return str;
}



template<typename Rng_, typename MHWalker_, typename MHRWStatsCollector_, typename LoggerType_,
         typename CountIntType_ = int>
class MHRandomWalk
  : public virtual Tools::NeedOwnOperatorNew<typename MHWalker_::PointType>::ProviderType
{
public:
  typedef Rng_ Rng;
  typedef MHWalker_ MHWalker;
  typedef MHRWStatsCollector_ MHRWStatsCollector;
  typedef LoggerType_ LoggerType;
  typedef CountIntType_ CountIntType;

  typedef typename MHWalker::PointType PointType;
  typedef typename MHWalker::StepRealType StepRealType;

  typedef MHRWParams<CountIntType, StepRealType> MHRWParamsType;

#ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN
  typedef typename tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,MHWalker::UseFnSyntaxType>::FnValueType
      FnValueType;
#else
  typedef _FnValueType FnValueType;
#endif

  enum {
    UseFnSyntaxType = MHWalker::UseFnSyntaxType
  };

private:
  const MHRWParamsType _n;

  Rng & _rng;
  MHWalker & _mhwalker;
  MHRWStatsCollector & _stats;
  Logger::LocalLogger<LoggerType> _logger;

  PointType curpt;
  FnValueType curptval;

  CountIntType num_accepted;
  CountIntType num_live_points;

public:

  MHRandomWalk(StepRealType step_size, CountIntType n_sweep, CountIntType n_therm, CountIntType n_run,
               MHWalker & mhwalker, MHRWStatsCollector & stats,
               Rng & rng, LoggerType & logger_)
    : _n(step_size, n_sweep, n_therm, n_run),
      _rng(rng),
      _mhwalker(mhwalker),
      _stats(stats),
      _logger(TOMO_ORIGIN, logger_),
      curpt(),
      curptval(),
      num_accepted(0),
      num_live_points(0)
  {
    _logger.debug([&](std::ostream & stream) {
        stream << "constructor(). n_sweep=" << n_sweep << ", step_size=" << step_size
               << "n_therm=" << n_therm << ", n_run=" << n_run;
      });
  }
  template<typename MHRWParamsType>
  MHRandomWalk(MHRWParamsType&& n_rw,
               MHWalker & mhwalker, MHRWStatsCollector & stats,
               Rng & rng, LoggerType & logger_)
    : _n(std::forward<MHRWParamsType>(n_rw)),
      _rng(rng),
      _mhwalker(mhwalker),
      _stats(stats),
      _logger(TOMO_ORIGIN, logger_),
      curpt(),
      curptval(),
      num_accepted(0),
      num_live_points(0)
  {
    _logger.debug([&](std::ostream & s) { s << "constructor(). mhrw parameters = " << _n; });
  }

  MHRandomWalk(const MHRandomWalk & other) = delete;


  inline MHRWParamsType mhrwParams() const { return _n; }

  inline StepRealType stepSize() const { return _n.step_size; }

  inline CountIntType nSweep() const { return _n.n_sweep; }
  inline CountIntType nTherm() const { return _n.n_therm; }
  inline CountIntType nRun() const { return _n.n_run; }



  inline bool hasAcceptanceRatio() const
  {
    return (num_live_points > 0);
  }
  template<typename RatioType = double>
  inline RatioType acceptanceRatio() const
  {
    return RatioType(num_accepted) / RatioType(num_live_points);
  }


  inline const PointType & getCurrentPoint() const
  {
    return curpt;
  }

  inline const FnValueType & getCurrentPointValue() const
  {
    return curptval;
  }

  inline void setCurrentPoint(const PointType& pt)
  {
    curpt = pt;
    curptval = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_ptval(_mhwalker, curpt);
    _logger.longdebug([&](std::ostream & s) {
        s << "setCurrentPoint: set internal state. Value = " << curptval << "; Point =\n" << pt << "\n";
      });
  }


private:

  inline void _init()
  {
    num_accepted = 0;
    num_live_points = 0;

    // starting point
    curpt = _mhwalker.startPoint();
    curptval = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_ptval(_mhwalker, curpt);

    _mhwalker.init();
    _stats.init();
    _logger.longdebug("_init() done.");
  }
  inline void _thermalizing_done()
  {
    _mhwalker.thermalizingDone();
    _stats.thermalizingDone();
    _logger.longdebug("_thermalizing_done() done.");
  }
  inline void _done()
  {
    _mhwalker.done();
    _stats.done();
    _logger.longdebug("_done() done.");
  }

  inline void _move(CountIntType k, bool is_thermalizing, bool is_live_iter)
  {
    _logger.longdebug("_move()");
    // The reason `step_size` is passed to jump_fn instead of leaving jump_fn itself
    // handle the step size, is that we might in the future want to dynamically adapt the
    // step size according to the acceptance ratio. That would have to be done in this
    // class.
    PointType newpt = _mhwalker.jumpFn(curpt, _n.step_size);

    FnValueType newptval;

    newptval = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_ptval(_mhwalker, newpt);

    double a = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_a_value(
        _mhwalker, newpt, newptval, curpt, curptval
        );

    // accept move?
    bool accept = true;
    if (a < 1.0) {
      accept = bool( _rng()-_rng.min() <= a*(_rng.max()-_rng.min()) );
    }

    // track acceptance ratio, except if we are thermalizing
    if (!is_thermalizing) {
      num_accepted += accept ? 1 : 0;
      ++num_live_points;
    }

    _stats.rawMove(k, is_thermalizing, is_live_iter, accept, a, newpt, newptval, curpt, curptval, *this);

    _logger.longdebug([&](std::ostream & stream) {
        stream << (is_thermalizing?"T":"#") << std::setw(3) << k << ": " << (accept?"AC":"RJ") << " "
               << std::setprecision(4)
               << "a=" << std::setw(5) << a << ", newptval=" << std::setw(5) << newptval
               << ", curptval=" << std::setw(5) << curptval << ", accept_ratio="
               << (!is_thermalizing ? Tools::fmts("%.2g", this->acceptanceRatio()) : std::string("N/A"))
               << Tools::streamIfPossible(curpt, "\ncurpt = ", "", "");
      });

    if (accept) {
      // update the internal state of the random walk
      curpt = newpt;
      curptval = newptval;
    }
    _logger.longdebug("_move() done.");
  }

  inline void _process_sample(CountIntType k, CountIntType n)
  {
    _stats.processSample(k, n, curpt, curptval, *this);
    _logger.longdebug("_process_sample() done.");
  }


public:

  void run()
  {
    _init();

    CountIntType k;

    _logger.longdebug([&](std::ostream & s) {
        s << "Starting random walk, sweep size = " << _n.n_sweep << ", step size = " << _n.step_size
          << ", # therm sweeps = " << _n.n_therm << ", # live sweeps = " << _n.n_run;
      });

    const CountIntType num_thermalize = _n.n_sweep * _n.n_therm;

    for (k = 0; k < num_thermalize; ++k) {
      // calculate a candidate jump point and see if we accept the move
      _move(k, true, false);
    }

    _thermalizing_done();

    _logger.longdebug("Thermalizing done, starting live runs.");

    const CountIntType num_run = _n.n_sweep * _n.n_run;

    CountIntType n = 0; // number of live samples

    for (k = 0; k < num_run; ++k) {

      bool is_live_iter = ((k+1) % _n.n_sweep == 0);
      
      // calculate a candidate jump point and see if we accept the move
      _move(k, false, is_live_iter);

      if (is_live_iter) {
        _process_sample(k, n);
        ++n;
      }

    }

    _done();

    _logger.longdebug("Random walk completed.");

    return;
  }
};



} // namespace Tomographer




#endif