/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) Statoil ASA
//  Copyright (C) Ceetron Solutions AS
//
//  ResInsight 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.
//
//  ResInsight 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 at <http://www.gnu.org/licenses/gpl.html>
//  for more details.
//
/////////////////////////////////////////////////////////////////////////////////

#include "RigStatisticsDataCache.h"

#include "RigStatisticsCalculator.h"
#include "RigStatisticsMath.h"

#include <cmath> // Needed for HUGE_VAL on Linux

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigStatisticsDataCache::RigStatisticsDataCache( RigStatisticsCalculator* statisticsCalculator )
    : m_statisticsCalculator( statisticsCalculator )
{
    CVF_ASSERT( m_statisticsCalculator.notNull() );

    clearAllStatistics();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::clearAllStatistics()
{
    m_statsAllTimesteps = StatisticsValues();
    m_statsPrTs.clear();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::minMaxCellScalarValues( double& min, double& max )
{
    if ( !m_statsAllTimesteps.m_isMaxMinCalculated )
    {
        min = HUGE_VAL;
        max = -HUGE_VAL;

        size_t i;
        for ( i = 0; i < m_statisticsCalculator->timeStepCount(); i++ )
        {
            double tsmin, tsmax;
            this->minMaxCellScalarValues( i, tsmin, tsmax );
            if ( tsmin < min ) min = tsmin;
            if ( tsmax > max ) max = tsmax;
        }

        m_statsAllTimesteps.m_minValue           = min;
        m_statsAllTimesteps.m_maxValue           = max;
        m_statsAllTimesteps.m_isMaxMinCalculated = true;
    }

    min = m_statsAllTimesteps.m_minValue;
    max = m_statsAllTimesteps.m_maxValue;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::minMaxCellScalarValues( size_t timeStepIndex, double& min, double& max )
{
    if ( timeStepIndex >= m_statsPrTs.size() )
    {
        m_statsPrTs.resize( timeStepIndex + 1 );
    }

    if ( !m_statsPrTs[timeStepIndex].m_isMaxMinCalculated )
    {
        double tsMin = HUGE_VAL;
        double tsMax = -HUGE_VAL;

        m_statisticsCalculator->minMaxCellScalarValues( timeStepIndex, tsMin, tsMax );

        m_statsPrTs[timeStepIndex].m_minValue = tsMin;
        m_statsPrTs[timeStepIndex].m_maxValue = tsMax;

        m_statsPrTs[timeStepIndex].m_isMaxMinCalculated = true;
    }

    min = m_statsPrTs[timeStepIndex].m_minValue;
    max = m_statsPrTs[timeStepIndex].m_maxValue;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::posNegClosestToZero( double& pos, double& neg )
{
    if ( !m_statsAllTimesteps.m_isClosestToZeroCalculated )
    {
        pos = HUGE_VAL;
        neg = -HUGE_VAL;

        size_t i;
        for ( i = 0; i < m_statisticsCalculator->timeStepCount(); i++ )
        {
            double tsNeg, tsPos;
            this->posNegClosestToZero( i, tsPos, tsNeg );
            if ( tsNeg > neg && tsNeg < 0 ) neg = tsNeg;
            if ( tsPos < pos && tsPos > 0 ) pos = tsPos;
        }

        m_statsAllTimesteps.m_posClosestToZero          = pos;
        m_statsAllTimesteps.m_negClosestToZero          = neg;
        m_statsAllTimesteps.m_isClosestToZeroCalculated = true;
    }

    pos = m_statsAllTimesteps.m_posClosestToZero;
    neg = m_statsAllTimesteps.m_negClosestToZero;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::posNegClosestToZero( size_t timeStepIndex, double& posNearZero, double& negNearZero )
{
    if ( timeStepIndex >= m_statsPrTs.size() )
    {
        m_statsPrTs.resize( timeStepIndex + 1 );
    }

    if ( !m_statsPrTs[timeStepIndex].m_isClosestToZeroCalculated )
    {
        double pos = HUGE_VAL;
        double neg = -HUGE_VAL;

        m_statisticsCalculator->posNegClosestToZero( timeStepIndex, pos, neg );

        m_statsPrTs[timeStepIndex].m_posClosestToZero = pos;
        m_statsPrTs[timeStepIndex].m_negClosestToZero = neg;

        m_statsPrTs[timeStepIndex].m_isClosestToZeroCalculated = true;
    }

    posNearZero = m_statsPrTs[timeStepIndex].m_posClosestToZero;
    negNearZero = m_statsPrTs[timeStepIndex].m_negClosestToZero;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::meanCellScalarValues( double& meanValue )
{
    if ( !m_statsAllTimesteps.m_isMeanCalculated )
    {
        m_statisticsCalculator->meanCellScalarValue( m_statsAllTimesteps.m_meanValue );

        m_statsAllTimesteps.m_isMeanCalculated = true;
    }

    meanValue = m_statsAllTimesteps.m_meanValue;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::meanCellScalarValues( size_t timeStepIndex, double& meanValue )
{
    if ( timeStepIndex >= m_statsPrTs.size() )
    {
        m_statsPrTs.resize( timeStepIndex + 1 );
    }

    if ( !m_statsPrTs[timeStepIndex].m_isMeanCalculated )
    {
        m_statisticsCalculator->meanCellScalarValue( timeStepIndex, m_statsPrTs[timeStepIndex].m_meanValue );
        m_statsPrTs[timeStepIndex].m_isMeanCalculated = true;
    }

    meanValue = m_statsPrTs[timeStepIndex].m_meanValue;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::sumCellScalarValues( double& sumValue )
{
    if ( !m_statsAllTimesteps.m_isValueSumCalculated )
    {
        double aggregatedSum = 0.0;
        for ( size_t i = 0; i < m_statisticsCalculator->timeStepCount(); i++ )
        {
            double valueSum = 0.0;
            this->sumCellScalarValues( i, valueSum );

            aggregatedSum += valueSum;
        }

        m_statsAllTimesteps.m_valueSum             = aggregatedSum;
        m_statsAllTimesteps.m_isValueSumCalculated = true;
    }

    sumValue = m_statsAllTimesteps.m_valueSum;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::sumCellScalarValues( size_t timeStepIndex, double& sumValue )
{
    if ( timeStepIndex >= m_statsPrTs.size() )
    {
        m_statsPrTs.resize( timeStepIndex + 1 );
    }

    if ( !m_statsPrTs[timeStepIndex].m_isValueSumCalculated )
    {
        double valueSum    = 0.0;
        size_t sampleCount = 0;
        m_statisticsCalculator->valueSumAndSampleCount( timeStepIndex, valueSum, sampleCount );
        m_statsPrTs[timeStepIndex].m_valueSum             = valueSum;
        m_statsPrTs[timeStepIndex].m_isValueSumCalculated = true;
    }

    sumValue = m_statsPrTs[timeStepIndex].m_valueSum;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<size_t>& RigStatisticsDataCache::cellScalarValuesHistogram()
{
    computeHistogramStatisticsIfNeeded();

    return m_statsAllTimesteps.m_histogram;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<size_t>& RigStatisticsDataCache::cellScalarValuesHistogram( size_t timeStepIndex )
{
    computeHistogramStatisticsIfNeeded( timeStepIndex );

    return m_statsPrTs[timeStepIndex].m_histogram;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<int>& RigStatisticsDataCache::uniqueCellScalarValues()
{
    computeUniqueValuesIfNeeded();

    return m_statsAllTimesteps.m_uniqueValues;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<int>& RigStatisticsDataCache::uniqueCellScalarValues( size_t timeStepIndex )
{
    computeUniqueValuesIfNeeded( timeStepIndex );

    return m_statsPrTs[timeStepIndex].m_uniqueValues;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::mobileVolumeWeightedMean( size_t timeStepIndex, double& mean )
{
    if ( timeStepIndex >= m_statsPrTs.size() )
    {
        m_statsPrTs.resize( timeStepIndex + 1 );
    }

    if ( !m_statsPrTs[timeStepIndex].m_isVolumeWeightedMeanCalculated )
    {
        m_statisticsCalculator->mobileVolumeWeightedMean( timeStepIndex, m_statsPrTs[timeStepIndex].m_volumeWeightedMean );
    }

    mean = m_statsPrTs[timeStepIndex].m_volumeWeightedMean;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::mobileVolumeWeightedMean( double& mean )
{
    if ( !m_statsAllTimesteps.m_isVolumeWeightedMeanCalculated )
    {
        m_statisticsCalculator->mobileVolumeWeightedMean( m_statsAllTimesteps.m_volumeWeightedMean );

        m_statsAllTimesteps.m_isVolumeWeightedMeanCalculated = true;
    }

    mean = m_statsAllTimesteps.m_volumeWeightedMean;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::p10p90CellScalarValues( double& p10, double& p90 )
{
    computeHistogramStatisticsIfNeeded();

    p10 = m_statsAllTimesteps.m_p10;
    p90 = m_statsAllTimesteps.m_p90;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::p10p90CellScalarValues( size_t timeStepIndex, double& p10, double& p90 )
{
    computeHistogramStatisticsIfNeeded( timeStepIndex );

    p10 = m_statsPrTs[timeStepIndex].m_p10;
    p90 = m_statsPrTs[timeStepIndex].m_p90;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::computeHistogramStatisticsIfNeeded()
{
    if ( m_statsAllTimesteps.m_histogram.size() == 0 )
    {
        double min;
        double max;
        size_t nBins = 100;
        this->minMaxCellScalarValues( min, max );

        RigHistogramCalculator histCalc( min, max, nBins, &m_statsAllTimesteps.m_histogram );

        m_statisticsCalculator->addDataToHistogramCalculator( histCalc );

        m_statsAllTimesteps.m_p10 = histCalc.calculatePercentil( 0.1, RigStatisticsMath::PercentileStyle::SWITCHED );
        m_statsAllTimesteps.m_p90 = histCalc.calculatePercentil( 0.9, RigStatisticsMath::PercentileStyle::SWITCHED );
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::computeHistogramStatisticsIfNeeded( size_t timeStepIndex )
{
    if ( m_statsPrTs[timeStepIndex].m_histogram.size() == 0 )
    {
        double min;
        double max;
        size_t nBins = 100;
        this->minMaxCellScalarValues( timeStepIndex, min, max );

        RigHistogramCalculator histCalc( min, max, nBins, &m_statsPrTs[timeStepIndex].m_histogram );

        m_statisticsCalculator->addDataToHistogramCalculator( timeStepIndex, histCalc );

        m_statsPrTs[timeStepIndex].m_p10 = histCalc.calculatePercentil( 0.1, RigStatisticsMath::PercentileStyle::SWITCHED );
        m_statsPrTs[timeStepIndex].m_p90 = histCalc.calculatePercentil( 0.9, RigStatisticsMath::PercentileStyle::SWITCHED );
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::computeUniqueValuesIfNeeded()
{
    if ( m_statsAllTimesteps.m_uniqueValues.size() == 0 )
    {
        std::set<int> setValues;
        m_statisticsCalculator->uniqueValues( 0, setValues ); // This is a Hack ! Only using first timestep. Ok for
                                                              // Static eclipse results but beware !

        for ( auto val : setValues )
        {
            m_statsAllTimesteps.m_uniqueValues.push_back( val );
        }
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigStatisticsDataCache::computeUniqueValuesIfNeeded( size_t timeStepIndex )
{
    if ( m_statsPrTs[timeStepIndex].m_uniqueValues.size() == 0 )
    {
        std::set<int> setValues;
        m_statisticsCalculator->uniqueValues( timeStepIndex, setValues );

        for ( auto val : setValues )
        {
            m_statsPrTs[timeStepIndex].m_uniqueValues.push_back( val );
        }
    }
}