mhrwstatscollectors.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_MHRWSTATSCOLLECTORS_H
#define TOMOGRAPHER_MHRWSTATSCOLLECTORS_H

#include <cstddef>

#include <limits>
#include <tuple>
#include <utility>
#include <type_traits>
#include <typeinfo>


#include <tomographer/tools/cxxutil.h>
#include <tomographer/tools/loggers.h>
#include <tomographer/histogram.h>
#include <tomographer/mhrw_bin_err.h>

namespace Tomographer {


template<typename... MHRWStatsCollectors>
class MultipleMHRWStatsCollectors
{
public:
  typedef std::tuple<MHRWStatsCollectors&...> MHRWStatsCollectorsRefTupleType;
  typedef std::tuple<MHRWStatsCollectors...> MHRWStatsCollectorsTupleType;

  static constexpr int NumStatColl = sizeof...(MHRWStatsCollectors);

private:
  MHRWStatsCollectorsRefTupleType statscollectors;

public:

  MultipleMHRWStatsCollectors(MHRWStatsCollectors&... statscollectors_)
    : statscollectors(statscollectors_...)
  {
  }

  // method to get a particular stats collector
  template<int I>
  inline const typename std::tuple_element<I, MHRWStatsCollectorsTupleType>::type & getStatsCollector() const
  {
    return std::get<I>(statscollectors);
  }

  // init() callback

  template<int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type init()
  {
    std::get<I>(statscollectors).init();
    init<I+1>();
  }
  template<int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type init()
  {
  }

  // thermalizingDone() callback

  template<int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type thermalizingDone()
  {
    std::get<I>(statscollectors).thermalizingDone();
    thermalizingDone<I+1>();
  }
  template<int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type thermalizingDone()
  {
  }

  // done() callback

  template<int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type done()
  {
    std::get<I>(statscollectors).done();
    done<I+1>();
  }
  template<int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type done()
  {
  }


  // rawMove() callback

  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type rawMove(
      CountIntType k, bool is_thermalizing, bool is_live_iter, bool accepted,
      double a, const PointType & newpt, FnValueType newptval,
      const PointType & curpt, FnValueType curptval,
      MHRandomWalk & rw
      )
  {
    std::get<I>(statscollectors).rawMove(
        k, is_thermalizing, is_live_iter, accepted, a,
        newpt, newptval, curpt, curptval, rw
        );
    rawMove<CountIntType, PointType, FnValueType, MHRandomWalk, I+1>(
        k, is_thermalizing, is_live_iter, accepted, a,
        newpt, newptval, curpt, curptval, rw
        );
  }
  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type rawMove(
      CountIntType, bool, bool, bool, double, const PointType &, FnValueType,
      const PointType &, FnValueType, MHRandomWalk &
      )
  {
  }


