mhrw_bin_err.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_MHRW_BIN_ERR_H
#define TOMOGRAPHER_MHRW_BIN_ERR_H


#include <tomographer/tools/loggers.h>
#include <tomographer/tools/eigenutil.h> // replicated(), powersOfTwo()
#include <tomographer/tools/cxxutil.h>
#include <tomographer/tools/needownoperatornew.h>

#include <boost/math/constants/constants.hpp>


namespace Tomographer {

namespace tomo_internal {
//
// internal helper to silence "left shift is negative" GCC warnings
//
template<int NumLevels, bool calculate>
struct helper_samples_size {
  enum {
    value = Eigen::Dynamic
  };
};
template<int NumLevels>
struct helper_samples_size<NumLevels,true> {
  enum {
    value = (1 << NumLevels)
  };
};
} // namespace tomo_internal



enum BinningConvergence {
  BINNING_UNKNOWN_CONVERGENCE = 0,
  BINNING_CONVERGED,
  BINNING_NOT_CONVERGED,
};






struct TOMOGRAPHER_EXPORT BinningErrorBarConvergenceSummary {
  BinningErrorBarConvergenceSummary(Eigen::Index n_bins_ = 0,
                                    Eigen::Index n_converged_ = 0,
                                    Eigen::Index n_unknown_ = 0,
                                    Eigen::Index n_unknown_isolated_ = 0,
                                    Eigen::Index n_not_converged_ = 0)
    : n_bins(n_bins_),
      n_converged(n_converged_),
      n_unknown(n_unknown_),
      n_unknown_isolated(n_unknown_isolated_),
      n_not_converged(n_not_converged_)
  {
  }

  Eigen::Index n_bins;
  Eigen::Index n_converged;
  Eigen::Index n_unknown;
  Eigen::Index n_unknown_isolated;
  Eigen::Index n_not_converged;

  static BinningErrorBarConvergenceSummary
    fromConvergedStatus(const Eigen::Ref<const Eigen::ArrayXi> & converged_status)
  {
    const Eigen::Index n_bins = converged_status.size();

    const Eigen::ArrayXi unkn_arr = 
      converged_status.cwiseEqual(BINNING_UNKNOWN_CONVERGENCE).cast<int>();
    const Eigen::ArrayXi conv_arr =
      converged_status.cwiseEqual(BINNING_CONVERGED).cast<int>();

    const Eigen::Index n_unknown = unkn_arr.cast<Eigen::Index>().count();
    const Eigen::Index n_converged = conv_arr.cast<Eigen::Index>().count();
    const Eigen::Index n_not_converged =
      converged_status.cwiseEqual(BINNING_NOT_CONVERGED).cast<Eigen::Index>().count();

    // Little heuristic to see whether the "unknown" converged error bars are isolated or
    // not. Use conv_arr shifted by one in each direction as a mask -- an unconverged
    // error bar is isolated if it is surrounded by converged error bars
    const Eigen::Index n_unknown_isol = (
        ( unkn_arr.segment(1,n_bins-2)
          .cwiseProduct(conv_arr.segment(0,n_bins-2))
          .cwiseProduct(conv_arr.segment(2,n_bins-2))
            ).count()
        // edges
        + unkn_arr(0)*conv_arr(1) + unkn_arr(n_bins-1)*conv_arr(n_bins-2)
        );

    return BinningErrorBarConvergenceSummary(n_bins, n_converged, n_unknown, n_unknown_isol, n_not_converged);
  }
};


inline std::ostream & operator<<(std::ostream & stream, const BinningErrorBarConvergenceSummary & s)
{
  return stream << s.n_converged << " converged / "
                << s.n_unknown << " maybe (" << s.n_unknown_isolated << " isolated) / "
                << s.n_not_converged << " not converged";
}





template<typename ValueType_, int NumTrackValues_ = Eigen::Dynamic, int NumLevels_ = Eigen::Dynamic,
         bool StoreBinSums_ = true, typename CountIntType_ = int>
struct TOMOGRAPHER_EXPORT BinningAnalysisParams
{
  typedef ValueType_ ValueType;
  typedef CountIntType_ CountIntType;

  static constexpr int NumTrackValuesCTime = NumTrackValues_;
  static constexpr int NumLevelsCTime = NumLevels_;
  static constexpr int NumLevelsPlusOneCTime = (NumLevelsCTime == Eigen::Dynamic
                                                ? Eigen::Dynamic
                                                : (NumLevelsCTime + 1));
  static constexpr int SamplesSizeCTime =
    tomo_internal::helper_samples_size<NumLevelsCTime, (NumLevelsCTime > 0 && NumLevelsCTime < 7)>::value;

