28 #ifndef TOMOGRAPHER_MHRW_BIN_ERR_H 29 #define TOMOGRAPHER_MHRW_BIN_ERR_H 37 #include <boost/math/constants/constants.hpp> 51 namespace tomo_internal {
55 template<
int NumLevels,
bool calculate>
56 struct helper_samples_size {
58 value = Eigen::Dynamic
61 template<
int NumLevels>
62 struct helper_samples_size<NumLevels,true> {
64 value = (1 << NumLevels)
109 Eigen::Index n_converged_ = 0,
110 Eigen::Index n_unknown_ = 0,
111 Eigen::Index n_unknown_isolated_ = 0,
112 Eigen::Index n_not_converged_ = 0)
114 n_converged(n_converged_),
115 n_unknown(n_unknown_),
116 n_unknown_isolated(n_unknown_isolated_),
117 n_not_converged(n_not_converged_)
122 Eigen::Index n_converged;
123 Eigen::Index n_unknown;
124 Eigen::Index n_unknown_isolated;
125 Eigen::Index n_not_converged;
130 static BinningErrorBarConvergenceSummary
133 const Eigen::Index n_bins = converged_status.size();
135 const Eigen::ArrayXi unkn_arr =
137 const Eigen::ArrayXi conv_arr =
140 const Eigen::Index n_unknown = unkn_arr.cast<Eigen::Index>().count();
141 const Eigen::Index n_converged = conv_arr.cast<Eigen::Index>().count();
142 const Eigen::Index n_not_converged =
148 const Eigen::Index n_unknown_isol = (
149 ( unkn_arr.segment(1,n_bins-2)
150 .cwiseProduct(conv_arr.segment(0,n_bins-2))
151 .cwiseProduct(conv_arr.segment(2,n_bins-2))
154 + unkn_arr(0)*conv_arr(1) + unkn_arr(n_bins-1)*conv_arr(n_bins-2)
157 return BinningErrorBarConvergenceSummary(n_bins, n_converged, n_unknown, n_unknown_isol, n_not_converged);
164 return stream << s.n_converged <<
" converged / " 165 << s.n_unknown <<
" maybe (" << s.n_unknown_isolated <<
" isolated) / " 166 << s.n_not_converged <<
" not converged";
210 template<
typename ValueType_,
int NumTrackValues_ = Eigen::Dynamic,
int NumLevels_ = Eigen::Dynamic,
211 bool StoreBinSums_ =
true,
typename CountIntType_ =
int>
224 static constexpr
int NumTrackValuesCTime = NumTrackValues_;
228 static constexpr
int NumLevelsCTime = NumLevels_;
232 static constexpr
int NumLevelsPlusOneCTime = (NumLevelsCTime == Eigen::Dynamic
234 : (NumLevelsCTime + 1));
257 static constexpr
int SamplesSizeCTime =
265 static constexpr
bool StoreBinSums = StoreBinSums_;
360 template<
typename Params,
typename LoggerType_>
376 static constexpr
int NumTrackValuesCTime = Params::NumTrackValuesCTime;
380 static constexpr
int NumLevelsCTime = Params::NumLevelsCTime;
384 static constexpr
int NumLevelsPlusOneCTime = Params::NumLevelsPlusOneCTime;
388 static constexpr
int SamplesSizeCTime = Params::SamplesSizeCTime;
392 static constexpr
bool StoreBinSums = Params::StoreBinSums;
464 SamplesArray samples;
472 CountIntType n_samples;
478 CountIntType n_flushes;
499 BinSumSqArray bin_sumsq;
519 : numTrackValues(num_track_values_),
520 numLevels(num_levels_),
521 samplesSize(1 << numLevels()),
522 samples((
Eigen::Index)numTrackValues(), (
Eigen::Index)samplesSize()),
524 bin_sum(BinSumArray::Zero((
Eigen::Index)numTrackValues())),
525 bin_sumsq(BinSumSqArray::Zero((
Eigen::Index)numTrackValues(), (
Eigen::Index)numLevels()+1)),
544 helper_reset_bin_sum();
545 bin_sumsq = BinSumSqArray::Zero((Eigen::Index)numTrackValues(), (Eigen::Index)numLevels()+1);
546 logger.longdebug(
"BinningAnalysis::reset()",
"ready to go.");
559 template<
typename Derived>
562 const CountIntType ninbin = n_samples % samplesSize();
570 samples.col(ninbin) = vals;
573 helper_update_bin_sum(samples.col(ninbin));
576 if ( ninbin == samplesSize() - 1 ) {
580 logger.longdebug(
"BinningAnalysis::processNewValues()", [&](
std::ostream & str) {
581 str <<
"n_samples is now " << n_samples <<
"; flushing bins. samplesSize() = " << samplesSize();
587 for (
int level = 0; level <= numLevels(); ++level) {
589 const CountIntType binnedsize = CountIntType(1) << (numLevels()-level);
591 logger.longdebug(
"BinningAnalysis::processNewValues()", [&](
std::ostream & str) {
592 str <<
"Processing binning level = " << level <<
": binnedsize="<<binnedsize
593 <<
"; n_flushes=" << n_flushes <<
"\n";
594 str <<
"\tbinned samples = \n" 595 << samples.block(0,0,(Eigen::Index)numTrackValues(),(Eigen::Index)binnedsize);
598 for (CountIntType ksample = 0; ksample < binnedsize; ++ksample) {
599 bin_sumsq.col(level) += samples.col(ksample).cwiseProduct(samples.col(ksample));
600 if (ksample % 2 == 0 && binnedsize > 1) {
601 samples.col(ksample/2) = boost::math::constants::half<ValueType>() *
602 (samples.col(ksample) + samples.col(ksample+1));
608 logger.longdebug(
"BinningAnalysis::processNewValues()", [&](
std::ostream & str) {
609 str <<
"Flushing #" << n_flushes <<
" done. bin_sum is = \n" << bin_sum <<
"\n" 610 <<
"\tbin_sumsq is = \n" << bin_sumsq <<
"\n";
623 template<
typename CalcValType, TOMOGRAPHER_ENABLED_IF_TMPL(NumTrackValuesCTime == 1)>
635 TOMOGRAPHER_ENABLED_IF(StoreBinSums)
636 inline void helper_reset_bin_sum()
638 bin_sum.value = BinSumArray::Zero((Eigen::Index)numTrackValues());
640 TOMOGRAPHER_ENABLED_IF(!StoreBinSums)
641 inline void helper_reset_bin_sum() { }
643 TOMOGRAPHER_ENABLED_IF(StoreBinSums)
647 bin_sum.value += new_samples;
649 template<
typename Derived,
650 TOMOGRAPHER_ENABLED_IF_TMPL(!StoreBinSums)>
668 TOMOGRAPHER_ENABLED_IF(StoreBinSums)
669 inline auto getBinMeans() const
670 #ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN 671 -> decltype(BinSumArray() / ValueType(n_samples))
674 return bin_sum.value / ValueType(n_samples);
684 #ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN 686 bin_sumsq.cwiseQuotient(n_flushes * Tools::replicated<NumTrackValuesCTime,1>(
688 .transpose().reverse(),
695 return bin_sumsq.cwiseQuotient(n_flushes * Tools::replicated<NumTrackValuesCTime,1>(
697 .transpose().reverse(),
708 TOMOGRAPHER_ENABLED_IF(StoreBinSums)
709 inline const BinSumArray & getBinSum()
const {
return bin_sum.value; }
732 template<
typename Derived>
737 const int n_levels_plus_one = numLevels()+1;
738 const int n_track_values = numTrackValues();
745 getBinSqmeans() - Tools::replicated<1,NumLevelsPlusOneCTime>(
750 ).cwiseMax(0).cwiseQuotient(
752 Tools::replicated<NumTrackValuesCTime,1>(
754 .transpose().reverse(),
759 (Eigen::Index)numTrackValues(), (Eigen::Index)n_levels_plus_one,
777 template<
typename Derived>
782 bin_sumsq.col((Eigen::Index)numLevels()) / ValueType(n_flushes)
784 ).cwiseMax(0).cwiseSqrt() /
std::sqrt(ValueType(n_flushes-1));
798 TOMOGRAPHER_ENABLED_IF(StoreBinSums)
799 inline BinSumSqArray calcErrorLevels()
const {
800 BinSumArray means = getBinMeans();
801 return calcErrorLevels(
std::move(means));
816 TOMOGRAPHER_ENABLED_IF(StoreBinSums)
817 inline BinSumArray calcErrorLastLevel()
const {
818 BinSumArray means = getBinMeans();
819 return calcErrorLastLevel(
std::move(means));
836 Eigen::ArrayXi converged_status(numTrackValues());
841 logger.longdebug(
"BinningAnalysis::determineErrorConvergence", [&](
std::ostream & str) {
842 str <<
"error_levels = \n" << error_levels <<
"\n";
849 if (numLevels() < range-1) {
851 converged_status = Eigen::ArrayXi::Constant((Eigen::Index)numTrackValues(), UNKNOWN_CONVERGENCE);
855 converged_status = Eigen::ArrayXi::Constant((Eigen::Index)numTrackValues(), CONVERGED);
857 const auto & errors = error_levels.col((Eigen::Index)numLevels());
859 for (
int level = numLevels()+1 - range; level < numLevels(); ++level) {
861 const auto & errors_thislevel = error_levels.col((Eigen::Index)level);
863 logger.longdebug(
"BinningAnalysis::determineErrorConvergence", [&](
std::ostream & str) {
864 str <<
"About to study level " << level <<
": at this point, converged_status = \n" 865 << converged_status <<
"\nand errors_thislevel = \n" << errors_thislevel;
868 for (Eigen::Index val_it = 0; val_it < (Eigen::Index)numTrackValues(); ++val_it) {
869 if (errors_thislevel(val_it) >= errors(val_it) &&
870 converged_status(val_it) != NOT_CONVERGED) {
871 converged_status(val_it) = CONVERGED;
872 }
else if (errors_thislevel(val_it) < 0.824 * errors(val_it)) {
873 converged_status(val_it) = NOT_CONVERGED;
874 }
else if ((errors_thislevel(val_it) < 0.9 * errors(val_it)) &&
875 converged_status(val_it) != NOT_CONVERGED) {
876 converged_status(val_it) = UNKNOWN_CONVERGENCE;
884 logger.longdebug(
"BinningAnalysis::determineErrorConvergence", [&](
std::ostream & str) {
885 str <<
"Done. converged_status [0=UNNOWN,1=CONVERGED,2=NOT CONVERGED] = \n" 889 return converged_status;
It is unknown whether the error bar has converged or not.
void processNewValue(const CalcValType val)
Process a new value (if we're tracking a single function only)
Params::BinSumSqArray BinSumSqArray
Type used to store the sum of squares of values at each binning level. See BinningAnalysisParams::Bin...
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...
auto getBinSqmeans() const
Get the raw average of the squared values observed, for each binning level.
Base namespace for the Tomographer project.
The error bar appears to have converged.
Group template parameters for BinningAnalysis.
Eigen::ArrayXi determineErrorConvergence(const Eigen::Ref< const BinSumSqArray > &error_levels) const
Attempt to determine if the error bars have converged.
Eigen::Array< ValueType, NumTrackValuesCTime, 1 > BinSumArray
Type used to store the sum of values.
Params::CountIntType CountIntType
Type used to count the number of samples. See BinningAnalysisParams::CountIntType.
void processNewValues(const Eigen::DenseBase< Derived > &vals)
Process new raw samples.
const Tools::StaticOrDynamic< int,(NumLevelsCTime==Eigen::Dynamic), NumLevelsCTime > numLevels
The number of levels in the binning analysis.
Params::ValueType ValueType
Type of the value(s) for which we are calculating error bars. See BinningAnalysisParams::ValueType.
const Tools::StaticOrDynamic< CountIntType,(SamplesSizeCTime==Eigen::Dynamic),(SamplesSizeCTime< 0 ? 0 :SamplesSizeCTime)> samplesSize
The size of our samples buffer. (avoid using, might change in the future)
CountIntType getNumFlushes() const
Return the number of times the collected samples were flushed.
const CwiseBinaryOp< internal::scalar_product_op< typename Derived::Scalar, typename OtherDerived::Scalar >, const Derived, const OtherDerived > cwiseProduct(const Eigen::ArrayBase< OtherDerived > &other) const
ValueType_ ValueType
Type of the value which we are calculating error bars for. Also the type of the error bars themselves...
void reset()
Reset this object and start again as if freshly constructed.
BinSumSqArray calcErrorLevels(const Eigen::ArrayBase< Derived > &means) const
Calculate the error bars of samples at different binning levels.
Basic utilities for dealing with Eigen matrices and other types.
Params::SamplesArray SamplesArray
Type of the internal buffer for the raw samples. See BinningAnalysisParams::SamplesArray.
Params::BinSumArray BinSumArray
Type used to store the sum of values. See BinningAnalysisParams::BinSumArray.
Some C++ utilities, with a tad of C++11 tricks.
Simple binning analysis for determining error bars.
LoggerType_ LoggerType
Type of the logger we will be logging debugging messages to.
BinningAnalysis(int num_track_values_, int num_levels_, LoggerType &logger_)
Constructor.
Managing the need for specific overrides to operator new() for some types (especially Eigen types) ...
The error bar appears NOT to have converged.
Eigen::Array< ValueType, NumTrackValuesCTime, NumLevelsPlusOneCTime > BinSumSqArray
Type used to store the sum of squares of values at each binning level.
VarValueDecoder< T >::RetType value(const Var &var)
Access the value of the given variable, as a C++ type.
Brief summary of convergence of error bars from a binning analysis.
const BinSumSqArray & getBinSumsq() const
Get the raw sums of the squared values observed, at each binning level.
Eigen::Array< ValueType, NumTrackValuesCTime, SamplesSizeCTime > SamplesArray
Type of the internal buffer for the raw samples.
BinSumArray calcErrorLastLevel(const Eigen::ArrayBase< Derived > &means) const
Calculate the error bar of samples (from the last binning level).
BinningConvergence
Convergence status of binning error bar.
CountIntType_ CountIntType
Type used to count the number of samples. Usually int is fine, except if you are really taking a LOT ...
#define tomographer_assert(...)
Assertion test macro.
const Tools::StaticOrDynamic< int,(NumTrackValuesCTime==Eigen::Dynamic), NumTrackValuesCTime > numTrackValues
The number of functions being tracked/analyzed.
Utilities for logging messages.