  // processSample() callback

  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type processSample(
      CountIntType k, CountIntType n, const PointType & curpt, FnValueType curptval, MHRandomWalk & rw
      )
  {
    std::get<I>(statscollectors).processSample(k, n, curpt, curptval, rw);
    processSample<CountIntType, PointType, FnValueType, MHRandomWalk, I+1>(k, n, curpt, curptval, rw);
  }

  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type processSample(
      CountIntType, CountIntType, const PointType &, FnValueType, MHRandomWalk &
      )
  {
  }

};





typedef MultipleMHRWStatsCollectors<> TrivialMHRWStatsCollector;




// -----------------



template<typename ValueCalculator_,
         typename LoggerType = Logger::VacuumLogger,
         typename HistogramType_ = UniformBinsHistogram<typename ValueCalculator_::ValueType>
         >
class ValueHistogramMHRWStatsCollector
  : public virtual Tools::NeedOwnOperatorNew<ValueCalculator_, HistogramType_>::ProviderType,
    public virtual Tools::EigenAlignedOperatorNewProvider // DEBUG:: REDUNDANT???
{
public:
  typedef ValueCalculator_ ValueCalculator;

  typedef typename ValueCalculator::ValueType ValueType;

  typedef HistogramType_ HistogramType;

  typedef HistogramType_ ResultType;

  typedef typename HistogramType::Params HistogramParams;

private:

  HistogramType _histogram;

  ValueCalculator _vcalc;

  LoggerType & _logger;

public:
  ValueHistogramMHRWStatsCollector(HistogramParams histogram_params,
                                   const ValueCalculator & vcalc,
                                   LoggerType & logger)
    : _histogram(histogram_params),
      _vcalc(vcalc),
      _logger(logger)
  {
  }

  inline const HistogramType & histogram() const
  {
    return _histogram;
  }

  inline const ResultType & getResult() const
  {
    return _histogram;
  }

  // stats collector part

  inline void init()
  {
    // reset our array
    _histogram.reset();
  }
  inline void thermalizingDone()
  {
  }
  template<bool PrintHistogram = true>
  inline void done()
  {
    if (PrintHistogram) {
      if (_logger.enabledFor(Logger::LONGDEBUG)) {
        // _logger.longdebug("ValueHistogramMHRWStatsCollector", "done()");
        _logger.longdebug("ValueHistogramMHRWStatsCollector",
                       "Done walking & collecting stats. Here's the histogram:\n"
                       + _histogram.prettyPrint());
      }
    }
  }

  template<typename CountIntType, typename PointType, typename LLHValueType, typename MHRandomWalk>
  void rawMove(CountIntType k, bool /*is_thermalizing*/, bool /*is_live_iter*/, bool /*accepted*/,
                double /*a*/, const PointType & /*newpt*/, LLHValueType /*newptval*/,
                const PointType & /*curpt*/, LLHValueType /*curptval*/, MHRandomWalk & /*mh*/)
  {
    _logger.longdebug("ValueHistogramMHRWStatsCollector", [&](std::ostream & stream) {
        stream << "rawMove(): k=" << k;
      });
  }

  template<typename CountIntType, typename PointType, typename LLHValueType, typename MHRandomWalk>
  std::size_t processSample(CountIntType k, CountIntType n, const PointType & curpt,
                             LLHValueType /*curptval*/, MHRandomWalk & /*mh*/)
  {
    ValueType val = _vcalc.getValue(curpt);

    _logger.longdebug("ValueHistogramMHRWStatsCollector", [&](std::ostream & stream) {
        stream << "in processSample(): "
               << "k=" << k << ", n=" << n << ", val=" << val
               << " [with ValueType=" << typeid(ValueType).name() << "]" ;
      });

    return _histogram.record(val);

    //_logger.longdebug("ValueHistogramMHRWStatsCollector", "processSample() finished");
  }
 

};




template<typename HistogramType_, typename BinningAnalysisParamsType_>
struct ValueHistogramWithBinningMHRWStatsCollectorResult
  : public virtual Tools::NeedOwnOperatorNew<
  HistogramType_,
  typename BinningAnalysisParamsType_::BinSumSqArray
  >::ProviderType
{
  typedef HistogramType_ HistogramType;
  typedef typename HistogramType::Params HistogramParams;
  typedef BinningAnalysisParamsType_ BinningAnalysisParamsType;

  explicit ValueHistogramWithBinningMHRWStatsCollectorResult()
    : hist(), error_levels(), converged_status()
  {
  }

  template<typename EigenDerived1, typename EigenDerived2>
  ValueHistogramWithBinningMHRWStatsCollectorResult(const HistogramType & hist_,
                                                    const Eigen::DenseBase<EigenDerived1> & error_levels_,
                                                    const Eigen::DenseBase<EigenDerived2> & converged_status_)
    : hist(hist_), error_levels(error_levels_), converged_status(converged_status_)
  {
  }

  template<typename BinningAnalysisType>
  ValueHistogramWithBinningMHRWStatsCollectorResult(HistogramParams p, const BinningAnalysisType & b)
    : hist(p),
      error_levels(b.numTrackValues(), b.numLevels()+1),
      converged_status(Eigen::ArrayXi::Constant(b.numTrackValues(), BinningAnalysisType::UNKNOWN_CONVERGENCE))
  {
    tomographer_assert(converged_status.rows() == b.numTrackValues() && converged_status.cols() == 1);
  }

  HistogramType hist;

  typename BinningAnalysisParamsType::BinSumSqArray error_levels;

  Eigen::ArrayXi converged_status;
  
  inline void dumpConvergenceAnalysis(std::ostream & str) const
  {
    for (int k = 0; k < converged_status.size(); ++k) {
      str << "\tval[" << std::setw(3) << k << "] = "
          << std::setw(12) << hist.bins(k)
          << " +- " << std::setw(12) << hist.delta(k);
      if (converged_status(k) == BinningAnalysisParamsType::CONVERGED) {
          str << "  [CONVERGED]";
      } else if (converged_status(k) == BinningAnalysisParamsType::NOT_CONVERGED) {
        str << "  [NOT CONVERGED]";
      } else if (converged_status(k) == BinningAnalysisParamsType::UNKNOWN_CONVERGENCE) {
        str << "  [UNKNOWN]";
      } else {
        str << "  [UNKNOWN CONVERGENCE STATUS: " << converged_status(k) << "]";
      }
      str << "\n";
    }
  }

  inline std::string dumpConvergenceAnalysis() const
  {
    std::stringstream ss;
    dumpConvergenceAnalysis(ss);
    return ss.str();
  }
  
};




template<typename ValueCalculator_,
         typename CountIntType_ = int,
         typename CountRealAvgType_ = double,
         int NumTrackValues_ = Eigen::Dynamic,
         int NumLevels_ = Eigen::Dynamic
         >
struct ValueHistogramWithBinningMHRWStatsCollectorParams
{
  typedef ValueCalculator_ ValueCalculator;
  typedef CountIntType_ CountIntType;
  typedef CountRealAvgType_ CountRealAvgType;

