RunningStatistics.hpp

Go to the documentation of this file.

/*
  Copyright 2024 Equinor ASA.
 
  This file is part of the Open Porous Media project (OPM).
 
  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
 
  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef OPM_RUNNING_STATISTICS_HPP
#define OPM_RUNNING_STATISTICS_HPP
 
#include <cmath>
#include <limits>
#include <optional>
 
namespace Opm {
 
template <typename Scalar>
class RunningStatistics
{
public:
    template <class Serializer>
    void serializeOp(Serializer& serializer)
    {
        serializer(this->sampleSize_);
        serializer(this->min_);
        serializer(this->max_);
        serializer(this->mean_);
        serializer(this->totalVariance_);
    }
 
    static RunningStatistics serializationTestObject()
    {
        auto stat = RunningStatistics{};
 
        stat.sampleSize_ = 12;
        stat.min_ = -static_cast<Scalar>(1);
        stat.max_ =  static_cast<Scalar>(2);
        stat.mean_ = static_cast<Scalar>(0.03);
        stat.totalVariance_ = static_cast<Scalar>(0.4);
 
        return stat;
    }
 
    bool operator==(const RunningStatistics& that) const
    {
        return (this->sampleSize_ == that.sampleSize_)
            && (this->min_ == that.min_)
            && (this->max_ == that.max_)
            && (this->mean_ == that.mean_)
            && (this->totalVariance_ == that.totalVariance_)
            ;
    }
 
    void reset()
    {
        this->sampleSize_ = 0;
        this->min_ = std::numeric_limits<Scalar>::max();
        this->max_ = std::numeric_limits<Scalar>::lowest();
        this->mean_ = Scalar{};
        this->totalVariance_ = Scalar{};
    }
 
    void addSamplePoint(const Scalar x)
    {
        if (x < this->min_) { this->min_ = x; }
        if (x > this->max_) { this->max_ = x; }
 
        // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
 
        ++this->sampleSize_;
 
        const auto d1 = x - this->mean();
        this->mean_ += d1 / this->sampleSize_;
 
        const auto d2 = x - this->mean();
        this->totalVariance_ += d1 * d2;
    }
 
    std::size_t sampleSize() const { return this->sampleSize_; }
 
    Scalar min() const { return this->min_; }
 
    Scalar max() const { return this->max_; }
 
    Scalar mean() const { return this->mean_; }
 
    std::optional<Scalar> stdev() const
    {
        if (this->sampleSize_ < 2) {
            return {};
        }
 
        using std::sqrt;
        return sqrt(this->totalVariance_ / (this->sampleSize_ - 1));
    }
 
private:
    std::size_t sampleSize_{};
 
    Scalar min_ { std::numeric_limits<Scalar>::max() };
 
    Scalar max_ { std::numeric_limits<Scalar>::lowest() };
 
    Scalar mean_{};
 
    Scalar totalVariance_{};
};
 
} // namespace Opm
 
#endif // OPM_RUNNING_STATISTICS_HPP