RunningStatistics.hpp

/*
  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