/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2015-     Statoil ASA
//  Copyright (C) 2015-     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 "RigWellLogCurveData.h"

#include "RiaCurveDataTools.h"

#include "RiaEclipseUnitTools.h"

#include "cvfAssert.h"
#include "cvfMath.h"

#include <cmath>

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigWellLogCurveData::RigWellLogCurveData()
{
    m_isExtractionCurve = false;
    m_depthUnit         = RiaDefines::UNIT_METER;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigWellLogCurveData::~RigWellLogCurveData() {}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigWellLogCurveData::setValuesAndMD( const std::vector<double>& xValues,
                                          const std::vector<double>& measuredDepths,
                                          RiaDefines::DepthUnitType  depthUnit,
                                          bool                       isExtractionCurve )
{
    CVF_ASSERT( xValues.size() == measuredDepths.size() );

    m_xValues        = xValues;
    m_measuredDepths = measuredDepths;
    m_tvDepths.clear();
    m_depthUnit = depthUnit;

    // Disable depth value filtering is intended to be used for
    // extraction curve data
    m_isExtractionCurve = isExtractionCurve;

    calculateIntervalsOfContinousValidValues();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigWellLogCurveData::setValuesWithTVD( const std::vector<double>& xValues,
                                            const std::vector<double>& measuredDepths,
                                            const std::vector<double>& tvDepths,
                                            RiaDefines::DepthUnitType  depthUnit,
                                            bool                       isExtractionCurve )
{
    CVF_ASSERT( xValues.size() == measuredDepths.size() );

    m_xValues        = xValues;
    m_measuredDepths = measuredDepths;
    m_tvDepths       = tvDepths;
    m_depthUnit      = depthUnit;

    m_isExtractionCurve = isExtractionCurve;

    calculateIntervalsOfContinousValidValues();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigWellLogCurveData::xValues() const
{
    return m_xValues;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigWellLogCurveData::measuredDepths() const
{
    return m_measuredDepths;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigWellLogCurveData::tvDepths() const
{
    return m_tvDepths;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RigWellLogCurveData::xPlotValues() const
{
    std::vector<double> filteredValues;
    RiaCurveDataTools::getValuesByIntervals( m_xValues, m_intervalsOfContinousValidValues, &filteredValues );

    return filteredValues;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RigWellLogCurveData::trueDepthPlotValues( RiaDefines::DepthUnitType destinationDepthUnit ) const
{
    std::vector<double> filteredValues;
    if ( m_tvDepths.size() )
    {
        if ( destinationDepthUnit == m_depthUnit )
        {
            RiaCurveDataTools::getValuesByIntervals( m_tvDepths, m_intervalsOfContinousValidValues, &filteredValues );
        }
        else
        {
            std::vector<double> convertedValues = convertDepthValues( destinationDepthUnit, m_tvDepths );
            RiaCurveDataTools::getValuesByIntervals( convertedValues, m_intervalsOfContinousValidValues, &filteredValues );
        }
    }

    return filteredValues;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RigWellLogCurveData::measuredDepthPlotValues( RiaDefines::DepthUnitType destinationDepthUnit ) const
{
    std::vector<double> filteredValues;

    if ( destinationDepthUnit == m_depthUnit )
    {
        RiaCurveDataTools::getValuesByIntervals( m_measuredDepths, m_intervalsOfContinousValidValues, &filteredValues );
    }
    else
    {
        std::vector<double> convertedValues = convertDepthValues( destinationDepthUnit, m_measuredDepths );
        RiaCurveDataTools::getValuesByIntervals( convertedValues, m_intervalsOfContinousValidValues, &filteredValues );
    }

    return filteredValues;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<size_t, size_t>> RigWellLogCurveData::polylineStartStopIndices() const
{
    return RiaCurveDataTools::computePolyLineStartStopIndices( m_intervalsOfContinousValidValues );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::ref<RigWellLogCurveData> RigWellLogCurveData::calculateResampledCurveData( double newMeasuredDepthStepSize ) const
{
    std::vector<double> xValues;
    std::vector<double> measuredDepths;

    bool                isTvDepthsAvailable = false;
    std::vector<double> tvDepths;

    if ( m_tvDepths.size() > 0 ) isTvDepthsAvailable = true;

    if ( m_measuredDepths.size() > 0 )
    {
        double currentMd = m_measuredDepths[0];

        size_t segmentStartIdx = 0;
        while ( segmentStartIdx < m_measuredDepths.size() - 1 )
        {
            double segmentStartMd = m_measuredDepths[segmentStartIdx];
            double segmentEndMd   = m_measuredDepths[segmentStartIdx + 1];
            double segmentStartX  = m_xValues[segmentStartIdx];
            double segmentEndX    = m_xValues[segmentStartIdx + 1];

            double segmentStartTvd = 0.0;
            double segmentEndTvd   = 0.0;
            if ( isTvDepthsAvailable )
            {
                segmentStartTvd = m_tvDepths[segmentStartIdx];
                segmentEndTvd   = m_tvDepths[segmentStartIdx + 1];
            }

            while ( currentMd <= segmentEndMd )
            {
                measuredDepths.push_back( currentMd );
                double endWeight = ( currentMd - segmentStartMd ) / ( segmentEndMd - segmentStartMd );
                xValues.push_back( ( 1.0 - endWeight ) * segmentStartX + endWeight * segmentEndX );

                // The tvd calculation is a simplification. We should use the wellpath, as it might have a better
                // resolution, and have a none-linear shape This is much simpler, and possibly accurate enough ?

                if ( isTvDepthsAvailable )
                {
                    tvDepths.push_back( ( 1.0 - endWeight ) * segmentStartTvd + endWeight * segmentEndTvd );
                }

                currentMd += newMeasuredDepthStepSize;
            }

            segmentStartIdx++;
        }
    }

    cvf::ref<RigWellLogCurveData> reSampledData = new RigWellLogCurveData;

    if ( isTvDepthsAvailable )
    {
        reSampledData->setValuesWithTVD( xValues, measuredDepths, tvDepths, m_depthUnit, true );
    }
    else
    {
        reSampledData->setValuesAndMD( xValues, measuredDepths, m_depthUnit, m_isExtractionCurve );
    }

    return reSampledData;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigWellLogCurveData::calculateIntervalsOfContinousValidValues()
{
    std::vector<std::pair<size_t, size_t>> intervalsOfValidValues =
        RiaCurveDataTools::calculateIntervalsOfValidValues( m_xValues, false );

    m_intervalsOfContinousValidValues.clear();

    if ( !m_isExtractionCurve )
    {
        m_intervalsOfContinousValidValues = intervalsOfValidValues;
    }
    else
    {
        size_t intervalsCount = intervalsOfValidValues.size();
        for ( size_t intIdx = 0; intIdx < intervalsCount; intIdx++ )
        {
            std::vector<std::pair<size_t, size_t>> depthValuesIntervals;
            splitIntervalAtEmptySpace( m_measuredDepths,
                                       intervalsOfValidValues[intIdx].first,
                                       intervalsOfValidValues[intIdx].second,
                                       &depthValuesIntervals );

            for ( size_t dvintIdx = 0; dvintIdx < depthValuesIntervals.size(); dvintIdx++ )
            {
                m_intervalsOfContinousValidValues.push_back( depthValuesIntervals[dvintIdx] );
            }
        }
    }
}

//--------------------------------------------------------------------------------------------------
/// Splits the start stop interval between cells that are not close enough.
//--------------------------------------------------------------------------------------------------
void RigWellLogCurveData::splitIntervalAtEmptySpace( const std::vector<double>&              depthValues,
                                                     size_t                                  startIdx,
                                                     size_t                                  stopIdx,
                                                     std::vector<std::pair<size_t, size_t>>* intervals )
{
    CVF_ASSERT( intervals );

    CVF_ASSERT( startIdx < stopIdx );

    if ( stopIdx - startIdx == 1 )
    {
        intervals->push_back( std::make_pair( startIdx, stopIdx ) );
        return;
    }

    // !! TODO: Find a reasonable tolerance
    const double depthDiffTolerance = 0.1;

    // Find intervals containing depth values that should be connected:
    //
    // vIdx = 0 is the first point of a well, usually outside of the model. Further depth values are
    // organized in pairs of depths (in and out of a cell), and sometimes the depths varies slightly. If
    // the distance between a depth pair is larger than the depthDiffTolerance, the two sections will be split
    // into two intervals.
    //
    // The first pair is located at vIdx = 1 & 2. If startIdx = 0, an offset of 1 is added to vIdx, to access
    // that pair in the loop. If startIdx = 1 (can happen if the start point is inside of the model and invalid),
    // the offset is not needed.

    size_t intervalStartIdx = startIdx;
    size_t offset           = 1 - startIdx % 2;
    for ( size_t vIdx = startIdx + offset; vIdx < stopIdx; vIdx += 2 )
    {
        if ( cvf::Math::abs( depthValues[vIdx + 1] - depthValues[vIdx] ) > depthDiffTolerance )
        {
            intervals->push_back( std::make_pair( intervalStartIdx, vIdx ) );
            intervalStartIdx = vIdx + 1;
        }
    }

    if ( intervalStartIdx <= stopIdx )
    {
        intervals->push_back( std::make_pair( intervalStartIdx, stopIdx ) );
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigWellLogCurveData::calculateMDRange( double* minimumDepth, double* maximumDepth ) const
{
    CVF_ASSERT( minimumDepth && maximumDepth );

    double minValue = HUGE_VAL;
    double maxValue = -HUGE_VAL;

    for ( size_t vIdx = 0; vIdx < m_measuredDepths.size(); vIdx++ )
    {
        double value = m_measuredDepths[vIdx];

        if ( value < minValue )
        {
            minValue = value;
        }

        if ( value > maxValue )
        {
            maxValue = value;
        }
    }

    if ( maxValue >= minValue )
    {
        *minimumDepth = minValue;
        *maximumDepth = maxValue;

        return true;
    }

    return false;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigWellLogCurveData::calculateTVDRange( double* minTVD, double* maxTVD ) const
{
    CVF_ASSERT( minTVD && maxTVD );

    double minValue = HUGE_VAL;
    double maxValue = -HUGE_VAL;

    for ( size_t vIdx = 0; vIdx < m_tvDepths.size(); vIdx++ )
    {
        double value = m_tvDepths[vIdx];

        if ( value < minValue )
        {
            minValue = value;
        }

        if ( value > maxValue )
        {
            maxValue = value;
        }
    }

    if ( maxValue >= minValue )
    {
        *minTVD = minValue;
        *maxTVD = maxValue;

        return true;
    }

    return false;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RiaDefines::DepthUnitType RigWellLogCurveData::depthUnit() const
{
    return m_depthUnit;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RigWellLogCurveData::convertFromMeterToFeet( const std::vector<double>& valuesInMeter )
{
    std::vector<double> valuesInFeet( valuesInMeter.size() );

    for ( size_t i = 0; i < valuesInMeter.size(); i++ )
    {
        valuesInFeet[i] = valuesInMeter[i] * RiaEclipseUnitTools::feetPerMeter();
    }

    return valuesInFeet;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RigWellLogCurveData::convertFromFeetToMeter( const std::vector<double>& valuesInFeet )
{
    std::vector<double> valuesInMeter( valuesInFeet.size() );

    for ( size_t i = 0; i < valuesInFeet.size(); i++ )
    {
        valuesInMeter[i] = valuesInFeet[i] / RiaEclipseUnitTools::feetPerMeter();
    }

    return valuesInMeter;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RigWellLogCurveData::convertDepthValues( RiaDefines::DepthUnitType  destinationDepthUnit,
                                                             const std::vector<double>& values ) const
{
    CVF_ASSERT( destinationDepthUnit != m_depthUnit );

    if ( destinationDepthUnit == RiaDefines::UNIT_METER )
    {
        return convertFromFeetToMeter( values );
    }
    else
    {
        return convertFromMeterToFeet( values );
    }
}