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) 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_BIN_ERR_H
#define TOMOGRAPHER_MHRW_BIN_ERR_H


#include <tomographer/qit/util.h>
#include <tomographer/tools/util.h>
#include <tomographer/tools/loggers.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



template<typename ValueType_, int NumTrackValues_ = Eigen::Dynamic, int NumLevels_ = Eigen::Dynamic,
         bool StoreBinSums_ = true, typename CountIntType_ = int>
struct 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;

  enum {
    UNKNOWN_CONVERGENCE = 0,
    CONVERGED,
    NOT_CONVERGED
  };
  
};

template<typename Params, typename LoggerType_>
class BinningAnalysis
{
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::static_or_dynamic<int, NumTrackValuesCTime> num_track_values;
  const Tools::static_or_dynamic<int, NumLevelsCTime> num_levels;
  const Tools::static_or_dynamic<CountIntType, SamplesSizeCTime> samples_size;

  enum {
    UNKNOWN_CONVERGENCE = Params::UNKNOWN_CONVERGENCE,
    CONVERGED = Params::CONVERGED,
    NOT_CONVERGED = Params::NOT_CONVERGED
  };

  typedef LoggerType_ LoggerType;

private:

  SamplesArray samples;

  // where we store the flushed values
  
  CountIntType n_samples;
  CountIntType n_flushes;
  Tools::store_if_enabled<BinSumArray, StoreBinSums> bin_sum;
  BinSumSqArray bin_sumsq;

  LoggerType & logger;

public:

  BinningAnalysis(int num_track_values_, int num_levels_, LoggerType & logger_)
    : num_track_values(num_track_values_),
      num_levels(num_levels_),
      samples_size(1 << num_levels()),
      samples(num_track_values(), samples_size()),
      n_flushes(0),
      bin_sum(BinSumArray::Zero(num_track_values())),
      bin_sumsq(BinSumSqArray::Zero(num_track_values(), num_levels()+1)),
      logger(logger_)
  {
    assert(Tools::is_positive(num_levels()));
    assert(Tools::is_power_of_two(samples_size()));
    assert( (1<<num_levels()) == samples_size() );

    reset();
  }

  inline void reset()
  {
    n_flushes = 0;
    n_samples = 0;
    helper_reset_bin_sum();
    bin_sumsq = BinSumSqArray::Zero(num_track_values(), num_levels()+1);
    logger.longdebug("BinningAnalysis::reset()", "ready to go.");
  }
  
  template<typename Derived>
  inline void process_new_values(const Eigen::DenseBase<Derived> & vals)
  {
    const int ninbin = n_samples % samples_size();

    ++n_samples;

    eigen_assert(vals.rows() == num_track_values());
    eigen_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 == samples_size() - 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::process_new_values()", [&](std::ostream & str) {
          str << "n_samples is now " << n_samples << "; flushing bins. samples_size() = " << samples_size();
        });

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

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

        const int binnedsize = 1 << (num_levels()-level);

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

        for (int 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::process_new_values()", [&](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 process_new_value(const CalcValType val)
  {
    // for a single value
    //process_new_values(Eigen::Map<const Eigen::Array<CalcValType,1,1> >(&val));
    process_new_values(Eigen::Array<CalcValType,1,1>::Constant(val));
  }



private:

  TOMOGRAPHER_ENABLED_IF(StoreBinSums)
  inline void helper_reset_bin_sum()
  {
    bin_sum.value = BinSumArray::Zero(num_track_values());
  }
  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 get_n_flushes() const { return n_flushes; }

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

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


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

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


  template<typename Derived>
  inline BinSumSqArray calc_error_levels(const Eigen::ArrayBase<Derived> & means) const
  {
    eigen_assert(means.rows() == num_track_values());
    eigen_assert(means.cols() == 1);
    const int n_levels_plus_one = num_levels()+1;
    const int n_track_values = num_track_values();

    return (
        get_bin_sqmeans() - 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.
            replicated<NumTrackValuesCTime,1>(
                powers_of_two<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(
                num_track_values(), num_levels()+1,
                1 // the constant...
                )
            ).cwiseSqrt();
  }


  template<typename Derived>
  inline BinSumArray calc_error_lastlevel(const Eigen::ArrayBase<Derived> & means) const {
    eigen_assert(means.rows() == num_track_values());
    eigen_assert(means.cols() == 1);
    return (
        bin_sumsq.col(num_levels()) / ValueType(n_flushes) - means.cwiseProduct(means).template cast<ValueType>()
        ).cwiseMax(0).cwiseSqrt() / std::sqrt(ValueType(n_flushes-1));
  }
  
  TOMOGRAPHER_ENABLED_IF(StoreBinSums)
  inline BinSumSqArray calc_error_levels() const {
    BinSumArray means = get_bin_means();
    return calc_error_levels(means);
  }

  TOMOGRAPHER_ENABLED_IF(StoreBinSums)
  inline BinSumArray calc_error_lastlevel() const {
    BinSumArray means = get_bin_means();
    return calc_error_lastlevel(means);
  }
  
  inline Eigen::ArrayXi determine_error_convergence(const Eigen::Ref<const BinSumSqArray> & error_levels) const
  {
    Eigen::ArrayXi converged_status(num_track_values()); // RVO will help

    eigen_assert(error_levels.rows() == num_track_values());
    eigen_assert(error_levels.cols() == num_levels() + 1);

    logger.longdebug("BinningAnalysis::determine_error_convergence", [&](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 (num_levels() < range-1) {

      converged_status = Eigen::ArrayXi::Constant(num_track_values(), UNKNOWN_CONVERGENCE);

    } else {

      converged_status = Eigen::ArrayXi::Constant(num_track_values(), CONVERGED);

      const auto & errors = error_levels.col(num_levels());

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

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

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

        for (int val_it = 0; val_it < num_track_values(); ++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::determine_error_convergence", [&](std::ostream & str) {
        str << "Done. converged_status [0=UNNOWN,1=CONVERGED,2=NOT CONVERGED] = \n"
            << converged_status;
      });

    return converged_status;
  }
};



} // namespace Tomographer















#endif