histogram.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_HISTOGRAM_H
#define THEORETICA_HISTOGRAM_H
 
#ifndef THEORETICA_NO_PRINT
#include <sstream>
#include <ostream>
#endif
 
#include <vector>
#include "../core/real_analysis.h"
#include "../core/dataset.h"
#include "./statistics.h"
 
 
namespace theoretica {
 
 
    class histogram {
        private:
 
            size_t N {0};
 
            std::vector<unsigned int> bin_counts {};
 
            real range_max {nan()};
 
            real range_min {nan()};
            
            real value_max {nan()};
 
            real value_min {nan()};
 
            real run_average {0};
 
            real run_tss {0};
 
        public:
 
            histogram() = default;
 
            histogram(unsigned int bin_count, real range_min, real range_max)
                : N(0), value_max(-inf()), value_min(-inf()), run_average(0), run_tss(0) {
 
                bin_counts.resize(bin_count);
                this->range_max = range_max;
                this->range_min = range_min;
            }
 
 
            template<typename Dataset, enable_vector<Dataset> = true>
            histogram(const Dataset& data, unsigned int bin_count = 0) {
 
                range_max = theoretica::max(data);
                range_min = theoretica::min(data);
                value_max = range_max;
                value_min = range_min;
                N = data.size();
 
                // Compute mean and TSS
                run_average = stats::mean(data);
                run_tss = stats::total_sum_squares(data);
 
                // Default bin count is sqrt(N)
                bin_counts.resize(
                    bin_count ? bin_count : floor(sqrt(N))
                );
 
                // The histogram contains all the data points by construction
                for (size_t i = 0; i < N; ++i)
                    bin_counts[index(data[i])]++;
            }
 
 
            inline void insert(real x) {
 
                if(x < range_min || x > range_max)
                    return;
 
                // Update average and TSS using Welford's method
                const real tmp = run_average;
                run_average = tmp + (x - tmp) / (N + 1);
                run_tss += (x - tmp) * (x - run_average);
 
                value_max = value_max < x ? x : value_max;
                value_min = value_min > x ? x : value_min;
 
                bin_counts[index(x)]++;
                N++;
            }
 
 
            inline unsigned int index(real x) const {
 
                if(abs(x - range_max) < MACH_EPSILON)
                    return bin_counts.size() - 1;
 
                return floor(
                    (x - range_min) / (range_max - range_min)
                    * bin_counts.size()
                );
            }
 
 
            inline vec2 range() const {
                return vec2({range_min, range_max});
            }
 
 
            // Statistical functions
 
 
            inline unsigned int number() const {
                return N;
            }
 
 
            inline std::vector<unsigned int> bins() const {
                return bin_counts;
            }
 
 
            inline real range_lower() const {
                return range_min;
            }
 
 
            inline real range_upper() const {
                return range_max;
            }
 
 
            inline real max() const {
                return value_max;
            }
 
 
            inline real min() const {
                return value_min;
            }
 
 
            inline real mean() const {
                return run_average;
            }
 
 
            inline real tss() const {
                return run_tss;
            }
 
 
            inline void rebuild(
                const std::vector<unsigned int>& bin_counts,
                const vec2& range, size_t N, real run_average, real run_tss,
                real value_min, real value_max) {
 
                this->bin_counts = bin_counts;
                this->range_min = range[0];
                this->range_max = range[1];
                this->N = N;
                this->run_average = run_average;
                this->run_tss = run_tss;
                this->value_min = value_min;
                this->value_max = value_max;
            }
 
 
            // Operators
 
 
            inline real operator()(real x) {
 
                if (x < range_min || x > range_max)
                    return 0.0;
 
                return bin_counts[index(x)];
            }
 
 
            inline unsigned int operator[](unsigned int i) const {
                return bin_counts[i];
            }
 
 
            // TO-DO Cumulative Distribution Function
 
 
#ifndef THEORETICA_NO_PRINT
 
            inline std::string to_string(
                const std::string& separator = " ",
                bool normalized = true,
                bool lower_extreme = false) const {
 
                if(N == 0)
                    return "";
 
                std::stringstream res;
                const real width = abs(range_max - range_min) / bin_counts.size();
                real mult = 0.5;
 
                if (lower_extreme)
                    mult = 0.0;
 
                for (size_t i = 0; i < bin_counts.size(); ++i) {
 
                    res << (range_min + (i + mult) * width) << separator;
 
                    if (normalized)
                        res << (bin_counts[i] / (real) N) << std::endl;
                    else
                        res << bin_counts[i] << std::endl;
                }
 
                return res.str();
            }
 
 
            inline operator std::string() {
                return to_string();
            }
 
 
            inline friend std::ostream& operator<<(
                std::ostream& out, const histogram& obj) {
                return out << obj.to_string();
            }
 
#endif
 
    };
 
 
    // Statistical functions over elements of a histogram
 
 
    namespace stats {
 
 
        inline real mean(const histogram& h) {
            return h.mean();
        }
 
 
        inline real tss(const histogram& h) {
            return h.tss();
        }
 
 
        inline real variance(const histogram& h) {
 
            if (h.number() <= 1) {
                TH_MATH_ERROR("variance", h.number(), MathError::DivByZero);
                return nan();
            }
 
            return h.tss() / (h.number() - 1);
        }
 
 
        inline real stdev(const histogram& h) {
            return sqrt(variance(h));
        }
    }
 
 
    inline real max(const histogram& h) {
        return h.max();
    }
 
 
    inline real min(const histogram& h) {
        return h.min();
    }
 
}
 
#endif