  static constexpr int NumTrackValues = NumTrackValues_;
  static constexpr int NumLevels = NumLevels_;

  typedef typename ValueCalculator::ValueType ValueType;
  
  typedef BinningAnalysisParams<ValueType,NumTrackValues,NumLevels,false/*StoreBinSums*/,CountIntType>
  BinningAnalysisParamsType;

  typedef UniformBinsHistogram<typename ValueCalculator::ValueType, CountIntType> BaseHistogramType;
  typedef UniformBinsHistogramWithErrorBars<typename ValueCalculator::ValueType, CountRealAvgType> HistogramType;
  typedef typename HistogramType::Params HistogramParams;

  typedef ValueHistogramWithBinningMHRWStatsCollectorResult<HistogramType,BinningAnalysisParamsType>  Result;

};

template<typename Params,
         typename LoggerType_ = Logger::VacuumLogger
         >
class ValueHistogramWithBinningMHRWStatsCollector
  : public virtual Tools::NeedOwnOperatorNew<ValueHistogramMHRWStatsCollector<
                                               typename Params::ValueCalculator,
                                               LoggerType_,
                                               typename Params::BaseHistogramType
                                               >,
                                             BinningAnalysis<typename Params::BinningAnalysisParamsType, LoggerType_>,
                                             typename Params::Result
                                             >::ProviderType
{
public:

  typedef typename Params::ValueCalculator ValueCalculator;
  typedef typename Params::CountIntType CountIntType;
  typedef typename Params::CountRealAvgType CountRealAvgType;
  
  typedef LoggerType_ LoggerType;

  typedef typename Params::BaseHistogramType BaseHistogramType;
  typedef typename Params::HistogramParams HistogramParams;
    
  typedef typename Params::ValueType ValueType;
  
  typedef typename Params::BinningAnalysisParamsType BinningAnalysisParamsType;
  typedef BinningAnalysis<BinningAnalysisParamsType, LoggerType> BinningAnalysisType;
    
  static constexpr int NumTrackValuesCTime = Params::NumTrackValues;
  static constexpr int NumLevelsCTime = Params::NumLevels;

  typedef typename Params::Result ResultType;

  typedef ValueHistogramMHRWStatsCollector<
    ValueCalculator,
    LoggerType,
    BaseHistogramType
    > ValueHistogramMHRWStatsCollectorType;
  
private:
    
  ValueHistogramMHRWStatsCollectorType value_histogram;

  BinningAnalysisType binning_analysis;

  LoggerType & logger;

  ResultType result;

public:
    
  ValueHistogramWithBinningMHRWStatsCollector(HistogramParams histogram_params,
                                              const ValueCalculator & vcalc,
                                              int num_levels,
                                              LoggerType & logger_)
    : value_histogram(histogram_params, vcalc, logger_),
      binning_analysis(histogram_params.num_bins, num_levels, logger_),
      logger(logger_),
      result(histogram_params, binning_analysis)
  {
    logger.longdebug("ValueHistogramWithBinningMHRWStatsCollector", "constructor()");
  }
  
  inline const BaseHistogramType & histogram() const
  {
    return value_histogram.histogram();
  }

  inline const BinningAnalysisType & getBinningAnalysis() const
  {
    return binning_analysis;
  }

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

  // stats collector part

  inline void init()
  {
    value_histogram.init();
  }
  inline void thermalizingDone()
  {
    value_histogram.thermalizingDone();
  }
  inline void done()
  {
    logger.longdebug("ValueHistogramWithBinningMHRWStatsCollector::done()", "finishing up ...");

    value_histogram.template done<false>();

    //
    // Determine the error bars from the binning analysis. Remember, the binning analysis was
    // applied to each of the indicator functions "chi(value) = (value in bin # i) ? 1 : 0"
    //
    // The total number of samples is simply h.bins.sum()+h.off_chart; this is the relevant
    // coefficient to calculate the bin means needed by binning_analysis.calcErrorLevels().  Indeed,
    // in this way we really get the averaged observed value of each indicator function for each
    // value interval.
    //
    const BaseHistogramType & h = value_histogram.histogram();
    result.hist.params = h.params;
    CountRealAvgType numsamples = h.bins.sum() + h.off_chart;
    result.hist.bins = h.bins.template cast<CountRealAvgType>() / numsamples;
    result.error_levels = binning_analysis.calcErrorLevels(result.hist.bins);
    result.hist.delta = result.error_levels.col(binning_analysis.numLevels()).template cast<CountRealAvgType>();
    result.hist.off_chart = h.off_chart / numsamples;

    result.converged_status = binning_analysis.determineErrorConvergence(result.error_levels);

    logger.debug("ValueHistogramWithBinningMHRWStatsCollector", [&,this](std::ostream & str) {
        str << "Binning analysis: bin sqmeans at different binning levels are:\n"
            << binning_analysis.getBinSqmeans() << "\n"
            << "\t-> so the error bars at different binning levels are:\n"
            << result.error_levels << "\n"
            << "\t-> convergence analysis: \n";
        result.dumpConvergenceAnalysis(str);
        str << "\t... and just for you, here is the final histogram:\n" << result.hist.prettyPrint() << "\n";
      });
  }

  template<typename CountIntType2, typename PointType, typename LLHValueType, typename MHRandomWalk>
  inline void rawMove(CountIntType2 k, bool is_thermalizing, bool is_live_iter, bool accepted,
                      double a, const PointType & newpt, LLHValueType newptval,
                      const PointType & curpt, LLHValueType curptval, MHRandomWalk & mh)
  {
    value_histogram.rawMove(k, is_thermalizing, is_live_iter, accepted, a, newpt, newptval, curpt, curptval, mh);
  }

  template<typename CountIntType2, typename PointType, typename LLHValueType, typename MHRandomWalk>
  inline void processSample(CountIntType2 k, CountIntType2 n, const PointType & curpt,
                            LLHValueType curptval, MHRandomWalk & mh)
  {
    std::size_t histindex = value_histogram.processSample(k, n, curpt, curptval, mh);
    binning_analysis.processNewValues(
        Tools::canonicalBasisVec<Eigen::Array<ValueType,Eigen::Dynamic,1> >(
            histindex,
            value_histogram.histogram().numBins()
            )
        );
  }

};






template<typename MHRWStatsCollector_>
struct MHRWStatsCollectorStatus
{
  typedef MHRWStatsCollector_ MHRWStatsCollector;

