Tomographer  v5.4
Tomographer C++ Framework Documentation
mhrwstatscollectors.h
Go to the documentation of this file.
1 /* This file is part of the Tomographer project, which is distributed under the
2  * terms of the MIT license.
3  *
4  * The MIT License (MIT)
5  *
6  * Copyright (c) 2016 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
7  * Copyright (c) 2017 Caltech, Institute for Quantum Information and Matter, Philippe Faist
8  *
9  * Permission is hereby granted, free of charge, to any person obtaining a copy
10  * of this software and associated documentation files (the "Software"), to deal
11  * in the Software without restriction, including without limitation the rights
12  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13  * copies of the Software, and to permit persons to whom the Software is
14  * furnished to do so, subject to the following conditions:
15  *
16  * The above copyright notice and this permission notice shall be included in
17  * all copies or substantial portions of the Software.
18  *
19  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25  * SOFTWARE.
26  */
27 
28 #ifndef TOMOGRAPHER_MHRWSTATSCOLLECTORS_H
29 #define TOMOGRAPHER_MHRWSTATSCOLLECTORS_H
30 
31 #include <cstddef>
32 
33 #include <limits>
34 #include <tuple>
35 #include <utility>
36 #include <type_traits>
37 #include <typeinfo>
38 
39 
42 #include <tomographer/tools/statusprovider.h>
43 #include <tomographer/histogram.h>
44 #include <tomographer/mhrw.h> // MHRWStatusReport
46 
57 namespace Tomographer {
58 
64 struct TOMOGRAPHER_EXPORT MHRWStatsCollectorNoResult
65 {
66  // nothing.
67 
68 private:
69  friend boost::serialization::access;
70  template<typename Archive>
71  void serialize(Archive & /* a */, unsigned int /* version */) { }
72 };
73 
74 
75 namespace tomo_internal {
76 
77 typedef MHRWStatsCollectorNoResult NoResult; // convenience alias
78 
79 template<typename MHRWStatsCollector, typename dummy = void>
80 struct multistatscoll_result_type_helper {
81  // static constexpr bool HasResultType = false;
82  typedef NoResult ResultType;
83  static inline ResultType getResult(const MHRWStatsCollector & ) { return NoResult(); }
84  static inline ResultType stealResult(MHRWStatsCollector & ) { return NoResult(); }
85 };
86 template<typename MHRWStatsCollector>
87 struct multistatscoll_result_type_helper<
88  MHRWStatsCollector,
89  // ignore specialization (SFINAE) if either ResultType is not defined or if ResultType is void
90  typename std::enable_if<
91  !std::is_same<typename MHRWStatsCollector::ResultType,void>::value,
92  void
93  >::type
94  >
95 {
96  // static constexpr bool HasResultType = true;
97  typedef typename MHRWStatsCollector::ResultType ResultType;
98  static inline ResultType getResult(const MHRWStatsCollector & s) { return s.getResult(); }
99  static inline ResultType stealResult(MHRWStatsCollector & s) { return s.stealResult(); }
100 };
101 
102 template<typename ResultType, typename... MHRWStatsCollectors>
103 struct multistatscoll_result_type_helper2 {
104  typedef std::tuple<MHRWStatsCollectors&...> TupleRefType;
105  typedef std::tuple<MHRWStatsCollectors...> TupleType;
106  template<typename T, int... Is>
107  static inline ResultType getResult(const TupleRefType & statscollectors, Tools::Sequence<T, Is...> )
108  {
109  return std::make_tuple(
110  multistatscoll_result_type_helper<typename std::tuple_element<Is, TupleType>::type>::getResult(
111  std::get<Is>(statscollectors)
112  ) ...
113  ) ;
114  }
115  template<typename T, int... Is>
116  static inline ResultType stealResult(TupleRefType & statscollectors, Tools::Sequence<T, Is...> )
117  {
118  return std::make_tuple(
119  std::move(
120  multistatscoll_result_type_helper<typename std::tuple_element<Is, TupleType>::type>::stealResult(
121  std::get<Is>(statscollectors)
122  )
123  ) ...
124  ) ;
125  }
126 };
127 
128 } // namespace tomo_internal
129 
130 
131 
166 template<typename... MHRWStatsCollectors>
167 class TOMOGRAPHER_EXPORT MultipleMHRWStatsCollectors
168 {
169 public:
170  typedef std::tuple<MHRWStatsCollectors&...> MHRWStatsCollectorsRefTupleType;
171  typedef std::tuple<MHRWStatsCollectors...> MHRWStatsCollectorsTupleType;
172 
174  static constexpr int NumStatColl = sizeof...(MHRWStatsCollectors);
175 
176 private:
177  MHRWStatsCollectorsRefTupleType statscollectors;
178 
179 public:
180 
181  MultipleMHRWStatsCollectors(MHRWStatsCollectors&... statscollectors_)
182  : statscollectors(statscollectors_...)
183  {
184  }
185 
188 
189  // method to get a particular stats collector
190  template<int I>
191  inline const typename std::tuple_element<I, MHRWStatsCollectorsTupleType>::type & getStatsCollector() const
192  {
193  return std::get<I>(statscollectors);
194  }
195 
197  ResultType;
198 
199  ResultType getResult() const
200  {
201  return tomo_internal::multistatscoll_result_type_helper2<ResultType, MHRWStatsCollectors...>::getResult(
202  statscollectors,
204  );
205  }
206  ResultType stealResult()
207  {
208  return tomo_internal::multistatscoll_result_type_helper2<ResultType, MHRWStatsCollectors...>::stealResult(
209  statscollectors,
211  );
212  }
213 
214 
215  // init() callback
216 
217  template<int I = 0>
218  inline typename std::enable_if<I < NumStatColl, void>::type init()
219  {
220  std::get<I>(statscollectors).init();
221  init<I+1>();
222  }
223  template<int I = 0>
224  inline typename std::enable_if<I == NumStatColl, void>::type init()
225  {
226  }
227 
228  // thermalizingDone() callback
229 
230  template<int I = 0>
231  inline typename std::enable_if<I < NumStatColl, void>::type thermalizingDone()
232  {
233  std::get<I>(statscollectors).thermalizingDone();
234  thermalizingDone<I+1>();
235  }
236  template<int I = 0>
237  inline typename std::enable_if<I == NumStatColl, void>::type thermalizingDone()
238  {
239  }
240 
241  // done() callback
242 
243  template<int I = 0>
244  inline typename std::enable_if<I < NumStatColl, void>::type done()
245  {
246  std::get<I>(statscollectors).done();
247  done<I+1>();
248  }
249  template<int I = 0>
250  inline typename std::enable_if<I == NumStatColl, void>::type done()
251  {
252  }
253 
254 
255  // rawMove() callback
256 
257  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
258  inline typename std::enable_if<I < NumStatColl, void>::type rawMove(
259  CountIntType k, bool is_thermalizing, bool is_live_iter, bool accepted,
260  double a, const PointType & newpt, FnValueType newptval,
261  const PointType & curpt, FnValueType curptval,
262  MHRandomWalk & rw
263  )
264  {
265  std::get<I>(statscollectors).rawMove(
266  k, is_thermalizing, is_live_iter, accepted, a,
267  newpt, newptval, curpt, curptval, rw
268  );
269  rawMove<CountIntType, PointType, FnValueType, MHRandomWalk, I+1>(
270  k, is_thermalizing, is_live_iter, accepted, a,
271  newpt, newptval, curpt, curptval, rw
272  );
273  }
274  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
275  inline typename std::enable_if<I == NumStatColl, void>::type rawMove(
276  CountIntType, bool, bool, bool, double, const PointType &, FnValueType,
277  const PointType &, FnValueType, MHRandomWalk &
278  )
279  {
280  }
281 
282 
283  // processSample() callback
284 
285  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
286  inline typename std::enable_if<I < NumStatColl, void>::type processSample(
287  CountIntType k, CountIntType n, const PointType & curpt, FnValueType curptval, MHRandomWalk & rw
288  )
289  {
290  std::get<I>(statscollectors).processSample(k, n, curpt, curptval, rw);
291  processSample<CountIntType, PointType, FnValueType, MHRandomWalk, I+1>(k, n, curpt, curptval, rw);
292  }
293 
294  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
295  inline typename std::enable_if<I == NumStatColl, void>::type processSample(
296  CountIntType, CountIntType, const PointType &, FnValueType, MHRandomWalk &
297  )
298  {
299  }
300 
301 };
302 
310 
311 
312 
318 template<typename... MHRWStatsCollectors>
319 inline MultipleMHRWStatsCollectors<MHRWStatsCollectors...>
320 mkMultipleMHRWStatsCollectors(MHRWStatsCollectors& ... mhrwstatscollectors)
321 {
322  return MultipleMHRWStatsCollectors<MHRWStatsCollectors...>(mhrwstatscollectors...) ;
323 }
324 
325 
326 
327 
328 
329 
330 // -----------------
331 
336 template<typename IterCountIntType_ = int>
338 {
339 public:
340  typedef IterCountIntType_ IterCountIntType;
341 
342 private:
343  Eigen::ArrayXi accept_buffer;
344  IterCountIntType pos;
345 
346 public:
347 
348  inline MHRWMovingAverageAcceptanceRatioStatsCollector(Eigen::Index num_samples_ = 2048)
349  : accept_buffer(Eigen::ArrayXi::Zero(num_samples_)), pos(0)
350  {
351  }
352 
353  inline Eigen::Index bufferSize() const {
354  return accept_buffer.size();
355  }
356 
358  inline double movingAverageAcceptanceRatio() const {
359  return (double)accept_buffer.count() / accept_buffer.size();
360  }
361 
363  inline bool hasMovingAverageAcceptanceRatio() const {
364  return ( pos >= (IterCountIntType)accept_buffer.size() );
365  }
366 
367 
368  // init() callback
369 
370  inline void init()
371  {
372  }
373 
374  // thermalizingDone() callback
375  inline void thermalizingDone()
376  {
377  }
378 
379  // done() callback
380  inline void done()
381  {
382  }
383 
384 
385  // rawMove() callback
386 
387  template<typename CountIntType, typename PointType, typename FnValueType, typename PointType2, typename MHRandomWalk>
388  inline void rawMove(
389  CountIntType /*k*/, bool /*is_thermalizing*/, bool /*is_live_iter*/, bool accepted,
390  double /*a*/, PointType && /*newpt*/, FnValueType /*newptval*/,
391  PointType2 && /*curpt*/, FnValueType /*curptval*/,
392  MHRandomWalk && /* rw */
393  )
394  {
395  // Strategy: update the item in the array at position pos%num_samples, and increment
396  // pos. This way we remove the oldest samples and replace them by the new ones.
397 
398  accept_buffer[ (Eigen::Index)(pos % (IterCountIntType)accept_buffer.size()) ] = (int)accepted;
399  ++pos;
400  }
401 
402 
403  // processSample() callback
404 
405  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk>
406  inline void processSample(
407  CountIntType /*k*/, CountIntType /*n*/, PointType && /*curpt*/, FnValueType /*curptval*/, MHRandomWalk && /*rw*/
408  )
409  {
410  }
411 
412 };
413 
414 
415 
416 
417 // -----------------
418 
419 
420 
435 template<typename ValueCalculator_,
436  typename LoggerType = Logger::VacuumLogger,
437  typename HistogramType_ = Histogram<typename ValueCalculator_::ValueType>
438  >
439 class TOMOGRAPHER_EXPORT ValueHistogramMHRWStatsCollector
440  : public virtual Tools::NeedOwnOperatorNew<ValueCalculator_, HistogramType_>::ProviderType
441 {
442 public:
446  typedef ValueCalculator_ ValueCalculator;
447 
449  typedef typename ValueCalculator::ValueType ValueType;
450 
452  typedef HistogramType_ HistogramType;
453 
455  typedef HistogramType_ ResultType;
456 
458  typedef typename HistogramType::Params HistogramParams;
459 
460 private:
461 
463  HistogramType _histogram;
464 
469  ValueCalculator _vcalc;
470 
471  LoggerType & _logger;
472 
473 public:
475  ValueHistogramMHRWStatsCollector(HistogramParams histogram_params,
476  const ValueCalculator & vcalc,
477  LoggerType & logger)
478  : _histogram(histogram_params),
479  _vcalc(vcalc),
480  _logger(logger)
481  {
482  }
483 
485  inline const HistogramType & histogram() const
486  {
487  return _histogram;
488  }
489 
499  inline const ResultType & getResult() const
500  {
501  return _histogram;
502  }
503 
511  inline ResultType stealResult()
512  {
513  return std::move(_histogram);
514  }
515 
516  // stats collector part
517 
519  inline void init()
520  {
521  // reset our array
522  _histogram.reset();
523  }
525  inline void thermalizingDone()
526  {
527  }
535  template<bool PrintHistogram = true>
536  inline void done()
537  {
538  if (PrintHistogram) {
539  if (_logger.enabledFor(Logger::LONGDEBUG)) {
540  // _logger.longdebug("ValueHistogramMHRWStatsCollector", "done()");
541  _logger.longdebug("ValueHistogramMHRWStatsCollector",
542  "Done walking & collecting stats. Here's the histogram:\n"
543  + _histogram.prettyPrint());
544  }
545  }
546  }
547 
549  template<typename CountIntType, typename PointType, typename LLHValueType, typename MHRandomWalk>
550  void rawMove(CountIntType k, bool /*is_thermalizing*/, bool /*is_live_iter*/, bool /*accepted*/,
551  double /*a*/, const PointType & /*newpt*/, LLHValueType /*newptval*/,
552  const PointType & /*curpt*/, LLHValueType /*curptval*/, MHRandomWalk & /*mh*/)
553  {
554  _logger.longdebug("ValueHistogramMHRWStatsCollector", [&](std::ostream & stream) {
555  stream << "rawMove(): k=" << k;
556  });
557  }
558 
560  template<typename CountIntType, typename PointType, typename LLHValueType, typename MHRandomWalk>
561  Eigen::Index processSample(CountIntType k, CountIntType n, const PointType & curpt,
562  LLHValueType /*curptval*/, MHRandomWalk & /*mh*/)
563  {
564  ValueType val = _vcalc.getValue(curpt);
565 
566  _logger.longdebug("ValueHistogramMHRWStatsCollector", [&](std::ostream & stream) {
567  stream << "in processSample(): "
568  << "k=" << k << ", n=" << n << ", val=" << val << " -> bin # "
569  << ( _histogram.params.isWithinBounds(val)
570  ? std::to_string(_histogram.params.binIndexUnsafe(val))
571  : std::string("<off-chart>") )
572  << " [with ValueType=" << typeid(ValueType).name() << "]" ;
573  });
574 
575  return _histogram.record(val);
576 
577  //_logger.longdebug("ValueHistogramMHRWStatsCollector", "processSample() finished");
578  }
579 
580 
581 };
582 
583 
584 
585 
594 template<typename HistogramType_, typename BinningAnalysisParamsType_>
596  : public virtual Tools::NeedOwnOperatorNew<
597  HistogramType_,
598  typename BinningAnalysisParamsType_::BinSumSqArray
599  >::ProviderType
600 {
601  typedef HistogramType_ HistogramType;
602  typedef typename HistogramType::Params HistogramParams;
603  typedef BinningAnalysisParamsType_ BinningAnalysisParamsType;
604 
607  : histogram(), error_levels(), converged_status()
608  {
609  }
610 
611  // can be moved
613  = default;
614 
615  // can be copied
618  )
619  : histogram(HistogramType::copy(copy.histogram)),
620  error_levels(copy.error_levels),
621  converged_status(copy.converged_status)
622  {
623  }
624 
626  template<typename EigenDerived1, typename EigenDerived2>
627  ValueHistogramWithBinningMHRWStatsCollectorResult(const HistogramType & histogram_,
628  const Eigen::DenseBase<EigenDerived1> & error_levels_,
629  const Eigen::DenseBase<EigenDerived2> & converged_status_)
630  : histogram(HistogramType::copy(histogram_)),
631  error_levels(error_levels_),
632  converged_status(converged_status_)
633  {
634  }
635 
637  template<typename BinningAnalysisType>
638  ValueHistogramWithBinningMHRWStatsCollectorResult(HistogramParams p, const BinningAnalysisType & b)
639  : histogram(p),
640  error_levels(b.numTrackValues(), b.numLevels()+1),
641  converged_status(Eigen::ArrayXi::Constant(b.numTrackValues(), BinningAnalysisType::UNKNOWN_CONVERGENCE))
642  {
643  tomographer_assert(converged_status.rows() == b.numTrackValues() && converged_status.cols() == 1);
644  }
645 
655  HistogramType histogram;
656 
659  typename BinningAnalysisParamsType::BinSumSqArray error_levels;
660 
664  Eigen::ArrayXi converged_status;
665 
668  {
670  }
671 
672 
674  inline void dumpConvergenceAnalysis(std::ostream & stream) const
675  {
676  for (int k = 0; k < converged_status.size(); ++k) {
677  stream << "\tval[" << std::setw(3) << k << "] = "
678  << std::setw(12) << histogram.bins(k)
679  << " +- " << std::setw(12) << histogram.delta(k);
680  if (converged_status(k) == BINNING_CONVERGED) {
681  stream << " [CONVERGED]";
682  } else if (converged_status(k) == BINNING_NOT_CONVERGED) {
683  stream << " [NOT CONVERGED]";
684  } else if (converged_status(k) == BINNING_UNKNOWN_CONVERGENCE) {
685  stream << " [UNKNOWN]";
686  } else {
687  stream << " [INVALID CONVERGENCE STATUS: " << converged_status(k) << "]";
688  }
689  stream << "\n";
690  }
691  }
692 
695  {
697  dumpConvergenceAnalysis(ss);
698  return ss.str();
699  }
700 
701 private:
702  friend boost::serialization::access;
703  template<typename Archive>
704  void serialize(Archive & a, unsigned int /* version */)
705  {
706  a & histogram;
707  a & error_levels;
708  a & converged_status;
709  }
710 };
711 
712 
713 
714 
722 template<typename ValueCalculator_,
723  typename CountIntType_ = int,
724  typename CountRealAvgType_ = double,
725  int NumTrackValues_ = Eigen::Dynamic,
726  int NumLevels_ = Eigen::Dynamic
727  >
729 {
732  typedef ValueCalculator_ ValueCalculator;
734  typedef CountIntType_ CountIntType;
736  typedef CountRealAvgType_ CountRealAvgType;
737 
739  static constexpr int NumTrackValues = NumTrackValues_;
741  static constexpr int NumLevels = NumLevels_;
742 
744  typedef typename ValueCalculator::ValueType ValueType;
745 
747  typedef BinningAnalysisParams<ValueType,NumTrackValues,NumLevels,false/*StoreBinSums*/,CountIntType>
749 
763 
767 
768 };
769 
777 template<typename Params,
778  typename LoggerType_ = Logger::VacuumLogger
779  >
781  : public virtual Tools::NeedOwnOperatorNew<ValueHistogramMHRWStatsCollector<
782  typename Params::ValueCalculator,
783  LoggerType_,
784  typename Params::BaseHistogramType
785  >,
786  BinningAnalysis<typename Params::BinningAnalysisParamsType, LoggerType_>,
787  typename Params::Result
788  >::ProviderType
789 {
790 public:
791 
793  typedef typename Params::ValueCalculator ValueCalculator;
795  typedef typename Params::CountIntType CountIntType;
797  typedef typename Params::CountRealAvgType CountRealAvgType;
798 
800  typedef LoggerType_ LoggerType;
801 
803  typedef typename Params::BaseHistogramType BaseHistogramType;
805  typedef typename Params::HistogramParams HistogramParams;
806 
808  typedef typename Params::ValueType ValueType;
809 
811  typedef typename Params::BinningAnalysisParamsType BinningAnalysisParamsType;
814 
816  static constexpr int NumTrackValuesCTime = Params::NumTrackValues;
818  static constexpr int NumLevelsCTime = Params::NumLevels;
819 
821  typedef typename Params::Result ResultType;
822 
827  ValueCalculator,
828  LoggerType,
829  BaseHistogramType
831 
832 private:
833 
834  ValueHistogramMHRWStatsCollectorType value_histogram;
835 
836  BinningAnalysisType binning_analysis;
837 
838  LoggerType & logger;
839 
840  ResultType result;
841 
842 public:
843 
844  ValueHistogramWithBinningMHRWStatsCollector(HistogramParams histogram_params,
845  const ValueCalculator & vcalc,
846  int num_levels,
847  LoggerType & logger_)
848  : value_histogram(histogram_params, vcalc, logger_),
849  binning_analysis((int)histogram_params.num_bins, num_levels, logger_),
850  logger(logger_),
851  result(histogram_params, binning_analysis)
852  {
853  logger.longdebug("ValueHistogramWithBinningMHRWStatsCollector", "constructor()");
854  }
855 
857  inline const BaseHistogramType & histogram() const
858  {
859  return value_histogram.histogram();
860  }
861 
862  inline const BinningAnalysisType & getBinningAnalysis() const
863  {
864  return binning_analysis;
865  }
866 
873  inline const ResultType & getResult() const
874  {
875  return result;
876  }
877 
882  {
883  const BaseHistogramType & hist = histogram();
884  return hist.bins.template cast<CountRealAvgType>() /
885  (CountRealAvgType)(hist.bins.sum() + hist.off_chart);
886  }
887 
898  inline ResultType stealResult()
899  {
900  return std::move(result);
901  }
902 
903  // stats collector part
904 
906  inline void init()
907  {
908  value_histogram.init();
909  }
911  inline void thermalizingDone()
912  {
913  value_histogram.thermalizingDone();
914  }
916  inline void done()
917  {
918  logger.longdebug("ValueHistogramWithBinningMHRWStatsCollector::done()", "finishing up ...");
919 
920  value_histogram.template done<false>();
921 
922  //
923  // Determine the error bars from the binning analysis. Remember, the binning analysis was
924  // applied to each of the indicator functions "chi(value) = (value in bin # i) ? 1 : 0"
925  //
926  // The total number of samples is simply h.bins.sum()+h.off_chart; this is the relevant
927  // coefficient to calculate the bin means needed by binning_analysis.calcErrorLevels(). Indeed,
928  // in this way we really get the averaged observed value of each indicator function for each
929  // value interval.
930  //
931  const BaseHistogramType & h = value_histogram.histogram();
932  result.histogram.params = h.params;
933  CountRealAvgType numsamples = h.bins.sum() + h.off_chart;
934  result.histogram.bins = h.bins.template cast<CountRealAvgType>() / numsamples;
935  result.error_levels = binning_analysis.calcErrorLevels(result.histogram.bins);
936  result.histogram.delta = result.error_levels.col(binning_analysis.numLevels()).template cast<CountRealAvgType>();
937  result.histogram.off_chart = h.off_chart / numsamples;
938 
939  result.converged_status = binning_analysis.determineErrorConvergence(result.error_levels);
940 
941  logger.debug("ValueHistogramWithBinningMHRWStatsCollector", [&,this](std::ostream & str) {
942  str << "Binning analysis: bin sqmeans at different binning levels are:\n"
943  << binning_analysis.getBinSqmeans() << "\n"
944  << "\t-> so the error bars at different binning levels are:\n"
945  << result.error_levels << "\n"
946  << "\t-> convergence analysis: \n";
947  result.dumpConvergenceAnalysis(str);
948  str << "\t... and just for you, here is the final histogram:\n" << result.histogram.prettyPrint() << "\n";
949  });
950  }
951 
953  template<typename CountIntType2, typename PointType, typename LLHValueType, typename MHRandomWalk>
954  inline void rawMove(CountIntType2 k, bool is_thermalizing, bool is_live_iter, bool accepted,
955  double a, const PointType & newpt, LLHValueType newptval,
956  const PointType & curpt, LLHValueType curptval, MHRandomWalk & mh)
957  {
958  value_histogram.rawMove(k, is_thermalizing, is_live_iter, accepted, a, newpt, newptval, curpt, curptval, mh);
959  }
960 
962  template<typename CountIntType2, typename PointType, typename LLHValueType, typename MHRandomWalk>
963  inline void processSample(CountIntType2 k, CountIntType2 n, const PointType & curpt,
964  LLHValueType curptval, MHRandomWalk & mh)
965  {
966  Eigen::Index histindex = value_histogram.processSample(k, n, curpt, curptval, mh);
967  binning_analysis.processNewValues(
969  histindex,
970  value_histogram.histogram().numBins()
971  )
972  );
973  }
974 
975 };
976 
977 
981 template<typename ValueCalculator_, typename LoggerType = Logger::VacuumLogger,
982  typename HistogramType_ = Histogram<typename ValueCalculator_::ValueType> >
983 inline
986  typename HistogramType_::Params hist_params,
987  ValueCalculator_ valcalc,
988  LoggerType & logger
989  )
990 {
992  std::move(hist_params),
993  std::move(valcalc),
994  logger
995  ) ;
996 }
997 
998 
1002 template<typename ValueCalculator_, typename CountIntType_ = int,
1003  typename CountRealAvgType_ = double, int NumTrackValues_ = Eigen::Dynamic,
1004  int NumLevels_ = Eigen::Dynamic,
1005  typename LoggerType = Tomographer::Logger::VacuumLogger>
1006 inline
1008  ValueHistogramWithBinningMHRWStatsCollectorParams<ValueCalculator_, CountIntType_,
1009  CountRealAvgType_, NumTrackValues_,
1010  NumLevels_>,
1011  LoggerType
1012  >
1015  ValueCalculator_ valcalc,
1016  int num_binning_levels,
1017  LoggerType & logger
1018  )
1019 {
1020  return
1022  ValueHistogramWithBinningMHRWStatsCollectorParams<ValueCalculator_, CountIntType_,
1023  CountRealAvgType_, NumTrackValues_,
1024  NumLevels_>,
1025  LoggerType
1026  >(std::move(hist_params), std::move(valcalc), num_binning_levels, logger) ;
1027 }
1028 
1029 
1030 
1034 template<typename MHRWParamsType_>
1035 class TOMOGRAPHER_EXPORT PredStatusReportMHRWStatsCollector
1036 {
1037 public:
1038  typedef MHRWParamsType_ MHRWParamsType;
1040 
1041 private:
1042  std::function<bool()> pred_fn;
1044 
1045 public:
1046  template<typename PredicateFn, typename SendStatusFn>
1047  PredStatusReportMHRWStatsCollector(PredicateFn && pred_fn_, SendStatusFn && send_status_fn_)
1048  : pred_fn(pred_fn_), send_status_fn(send_status_fn_)
1049  { }
1050 
1051  inline void init() { }
1052  inline void thermalizingDone() { }
1053  inline void done() { }
1054 
1055  template<typename IterCountIntType, typename PointType, typename FnValueType, typename MHRandomWalk>
1056  inline void rawMove(
1057  IterCountIntType k, bool is_thermalizing, bool, bool, double, const PointType &, FnValueType,
1058  const PointType &, FnValueType, MHRandomWalk & rw
1059  )
1060  {
1061  // see if we should provide a status report
1062  // fprintf(stderr, "StatusReportCheck::rawMove(): testing for status report requested!\n");
1063  // only check once per sweep, to speed up things
1064  if ((k % rw.nSweep() == 0) && pred_fn()) {
1065  send_status_fn(
1066  MHRWStatusReportType::createFromRandWalkStatInfo(k, is_thermalizing, rw,
1067  rw.statsCollector(),
1068  rw.mhrwController())
1069  );
1070  }
1071  // fprintf(stderr, "StatusReportCheck::rawMove(): done\n");
1072  }
1073 
1074  template<typename IterCountIntType, typename PointType, typename FnValueType, typename MHRandomWalk>
1075  inline void processSample(IterCountIntType, IterCountIntType, PointType&&, FnValueType, MHRandomWalk &)
1076  {
1077  }
1078 };
1079 
1080 
1084 template<typename MHRWParamsType_, typename ClockType_ =
1085  // exception for gcc/g++ 4.6, which doesn't have steady_clock but has monotonic_clock instead
1086 #if defined(__GNUC__) && __GNUC__ == 4 && __GNUC_MINOR__ <= 6 && !defined(__clang__)
1087  std::chrono::monotonic_clock
1088 #else
1090 #endif
1091  >
1093 {
1094 public:
1095  typedef MHRWParamsType_ MHRWParamsType;
1097 
1098  typedef ClockType_ ClockType;
1099  typedef typename ClockType::time_point TimePoint;
1100  typedef typename ClockType::duration DurationType;
1101 
1102 private:
1103  TimePoint last_status_report;
1104  DurationType interval;
1106 
1107 public:
1108  template<typename DurationInit, typename SendStatusFn>
1109  PeriodicStatusReportMHRWStatsCollector(DurationInit && interval_, SendStatusFn && send_status_fn_)
1110  : interval(interval_), send_status_fn(send_status_fn_)
1111  { }
1112 
1113  inline void init() { }
1114  inline void thermalizingDone() { }
1115  inline void done() { }
1116 
1117  template<typename IterCountIntType, typename PointType, typename FnValueType, typename MHRandomWalk>
1118  inline void rawMove(
1119  IterCountIntType k, bool is_thermalizing, bool, bool, double, const PointType &, FnValueType,
1120  const PointType &, FnValueType, MHRandomWalk & rw
1121  )
1122  {
1123  // only check once per sweep, to speed up things
1124  if ((k % rw.nSweep() == 0)) {
1125 
1126  auto now = ClockType::now();
1127 
1128  if ( interval.count() > 0 && (now - (last_status_report + interval)).count() > 0 ) {
1129  // send in a status report
1130  send_status_fn(
1131  MHRWStatusReportType::createFromRandWalkStatInfo(k, is_thermalizing, rw,
1132  rw.statsCollector(),
1133  rw.mhrwController())
1134  );
1135  last_status_report = now;
1136  }
1137 
1138  }
1139  // fprintf(stderr, "StatusReportCheck::rawMove(): done\n");
1140  }
1141 
1142  template<typename IterCountIntType, typename PointType, typename FnValueType, typename MHRandomWalk>
1143  inline void processSample(IterCountIntType, IterCountIntType, PointType&&, FnValueType, MHRandomWalk &)
1144  {
1145  }
1146 };
1147 
1148 
1149 
1150 
1151 
1152 
1153 
1154 
1155 
1156 //
1157 // Specialize Tomographer::Tools::StatusProvider for our stats collector classes
1158 //
1159 
1160 
1161 namespace Tools {
1162 
1166 template<typename... Args>
1167 struct TOMOGRAPHER_EXPORT StatusProvider<MultipleMHRWStatsCollectors<Args... > >
1168 {
1171 
1172  static constexpr int NumStatColl = MHRWStatsCollector::NumStatColl;
1173 
1174  static constexpr bool CanProvideStatusLine = true;
1175 
1176  template<int I = 0, typename std::enable_if<(I < NumStatColl), bool>::type dummy = true>
1177  static inline std::string getStatusLine(const MHRWStatsCollector * stats)
1178  {
1179  typedef typename std::tuple_element<I, typename MHRWStatsCollector::MHRWStatsCollectorsTupleType>::type
1180  ThisStatsCollector;
1181  return
1183  ? StatusQuery<ThisStatsCollector>::getStatusLine(& stats->template getStatsCollector<I>())
1184  : std::string()),
1185  getStatusLine<I+1>(stats) );
1186  };
1187 
1188  template<int I = 0, typename std::enable_if<(I == NumStatColl), bool>::type dummy = true>
1189  static inline std::string getStatusLine(const MHRWStatsCollector * stats)
1190  {
1191  (void)stats;
1192  return std::string();
1193  }
1194 
1195 private:
1196  static inline std::string _joinnl(std::string a, std::string b) {
1197  if (a.size() && b.size()) {
1198  return std::move(a) + "\n" + std::move(b);
1199  }
1200  return std::move(a) + std::move(b); // one of these guys is empty
1201  }
1202 };
1203 // static members:
1204 template<typename... Args>
1205 constexpr int StatusProvider<MultipleMHRWStatsCollectors<Args... > >::NumStatColl;
1206 template<typename... Args>
1207 constexpr bool StatusProvider<MultipleMHRWStatsCollectors<Args... > >::CanProvideStatusLine;
1208 
1209 
1210 
1214 template<typename ValueCalculator_,
1215  typename LoggerType_,
1216  typename HistogramType_
1217  >
1218 struct TOMOGRAPHER_EXPORT StatusProvider<ValueHistogramMHRWStatsCollector<ValueCalculator_, LoggerType_, HistogramType_> >
1219 {
1221 
1222  static constexpr bool CanProvideStatusLine = true;
1223 
1224  static inline std::string getStatusLine(const MHRWStatsCollector * stats)
1225  {
1226  const int maxbarwidth = 50;
1227 
1228  typedef typename MHRWStatsCollector::HistogramType HistogramType;
1229 
1230  return "Histogram: " + histogramShortBar<HistogramType>(stats->histogram(), true, maxbarwidth);
1231  }
1232 };
1233 // static members:
1234 template<typename ValueCalculator_,
1235  typename LoggerType_,
1236  typename HistogramType_
1237  >
1238 constexpr bool
1240 
1241 
1245 template<typename Params_,
1246  typename LoggerType_
1247  >
1248 struct TOMOGRAPHER_EXPORT StatusProvider<ValueHistogramWithBinningMHRWStatsCollector<Params_, LoggerType_> >
1249 {
1251 
1252  static constexpr bool CanProvideStatusLine = true;
1253 
1254  static inline std::string getStatusLine(const MHRWStatsCollector * stats)
1255  {
1256  const int maxbarwidth = 50;
1257 
1258  typedef typename MHRWStatsCollector::BaseHistogramType BaseHistogramType;
1259  const BaseHistogramType & histogram = stats->histogram();
1260 
1261  // calculate the error bars at different levels, to determine convergence status.
1262  //typedef typename MHRWStatsCollector::BinningAnalysisType BinningAnalysisType;
1263  //typedef typename MHRWStatsCollector::CountRealAvgType CountRealAvgType;
1264  //typedef typename BinningAnalysisType::ValueType ValueType;
1265  const auto& binning_analysis = stats->getBinningAnalysis();
1266  const auto binmeans = stats->binMeans();
1267 
1268  auto error_levels = binning_analysis.calcErrorLevels(binmeans);
1269  auto conv_status = binning_analysis.determineErrorConvergence(error_levels);
1270 
1271  auto conv_summary = BinningErrorBarConvergenceSummary::fromConvergedStatus(conv_status);
1272 
1273  return histogramShortBar<BaseHistogramType>(histogram, true, maxbarwidth)
1274  + Tools::fmts(" err(cnvg/?/x): %d/%d/%d",
1275  (int)conv_summary.n_converged, (int)conv_summary.n_unknown,
1276  (int)conv_summary.n_not_converged);
1277  }
1278 };
1279 // static members:
1280 template<typename Params_,
1281  typename LoggerType_
1282  >
1283 constexpr bool
1285 
1286 
1290 template<typename CountIntType_>
1291 struct TOMOGRAPHER_EXPORT StatusProvider<MHRWMovingAverageAcceptanceRatioStatsCollector<CountIntType_> >
1292 {
1293  static constexpr bool CanProvideStatusLine = true;
1294 
1295  static inline std::string getStatusLine(const MHRWMovingAverageAcceptanceRatioStatsCollector<CountIntType_> * obj)
1296  {
1297  if (obj->hasMovingAverageAcceptanceRatio()) {
1298  return Tools::fmts("acceptance ratio = %.2g (over last %d iter)",
1299  (double)obj->movingAverageAcceptanceRatio(),
1300  (int)obj->bufferSize());
1301  }
1302  return std::string();
1303  }
1304 };
1305 
1306 
1307 
1308 } // namespace Tools
1309 
1310 } // namespace Tomographer
1311 
1312 
1313 
1314 #endif
A StatsCollector which builds a histogram of values calculated with a ValueCalculator for each data s...
A Metropolis-Hastings Random Walk.
Definition: mhrw.h:643
HistogramType histogram
Histogram, already with error bars.
It is unknown whether the error bar has converged or not.
Definition: mhrw_bin_err.h:79
auto canonicalBasisVec(IndexType k, IndexType size) -> const Eigen::CwiseNullaryOp< tomo_internal::can_basis_vec_generator< typename Eigen::internal::traits< Der >::Scalar, IndexType >, Der >
Expression for the k-th canonical basis vector of given dimension.
Definition: eigenutil.h:214
Eigen::Array< CountRealAvgType, Eigen::Dynamic, 1 > binMeans() const
Get the current bin means collected so far.
static BinningErrorBarConvergenceSummary fromConvergedStatus(const Eigen::Ref< const Eigen::ArrayXi > &converged_status)
Construct a summary object from a list of converged status obtained from BinningAnalysis::determineEr...
Definition: mhrw_bin_err.h:131
Params::BaseHistogramType BaseHistogramType
See ValueHistogramWithBinningMHRWStatsCollectorParams::BaseHistogramType .
const MHRWStatsCollector & statsCollector() const
Access the stats collector.
Definition: mhrw.h:766
void done()
Finalize the data collection. Part of the MHRWStatsCollector Interface.
BinningAnalysis< BinningAnalysisParamsType, LoggerType > BinningAnalysisType
The corresponding BinningAnalysis type for this value histogram stats collector.
auto getBinSqmeans() const
Get the raw average of the squared values observed, for each binning level.
Definition: mhrw_bin_err.h:703
The parameters of a Histogram.
Definition: histogram.h:72
void done()
Part of the MHRWStatsCollector Interface.
Base namespace for the Tomographer project.
Definition: densellh.h:45
CountRealAvgType_ CountRealAvgType
Type used to store the averages of the histogram bins.
ValueHistogramMHRWStatsCollector(HistogramParams histogram_params, const ValueCalculator &vcalc, LoggerType &logger)
Simple constructor, initializes with the given values.
const BaseHistogramType & histogram() const
Get the histogram data collected so far. See BaseHistogramType .
MultipleMHRWStatsCollectors< MHRWStatsCollectors... > mkMultipleMHRWStatsCollectors(MHRWStatsCollectors &... mhrwstatscollectors)
Convenience function to create a MultipleMHRWStatsCollector (using template argument deduction) ...
ValueCalculator::ValueType ValueType
The type to use to represent a calculated distance.
HistogramWithErrorBars< typename ValueCalculator::ValueType, CountRealAvgType > HistogramType
The Final Histogram Type (with error bars).
The error bar appears to have converged.
Definition: mhrw_bin_err.h:81
ValueHistogramMHRWStatsCollector< ValueCalculator_, LoggerType, HistogramType_ > mkValueHistogramMHRWStatsCollector(typename HistogramType_::Params hist_params, ValueCalculator_ valcalc, LoggerType &logger)
Helper to easily instantiate a ValueHistogramMHRWStatsCollector.
Eigen::ArrayXi converged_status
Information of convergence status of the error bars (see e.g. BINNING_CONVERGED)
Group template parameters for BinningAnalysis.
Definition: mhrw_bin_err.h:212
Eigen::ArrayXi determineErrorConvergence(const Eigen::Ref< const BinSumSqArray > &error_levels) const
Attempt to determine if the error bars have converged.
Definition: mhrw_bin_err.h:854
Result type of a ValueHistogramWithBinningMHRWStatsCollector.
HistogramType::Params HistogramParams
The corresponding histogram params type.
ValueHistogramWithBinningMHRWStatsCollectorResult(HistogramParams p, const BinningAnalysisType &b)
Constructor which initializes the fields from the histogram and binning analysis type.
Stores a histogram along with error bars.
Definition: histogram.h:566
Params::HistogramParams HistogramParams
See ValueHistogramWithBinningMHRWStatsCollectorParams::HistogramParams .
void processNewValues(const Eigen::DenseBase< Derived > &vals)
Process new raw samples.
Definition: mhrw_bin_err.h:562
T to_string(T... args)
Provide appropriate operator new() definitions for a structure which has a member of the given stored...
ValueHistogramWithBinningMHRWStatsCollectorResult(const HistogramType &histogram_, const Eigen::DenseBase< EigenDerived1 > &error_levels_, const Eigen::DenseBase< EigenDerived2 > &converged_status_)
Simple constructor with direct initialization of fields.
ValueCalculator_ ValueCalculator
The type which calculates the interesting value. Should be of type interface ValueCalculator Interfac...
const Tools::StaticOrDynamic< int,(NumLevelsCTime==Eigen::Dynamic), NumLevelsCTime > numLevels
The number of levels in the binning analysis.
Definition: mhrw_bin_err.h:428
void thermalizingDone()
Part of the MHRWStatsCollector Interface. No-op.
Logger that discards all messages.
Definition: loggers.h:1149
T make_tuple(T... args)
bool hasMovingAverageAcceptanceRatio() const
Returns true after enough samples have been collected to fill the internal buffer.
BinningErrorBarConvergenceSummary errorBarConvergenceSummary() const
A summary of the convergence status of the binning error bars.
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)
Part of the MHRWStatsCollector Interface. No-op.
A "stats collector" which produces status reports whenever a predicate evaluates to true...
HistogramType::Params HistogramParams
Structure which holds the parameters of the histogram we&#39;re recording.
const MHRWController & mhrwController() const
Access the random walk controller.
Definition: mhrw.h:769
T setw(T... args)
Query status from different objects.
STL class.
HistogramType_ ResultType
Required for compliance with Resultable Interface type.
const ResultType & getResult() const
Get the histogram data collected. This method is needed for Resultable Interface compliance.
Traits-like class for ValueHistogramWithBinningMHRWStatsCollector.
double movingAverageAcceptanceRatio() const
Returns the average acceptance ratio as stored in the buffer.
An empty struct used as a ResultType in MultipleMHRWStatsCollector for stats collectors which don&#39;t r...
Collect a histogram of values from a MH random walk, with binning analysis.
BinSumSqArray calcErrorLevels(const Eigen::ArrayBase< Derived > &means) const
Calculate the error bars of samples at different binning levels.
Definition: mhrw_bin_err.h:753
const HistogramType & histogram() const
Get the histogram data collected so far. See HistogramType.
Params::ValueType ValueType
See ValueHistogramWithBinningMHRWStatsCollectorParams::ValueType .
Params::BinningAnalysisParamsType BinningAnalysisParamsType
See ValueHistogramWithBinningMHRWStatsCollectorParams::BinningAnalysisParamsType .
Params::ValueCalculator ValueCalculator
See ValueHistogramWithBinningMHRWStatsCollectorParams::ValueCalculator .
void thermalizingDone()
Part of the MHRWStatsCollector Interface. No-op.
Stores a histogram.
Definition: histogram.h:230
void rawMove(CountIntType k, bool, bool, bool, double, const PointType &, LLHValueType, const PointType &, LLHValueType, MHRandomWalk &)
Part of the MHRWStatsCollector Interface. No-op.
ResultType stealResult()
Retrieve the collected histogram data, in compliance with Resultable Interface compliance.
ValueHistogramWithBinningMHRWStatsCollectorResult< HistogramType, BinningAnalysisParamsType > Result
Result type of the corresponding ValueHistogramWithBinningMHRWStatsCollector.
T str(T... args)
Params::CountRealAvgType CountRealAvgType
See ValueHistogramWithBinningMHRWStatsCollectorParams::CountRealAvgType .
A "stats collector" which produces status reports periodically.
Some C++ utilities, with a tad of C++11 tricks.
T move(T... args)
void dumpConvergenceAnalysis(std::ostream &stream) const
Dump values, error bars and convergence status in verbose, debug-human-readable form into ostream...
HistogramType_ HistogramType
The type of the histogram. Usually a Histogram with ValueType range type.
CountIntType nSweep() const
Number of iterations in a sweep.
Definition: mhrw.h:779
The error bar appears NOT to have converged.
Definition: mhrw_bin_err.h:83
T size(T... args)
Status Report structure representing the status of a MHRandomWalk.
Definition: mhrw.h:1142
void init()
Part of the MHRWStatsCollector Interface. Initializes the histogram to zeros.
std::string fmts(const char *fmt,...)
printf- format to a std::string
Definition: fmt.h:128
Template, specializable class to provide status reports for different objects.
std::string dumpConvergenceAnalysis() const
Dump values, error bars and convergence status in verbose, debug-human-readable form as string...
Params::CountIntType CountIntType
See ValueHistogramWithBinningMHRWStatsCollectorParams::CountIntType .
BinningAnalysisParams< ValueType, NumTrackValues, NumLevels, false, CountIntType > BinningAnalysisParamsType
The relevant BinningAnalysis parameters for us.
Brief summary of convergence of error bars from a binning analysis.
Definition: mhrw_bin_err.h:107
ValueCalculator_ ValueCalculator
The type of the ValueCalculator Interface which calculates the value of which we&#39;re collecting a hist...
ResultType stealResult()
Retrieve the final histogram data, in compliance with Resultable Interface compliance.
void processSample(CountIntType2 k, CountIntType2 n, const PointType &curpt, LLHValueType curptval, MHRandomWalk &mh)
Part of the MHRWStatsCollector Interface. Records the sample in the histogram.
LoggerType_ LoggerType
Somewhere where this object may log what it&#39;s doing.
ValueHistogramWithBinningMHRWStatsCollectorResult()
Simple default constructor (e.g. to use as std::vector<Result>).
ValueHistogramWithBinningMHRWStatsCollector< ValueHistogramWithBinningMHRWStatsCollectorParams< ValueCalculator_, CountIntType_, CountRealAvgType_, NumTrackValues_, NumLevels_ >, LoggerType > mkValueHistogramWithBinningMHRWStatsCollector(HistogramParams< typename ValueCalculator_::ValueType > hist_params, ValueCalculator_ valcalc, int num_binning_levels, LoggerType &logger)
Helper to easily instantiate a ValueHistogramWithBinningMHRWStatsCollector.
Definitions for Histogram Types.
Long Debug logging level.
Definition: loggers.h:122
Histogram< typename ValueCalculator::ValueType, CountIntType > BaseHistogramType
The Base Histogram Type.
Params::Result ResultType
See ValueHistogramWithBinningMHRWStatsCollectorParams::Result .
Routines for performing a Metropolis-Hastings random walk.
Binning Analysis in a Metropolis-Hastings random walk.
STL class.
void init()
Part of the MHRWStatsCollector Interface. Initializes the histogram to zeros.
const ResultType & getResult() const
Get the final histogram data. This method is needed for Resultable Interface compliance.
#define tomographer_assert(...)
Assertion test macro.
Definition: cxxdefs.h:84
BinningAnalysisParamsType::BinSumSqArray error_levels
Detailed error bars for all binning levels.
ValueHistogramMHRWStatsCollector< ValueCalculator, LoggerType, BaseHistogramType > ValueHistogramMHRWStatsCollectorType
This is the natural ValueHistogramMHRWStatsCollector type on which we&#39;re adding error bars...
Utilities for logging messages.
ValueCalculator::ValueType ValueType
The type of a value calculated by the ValueCalculator Interface.
A simple MHRWStatsCollector interface which combines several stats collectors.
Eigen::Index processSample(CountIntType k, CountIntType n, const PointType &curpt, LLHValueType, MHRandomWalk &)
Part of the MHRWStatsCollector Interface. Records the sample in the histogram.
MultipleMHRWStatsCollectors TrivialMHRWStatsCollector
Trivial, NO-OP stats collector.
CountIntType_ CountIntType
Type used to count the number of hits in each bin.