  static constexpr bool StoreBinSums = StoreBinSums_;

  typedef Eigen::Array<ValueType, NumTrackValuesCTime, SamplesSizeCTime> SamplesArray;

  typedef Eigen::Array<ValueType, NumTrackValuesCTime, 1> BinSumArray;
  typedef Eigen::Array<ValueType, NumTrackValuesCTime, NumLevelsPlusOneCTime> BinSumSqArray;

  // //! Constants for error bar convergence analysis.  
  // enum {
  //   //! Unable to determine whether the error bars have converged.
  //   UNKNOWN_CONVERGENCE = BINNING_UNKNOWN_CONVERGENCE,
  //   //! The error bars appear to have converged.
  //   CONVERGED = BINNING_CONVERGED,
  //   //! The error bars don't seem to have converged.
  //   NOT_CONVERGED = BINNING_NOT_CONVERGED
  // };
  
};

template<typename Params, typename LoggerType_>
class TOMOGRAPHER_EXPORT BinningAnalysis
  // inheriting from this has some advantages over EIGEN_MAKE_ALIGNED_OPERATOR_NEW, such
  // as not needing to explicitly declare the specialization
  // NeedOwnOperatorNew<BinningAnalysis<...> >
  : public Tools::EigenAlignedOperatorNewProvider
{
public:
  typedef typename Params::ValueType ValueType;
  typedef typename Params::CountIntType CountIntType;

  static constexpr int NumTrackValuesCTime = Params::NumTrackValuesCTime;
  static constexpr int NumLevelsCTime = Params::NumLevelsCTime;
  static constexpr int NumLevelsPlusOneCTime = Params::NumLevelsPlusOneCTime;
  static constexpr int SamplesSizeCTime = Params::SamplesSizeCTime;
  static constexpr bool StoreBinSums = Params::StoreBinSums;

  typedef typename Params::SamplesArray SamplesArray;
  typedef typename Params::BinSumArray BinSumArray;
  typedef typename Params::BinSumSqArray BinSumSqArray;

  const Tools::StaticOrDynamic<int, (NumTrackValuesCTime==Eigen::Dynamic), NumTrackValuesCTime> numTrackValues;
  const Tools::StaticOrDynamic<int, (NumLevelsCTime==Eigen::Dynamic), NumLevelsCTime> numLevels;
  const Tools::StaticOrDynamic<CountIntType, (SamplesSizeCTime==Eigen::Dynamic),
                               (SamplesSizeCTime < 0 ? 0 : SamplesSizeCTime)> samplesSize;

  enum {
    UNKNOWN_CONVERGENCE = BINNING_UNKNOWN_CONVERGENCE,
    CONVERGED = BINNING_CONVERGED,
    NOT_CONVERGED = BINNING_NOT_CONVERGED
  };

  typedef LoggerType_ LoggerType;

private:

  // private not picked up by doxygen anyway

  /*   * \brief The array in which we store samples that arrive from the simulation.
   *
   * This array has size \a samplesSize() (for each tracking value). Once this array is
   * filled, it is <em>flushed</em>, i.e. the values are processed and stored as
   * appropriate in \ref bin_sum and \ref bin_sumsq.
   *
   * This array has \a numTrackingValues() rows and \a samplesSize() columns.
   */
  SamplesArray samples;

  // where we store the flushed values
  
  /*  * \brief Number of samples seen.
   *
   * This is equal to the number of times \ref processNewValues() was called.
   */
  CountIntType n_samples;
  /*  * \brief Number of flushes.
   *
   * A flush corresponds to having filled all the samples in the sample vector (of size \a
   * samplesSize()), and pushing new values into \ref bin_sum and \ref bin_sumsq.
   */
  CountIntType n_flushes;
  /*  * \brief Sum of all values seen.
   *
   * This is a column vector of \a numTrackingValues() entries.
   *
   * \note This member is only available if the template parameter \a StoreBinSums is set
   * to \c true.
   *
   * \note values are added to this array as soon as they are seen, not when the samples
   * array is flushed. In particular, if the total number of values is not a multiple of
   * \a samplesSize(), then there will be samples counted into \a bin_sum but not into \a
   * bin_sumsq.
   */
  Tools::StoreIfEnabled<BinSumArray, StoreBinSums> bin_sum;
  /*  * \brief Sum of the squares of all flushed & processed values, at different binning
   * levels.
   *
   * This is a matrix of \a numTrackingValues() rows and <em>numLevels()+1</em>
   * columns.
   *
   */
  BinSumSqArray bin_sumsq;

  //  ! Just a boring logger...
  LoggerType & logger;


public:

  BinningAnalysis(int num_track_values_, int num_levels_, LoggerType & logger_)
    : numTrackValues(num_track_values_),
      numLevels(num_levels_),
      samplesSize(1 << numLevels()),
      samples((Eigen::Index)numTrackValues(), (Eigen::Index)samplesSize()),
      n_flushes(0),
      bin_sum(BinSumArray::Zero((Eigen::Index)numTrackValues())),
      bin_sumsq(BinSumSqArray::Zero((Eigen::Index)numTrackValues(), (Eigen::Index)numLevels()+1)),
      logger(logger_)
  {
    tomographer_assert(Tools::isPositive(numLevels()));
    tomographer_assert(Tools::isPowerOfTwo(samplesSize()));
    tomographer_assert( (1<<numLevels()) == samplesSize() );

    reset();
  }

  inline void reset()
  {
    n_flushes = 0;
    n_samples = 0;
    helper_reset_bin_sum();
    bin_sumsq = BinSumSqArray::Zero((Eigen::Index)numTrackValues(), (Eigen::Index)numLevels()+1);
    logger.longdebug("BinningAnalysis::reset()", "ready to go.");
  }
  
  template<typename Derived>
  inline void processNewValues(const Eigen::DenseBase<Derived> & vals)
  {
    const CountIntType ninbin = n_samples % samplesSize();

    ++n_samples;

    tomographer_assert(vals.rows() == numTrackValues());
    tomographer_assert(vals.cols() == 1);

    // store the new values in the bins  [also if ninbin == 0]
    samples.col(ninbin) = vals;

    // add to our sum of values, if applicable.
    helper_update_bin_sum(samples.col(ninbin));

    // see if we have to flush the bins (equivalent to `ninbin == samples_size()-1`)
    if ( ninbin == samplesSize() - 1 ) {
      
      // we have filled all bins. Flush them. Re-use the beginning of the samples[] array
      // to store the reduced bins while flushing them.
      logger.longdebug("BinningAnalysis::processNewValues()", [&](std::ostream & str) {
          str << "n_samples is now " << n_samples << "; flushing bins. samplesSize() = " << samplesSize();
        });

      // the size of the samples at the current level of binning. Starts at samplesSize,
      // and decreases by half at each higher level.

      for (int level = 0; level <= numLevels(); ++level) {

        const CountIntType binnedsize = CountIntType(1) << (numLevels()-level);

        logger.longdebug("BinningAnalysis::processNewValues()", [&](std::ostream & str) {
            str << "Processing binning level = " << level << ": binnedsize="<<binnedsize
                << "; n_flushes=" << n_flushes << "\n";
            str << "\tbinned samples = \n"
                << samples.block(0,0,(Eigen::Index)numTrackValues(),(Eigen::Index)binnedsize);
          });

        for (CountIntType ksample = 0; ksample < binnedsize; ++ksample) {
          bin_sumsq.col(level) += samples.col(ksample).cwiseProduct(samples.col(ksample));
          if (ksample % 2 == 0 && binnedsize > 1) {
            samples.col(ksample/2) = boost::math::constants::half<ValueType>() *
              (samples.col(ksample) + samples.col(ksample+1));
          }
        }

      }

      logger.longdebug("BinningAnalysis::processNewValues()", [&](std::ostream & str) {
          str << "Flushing #" << n_flushes << " done. bin_sum is = \n" << bin_sum << "\n"
              << "\tbin_sumsq is = \n" << bin_sumsq << "\n";
        });

      ++n_flushes;
    }

  }

  template<typename CalcValType, TOMOGRAPHER_ENABLED_IF_TMPL(NumTrackValuesCTime == 1)>
  inline void processNewValue(const CalcValType val)
  {
    // for a single value
    //processNewValues(Eigen::Map<const Eigen::Array<CalcValType,1,1> >(&val));
    processNewValues(Eigen::Array<CalcValType,1,1>::Constant(val));
  }


  inline Eigen::Index effectiveSampleSize() const
  {
    return samplesSize();
  }


private:

  TOMOGRAPHER_ENABLED_IF(StoreBinSums)
  inline void helper_reset_bin_sum()
  {
    bin_sum.value = BinSumArray::Zero((Eigen::Index)numTrackValues());
  }
  TOMOGRAPHER_ENABLED_IF(!StoreBinSums)
  inline void helper_reset_bin_sum() { }

  TOMOGRAPHER_ENABLED_IF(StoreBinSums)
  inline void helper_update_bin_sum(const Eigen::Ref<const Eigen::Array<ValueType, NumTrackValuesCTime, 1> > &
                                    new_samples)
  {
    bin_sum.value += new_samples;
  }
  template<typename Derived,
           TOMOGRAPHER_ENABLED_IF_TMPL(!StoreBinSums)>
  inline void helper_update_bin_sum(const Eigen::DenseBase<Derived> & ) { }


public:

  // retrieve results.

  inline CountIntType getNumFlushes() const { return n_flushes; }

  TOMOGRAPHER_ENABLED_IF(StoreBinSums)
  inline auto getBinMeans() const
#ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN
    -> decltype(BinSumArray() / ValueType(n_samples))
#endif
  {
    return bin_sum.value / ValueType(n_samples);
  }

  inline auto getBinSqmeans() const
#ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN
    -> decltype(
        bin_sumsq.cwiseQuotient(n_flushes * Tools::replicated<NumTrackValuesCTime,1>(
                                    Tools::powersOfTwo<Eigen::Array<ValueType, NumLevelsPlusOneCTime, 1> >(numLevels()+1)
                                    .transpose().reverse(),
                                    // replicated by:
                                    numTrackValues(), 1
                                    ))
        )
#endif
  {
    return bin_sumsq.cwiseQuotient(n_flushes * Tools::replicated<NumTrackValuesCTime,1>(
                                  Tools::powersOfTwo<Eigen::Array<ValueType, NumLevelsPlusOneCTime, 1> >(numLevels()+1)
                                  .transpose().reverse(),
                                  // replicated by:
                                  numTrackValues(), 1
                                  ));
  }


  TOMOGRAPHER_ENABLED_IF(StoreBinSums)
  inline const BinSumArray & getBinSum() const { return bin_sum.value; }

  inline const BinSumSqArray & getBinSumsq() const {
    return bin_sumsq;
  }


  template<typename Derived>
  inline BinSumSqArray calcErrorLevels(const Eigen::ArrayBase<Derived> & means) const
  {
    tomographer_assert(means.rows() == (Eigen::Index)numTrackValues());
    tomographer_assert(means.cols() == 1);
    const int n_levels_plus_one = numLevels()+1;
    const int n_track_values = numTrackValues();

    return (
        getBinSqmeans() - Tools::replicated<1,NumLevelsPlusOneCTime>(
            means.cwiseProduct(means).template cast<ValueType>(),
            // replicated by:
            1, n_levels_plus_one
            )
        ).cwiseMax(0).cwiseQuotient(
            // divide by the number of samples from which these bin-means were obtained, minus one.
            Tools::replicated<NumTrackValuesCTime,1>(
                Tools::powersOfTwo<Eigen::Array<ValueType, NumLevelsPlusOneCTime, 1> >(n_levels_plus_one)
                .transpose().reverse(),
                // replicated by:
                n_track_values, 1
                ) * n_flushes
            - Eigen::Array<ValueType, NumTrackValuesCTime, NumLevelsPlusOneCTime>::Constant(
                (Eigen::Index)numTrackValues(), (Eigen::Index)n_levels_plus_one,
                1 // the constant...
                )
            ).cwiseSqrt();
  }


  template<typename Derived>
  inline BinSumArray calcErrorLastLevel(const Eigen::ArrayBase<Derived> & means) const {
    tomographer_assert(means.rows() == (Eigen::Index)numTrackValues());
    tomographer_assert(means.cols() == 1);
    return (
        bin_sumsq.col((Eigen::Index)numLevels()) / ValueType(n_flushes)
        - means.cwiseProduct(means).template cast<ValueType>()
        ).cwiseMax(0).cwiseSqrt() / std::sqrt(ValueType(n_flushes-1));
  }
  
  TOMOGRAPHER_ENABLED_IF(StoreBinSums)
  inline BinSumSqArray calcErrorLevels() const {
    BinSumArray means = getBinMeans();
    return calcErrorLevels(std::move(means));
  }

  TOMOGRAPHER_ENABLED_IF(StoreBinSums)
  inline BinSumArray calcErrorLastLevel() const {
    BinSumArray means = getBinMeans();
    return calcErrorLastLevel(std::move(means));
  }
  
  inline Eigen::ArrayXi determineErrorConvergence(const Eigen::Ref<const BinSumSqArray> & error_levels) const
  {
    Eigen::ArrayXi converged_status(numTrackValues()); // RVO will help

    tomographer_assert(error_levels.rows() == (Eigen::Index)numTrackValues());
    tomographer_assert(error_levels.cols() == (Eigen::Index)numLevels() + 1);

    logger.longdebug("BinningAnalysis::determineErrorConvergence", [&](std::ostream & str) {
        str << "error_levels = \n" << error_levels << "\n";
      });

    // verify that indeed the errors have converged. Inspired from ALPS code, see
    // https://alps.comp-phys.org/svn/alps1/trunk/alps/src/alps/alea/simplebinning.h

    const int range = 4;
    if (numLevels() < range-1) {

      converged_status = Eigen::ArrayXi::Constant((Eigen::Index)numTrackValues(), UNKNOWN_CONVERGENCE);

    } else {

      converged_status = Eigen::ArrayXi::Constant((Eigen::Index)numTrackValues(), CONVERGED);

      const auto & errors = error_levels.col((Eigen::Index)numLevels());

      for (int level = numLevels()+1 - range; level < numLevels(); ++level) {

        const auto & errors_thislevel = error_levels.col((Eigen::Index)level);

        logger.longdebug("BinningAnalysis::determineErrorConvergence", [&](std::ostream & str) {
            str << "About to study level " << level << ": at this point, converged_status = \n"
                << converged_status << "\nand errors_thislevel = \n" << errors_thislevel;
          });

        for (Eigen::Index val_it = 0; val_it < (Eigen::Index)numTrackValues(); ++val_it) {
          if (errors_thislevel(val_it) >= errors(val_it) &&
              converged_status(val_it) != NOT_CONVERGED) {
            converged_status(val_it) = CONVERGED;
          } else if (errors_thislevel(val_it) < 0.824 * errors(val_it)) {
            converged_status(val_it) = NOT_CONVERGED;
          } else if ((errors_thislevel(val_it) < 0.9 * errors(val_it)) &&
                     converged_status(val_it) != NOT_CONVERGED) {
            converged_status(val_it) = UNKNOWN_CONVERGENCE;
          }
        }

      }

    }

    logger.longdebug("BinningAnalysis::determineErrorConvergence", [&](std::ostream & str) {
        str << "Done. converged_status [0=UNNOWN,1=CONVERGED,2=NOT CONVERGED] = \n"
            << converged_status;
      });

    return converged_status;
  }


};



template<typename IterCountIntType, typename LocalLoggerType>
inline int sanitizeBinningLevels(int binning_num_levels,
                                 IterCountIntType n_run,
                                 IterCountIntType samples_last_level,
                                 LocalLoggerType & logger)
{
  if (binning_num_levels > 0) {
    // provided manual value
    if (binning_num_levels < 4) {
      logger.warning([&](std::ostream & stream) {
          stream << "You are using few binning levels = " << binning_num_levels
                 << " and the resulting error bars might not be reliable." ;
      });
    }
  } else {
    // choose automatically. Make sure that the last level has
    // ~samples_last_level samples to calculate std deviation.
    binning_num_levels = (int)(
        std::floor(std::log(n_run/samples_last_level)
                   / std::log(2)) + 1e-3
        ) ;
    if (binning_num_levels < 1) {
      binning_num_levels = 1;
    }
    if (binning_num_levels < 4) {
      logger.warning([&](std::ostream & stream) {
          stream << "Because n_run is small, you will be using few binning levels = "
                 << binning_num_levels << ", and the resulting error bars "
                 << "might not be reliable." ;
        });
    }
  }
  const IterCountIntType binning_last_level_num_samples =
    (IterCountIntType) std::ldexp((double)n_run, - binning_num_levels);

  logger.debug([&](std::ostream & stream) {
      stream << "Binning analysis: " << binning_num_levels << " levels, with "
             << binning_last_level_num_samples << " samples at last level";
    });

  // warn if number of samples @ last level is below recommended value
  if (binning_last_level_num_samples < samples_last_level) {
    logger.warning([&](std::ostream & stream) {
        stream << "The number of samples (" << binning_last_level_num_samples
               << ") at the last binning level is below the recommended value ("
               << samples_last_level << ").  Consider increasing n_run "
               << "or decreasing binning_num_levels.";
      });
  }
  return binning_num_levels;
}






} // namespace Tomographer















#endif