  static constexpr bool CanProvideStatus = false;

  static inline std::string getStatus(const MHRWStatsCollector * /*stats*/)
  {
    return std::string();
  }
};
// static members:
template<typename MHRWStatsCollector_>
constexpr bool MHRWStatsCollectorStatus<MHRWStatsCollector_>::CanProvideStatus;



template<typename... Args>
struct MHRWStatsCollectorStatus<MultipleMHRWStatsCollectors<Args... > >
{
  typedef MultipleMHRWStatsCollectors<Args... > MHRWStatsCollector;

  static constexpr int NumStatColl = MHRWStatsCollector::NumStatColl;
  
  static constexpr bool CanProvideStatus = true;

  template<int I = 0, typename std::enable_if<(I < NumStatColl), bool>::type dummy = true>
  static inline std::string getStatus(const MHRWStatsCollector * stats)
  {
    typedef typename std::tuple_element<I, typename MHRWStatsCollector::MHRWStatsCollectorsTupleType>::type
      ThisStatsCollector;
    return
      (MHRWStatsCollectorStatus<ThisStatsCollector>::CanProvideStatus
       ? (MHRWStatsCollectorStatus<ThisStatsCollector>::getStatus(& stats->template getStatsCollector<I>())
          + ((I < (NumStatColl-1)) ? std::string("\n") : std::string()))
       : std::string())
      + getStatus<I+1>(stats);
  };

  template<int I = 0, typename std::enable_if<(I == NumStatColl), bool>::type dummy = true>
  static inline std::string getStatus(const MHRWStatsCollector * stats)
  {
    (void)stats;
    return std::string();
  }

};
// static members:
template<typename... Args>
constexpr int MHRWStatsCollectorStatus<MultipleMHRWStatsCollectors<Args... > >::NumStatColl;
template<typename... Args>
constexpr bool MHRWStatsCollectorStatus<MultipleMHRWStatsCollectors<Args... > >::CanProvideStatus;



template<typename ValueCalculator_,
         typename LoggerType_,
         typename HistogramType_
         >
struct MHRWStatsCollectorStatus<ValueHistogramMHRWStatsCollector<ValueCalculator_, LoggerType_, HistogramType_> >
{
  typedef ValueHistogramMHRWStatsCollector<ValueCalculator_, LoggerType_, HistogramType_> MHRWStatsCollector;
  
  static constexpr bool CanProvideStatus = true;

  static inline std::string getStatus(const MHRWStatsCollector * stats)
  {
    const int maxbarwidth = 50;

    typedef typename MHRWStatsCollector::HistogramType HistogramType;

    return "Histogram: " + histogramShortBar<HistogramType>(stats->histogram(), true, maxbarwidth);
  }
};
// static members:
template<typename ValueCalculator_,
         typename LoggerType_,
         typename HistogramType_
         >
constexpr bool
MHRWStatsCollectorStatus<ValueHistogramMHRWStatsCollector<ValueCalculator_, LoggerType_, HistogramType_> >::CanProvideStatus;


template<typename Params_,
         typename LoggerType_
         >
struct MHRWStatsCollectorStatus<ValueHistogramWithBinningMHRWStatsCollector<Params_, LoggerType_> >
{
  typedef ValueHistogramWithBinningMHRWStatsCollector<Params_, LoggerType_> MHRWStatsCollector;
  
  static constexpr bool CanProvideStatus = true;

  static inline std::string getStatus(const MHRWStatsCollector * stats)
  {
    const int maxbarwidth = 50;

    typedef typename MHRWStatsCollector::BaseHistogramType BaseHistogramType;
    const BaseHistogramType & histogram = stats->histogram();

    // calculate the error bars at different levels, to determine convergence status.
    typedef typename MHRWStatsCollector::BinningAnalysisType BinningAnalysisType;
    //typedef typename MHRWStatsCollector::CountRealAvgType CountRealAvgType;
    typedef typename BinningAnalysisType::ValueType ValueType;
    const auto& binning_analysis = stats->getBinningAnalysis();
    Eigen::Array<ValueType, Eigen::Dynamic, 1> binmeans(histogram.numBins());
    binmeans = histogram.bins.template cast<ValueType>() /
      (ValueType)(histogram.bins.sum() + histogram.off_chart);

    auto error_levels = binning_analysis.calcErrorLevels(binmeans);
    auto conv_status = binning_analysis.determineErrorConvergence(error_levels);

    int n_cnvg = 0;
    int n_unknown = 0;
    int n_not_cnvg = 0;
    for (std::size_t k = 0; k < (std::size_t)histogram.numBins(); ++k) {
      if (conv_status(k) == BinningAnalysisType::CONVERGED) {
        ++n_cnvg;
      } else if (conv_status(k) == BinningAnalysisType::NOT_CONVERGED) {
        ++n_not_cnvg;
      } else {
        ++n_unknown;
      }
    }

    return tomo_internal::histogram_short_bar_fmt<BaseHistogramType>(histogram, "", maxbarwidth)
      + Tools::fmts("   err: (cnvg/?/fail) %d/%d/%d", n_cnvg, n_unknown, n_not_cnvg);
  }
};
// static members:
template<typename Params_,
         typename LoggerType_
         >
constexpr bool
MHRWStatsCollectorStatus<ValueHistogramWithBinningMHRWStatsCollector<Params_, LoggerType_> >::CanProvideStatus;


} // namespace Tomographer



#endif