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

#include "RiaLogging.h"
#include "RiaWeightedMeanCalculator.h"
#include "RigActiveCellInfo.h"
#include "RigFractureTransmissibilityEquations.h"

#include "cvfAssert.h"

#include <Eigen/Core>
#include <Eigen/LU>

#include <iomanip>

#include <QDebug>
#include <fstream>

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigTransmissibilityCondenser::RigTransmissibilityCondenser()
    : m_transmissibilityThreshold( 1.0e-9 )
{
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigTransmissibilityCondenser::RigTransmissibilityCondenser( const RigTransmissibilityCondenser& copyFrom )
    : m_neighborTransmissibilities( copyFrom.m_neighborTransmissibilities )
    , m_condensedTransmissibilities( copyFrom.m_condensedTransmissibilities )
    , m_externalCellAddrSet( copyFrom.m_externalCellAddrSet )
    , m_TiiInv( copyFrom.m_TiiInv )
    , m_Tie( copyFrom.m_Tie )
    , m_transmissibilityThreshold( copyFrom.m_transmissibilityThreshold )
{
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigTransmissibilityCondenser& RigTransmissibilityCondenser::operator=( const RigTransmissibilityCondenser& rhs )
{
    m_neighborTransmissibilities  = rhs.m_neighborTransmissibilities;
    m_condensedTransmissibilities = rhs.m_condensedTransmissibilities;
    m_externalCellAddrSet         = rhs.m_externalCellAddrSet;
    m_TiiInv                      = rhs.m_TiiInv;
    m_Tie                         = rhs.m_Tie;
    m_transmissibilityThreshold   = rhs.m_transmissibilityThreshold;

    return *this;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigTransmissibilityCondenser::setTransmissibilityThreshold( double threshold )
{
    m_transmissibilityThreshold = threshold;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigTransmissibilityCondenser::transmissibilityThreshold() const
{
    return m_transmissibilityThreshold;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigTransmissibilityCondenser::addNeighborTransmissibility( CellAddress cell1, CellAddress cell2, double transmissibility )
{
    if ( transmissibility < m_transmissibilityThreshold ) return;

    m_condensedTransmissibilities.clear();
    m_externalCellAddrSet.clear();
    if ( cell1 < cell2 )
        m_neighborTransmissibilities[cell1][cell2] += transmissibility;
    else
        m_neighborTransmissibilities[cell2][cell1] += transmissibility;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::set<RigTransmissibilityCondenser::CellAddress> RigTransmissibilityCondenser::externalCells()
{
    if ( m_externalCellAddrSet.empty() )
    {
        calculateCondensedTransmissibilities();
    }

    return m_externalCellAddrSet;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigTransmissibilityCondenser::condensedTransmissibility( CellAddress externalCell1, CellAddress externalCell2 )
{
    CAF_ASSERT( !( externalCell1 == externalCell2 ) );

    if ( m_condensedTransmissibilities.empty() )
    {
        calculateCondensedTransmissibilities();
    }

    if ( externalCell2 < externalCell1 ) std::swap( externalCell1, externalCell2 );

    const auto& adrToAdrTransMapPair = m_condensedTransmissibilities.find( externalCell1 );
    if ( adrToAdrTransMapPair != m_condensedTransmissibilities.end() )
    {
        const auto& adrTransPair = adrToAdrTransMapPair->second.find( externalCell2 );
        if ( adrTransPair != adrToAdrTransMapPair->second.end() )
        {
            return adrTransPair->second;
        }
    }
    return 0.0;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::map<size_t, double>
    RigTransmissibilityCondenser::scaleMatrixToFracTransByMatrixWellDP( const RigActiveCellInfo*   actCellInfo,
                                                                        double                     currentWellPressure,
                                                                        const std::vector<double>& currentMatrixPressures,
                                                                        double*                    minPressureDrop,
                                                                        double*                    maxPressureDrop )
{
    std::map<size_t, double> originalLumpedMatrixToFractureTrans; // Sum(T_mf)

    double epsilonDeltaPressure = 1.0e-6;

    double minNonZeroDeltaPressure = std::numeric_limits<double>::infinity();
    double maxNonZeroDeltaPressure = -std::numeric_limits<double>::infinity();

    for ( auto it = m_neighborTransmissibilities.begin(); it != m_neighborTransmissibilities.end(); ++it )
    {
        if ( it->first.m_cellIndexSpace == CellAddress::STIMPLAN )
        {
            for ( auto jt = it->second.begin(); jt != it->second.end(); ++jt )
            {
                if ( jt->first.m_cellIndexSpace == CellAddress::ECLIPSE )
                {
                    size_t globalMatrixCellIdx = jt->first.m_globalCellIdx;
                    size_t eclipseResultIndex  = actCellInfo->cellResultIndex( globalMatrixCellIdx );
                    CVF_ASSERT( eclipseResultIndex < currentMatrixPressures.size() );
                    double unsignedDeltaPressure =
                        std::abs( currentMatrixPressures[eclipseResultIndex] - currentWellPressure );
                    double nonZeroDeltaPressure = std::max( epsilonDeltaPressure, unsignedDeltaPressure );
                    maxNonZeroDeltaPressure     = std::max( maxNonZeroDeltaPressure, nonZeroDeltaPressure );
                    minNonZeroDeltaPressure     = std::min( minNonZeroDeltaPressure, nonZeroDeltaPressure );
                }
            }
        }
    }

    for ( auto it = m_neighborTransmissibilities.begin(); it != m_neighborTransmissibilities.end(); ++it )
    {
        if ( it->first.m_cellIndexSpace == CellAddress::STIMPLAN )
        {
            for ( auto jt = it->second.begin(); jt != it->second.end(); ++jt )
            {
                if ( jt->first.m_cellIndexSpace == CellAddress::ECLIPSE )
                {
                    size_t globalMatrixCellIdx = jt->first.m_globalCellIdx;
                    size_t eclipseResultIndex  = actCellInfo->cellResultIndex( globalMatrixCellIdx );
                    CVF_ASSERT( eclipseResultIndex < currentMatrixPressures.size() );

                    originalLumpedMatrixToFractureTrans[globalMatrixCellIdx] += jt->second;

                    double unsignedDeltaPressure =
                        std::abs( currentMatrixPressures[eclipseResultIndex] - currentWellPressure );
                    double nonZeroDeltaPressure = std::max( epsilonDeltaPressure, unsignedDeltaPressure );

                    jt->second *= nonZeroDeltaPressure / maxNonZeroDeltaPressure;
                }
            }
        }
    }

    if ( minPressureDrop && minNonZeroDeltaPressure != std::numeric_limits<double>::infinity() )
    {
        *minPressureDrop = minNonZeroDeltaPressure;
    }
    if ( maxPressureDrop && maxNonZeroDeltaPressure != std::numeric_limits<double>::infinity() )
    {
        *maxPressureDrop = maxNonZeroDeltaPressure;
    }

    return originalLumpedMatrixToFractureTrans;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::map<size_t, double> RigTransmissibilityCondenser::calculateFicticiousFractureToWellTransmissibilities()
{
    std::map<size_t, double> matrixToAllFracturesTrans;
    for ( auto it = m_neighborTransmissibilities.begin(); it != m_neighborTransmissibilities.end(); ++it )
    {
        if ( it->first.m_cellIndexSpace == CellAddress::STIMPLAN )
        {
            for ( auto jt = it->second.begin(); jt != it->second.end(); ++jt )
            {
                if ( jt->first.m_cellIndexSpace == CellAddress::ECLIPSE )
                {
                    size_t globalMatrixCellIdx = jt->first.m_globalCellIdx;
                    // T'_mf
                    double matrixToFractureTrans = jt->second;
                    // Sum(T'_mf)
                    matrixToAllFracturesTrans[globalMatrixCellIdx] += matrixToFractureTrans;
                }
            }
        }
    }

    std::map<size_t, double> fictitiousFractureToWellTrans; // T'_fjw
    for ( const CellAddress& externalCell : m_externalCellAddrSet )
    {
        if ( externalCell.m_cellIndexSpace == CellAddress::ECLIPSE )
        {
            size_t globalMatrixCellIdx = externalCell.m_globalCellIdx;
            // Sum(T'_mf)
            double scaledMatrixToFractureTrans = matrixToAllFracturesTrans[globalMatrixCellIdx];
            // T'mw
            double scaledMatrixToWellTrans =
                condensedTransmissibility( externalCell, { true, RigTransmissibilityCondenser::CellAddress::WELL, 1 } );
            // T'_fjw
            fictitiousFractureToWellTrans[globalMatrixCellIdx] =
                RigFractureTransmissibilityEquations::effectiveInternalFractureToWellTransPDDHC( scaledMatrixToFractureTrans,
                                                                                                 scaledMatrixToWellTrans );
        }
    }
    return fictitiousFractureToWellTrans;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::map<size_t, double> RigTransmissibilityCondenser::calculateEffectiveMatrixToWellTransmissibilities(
    const std::map<size_t, double>& originalLumpedMatrixToFractureTrans,
    const std::map<size_t, double>& ficticuousFractureToWellTransMap )
{
    std::map<size_t, double> effectiveMatrixToWellTrans;
    for ( const CellAddress& externalCell : m_externalCellAddrSet )
    {
        if ( externalCell.m_cellIndexSpace == CellAddress::ECLIPSE )
        {
            size_t globalMatrixCellIdx = externalCell.m_globalCellIdx;

            auto matrixToFractureIt = originalLumpedMatrixToFractureTrans.find( globalMatrixCellIdx );
            CVF_ASSERT( matrixToFractureIt != originalLumpedMatrixToFractureTrans.end() );
            // Sum(T_mf)
            double lumpedOriginalMatrixToFractureT = matrixToFractureIt->second;
            // T'_fjw
            auto fictitiousFractureToWellIt = ficticuousFractureToWellTransMap.find( globalMatrixCellIdx );
            CVF_ASSERT( fictitiousFractureToWellIt != ficticuousFractureToWellTransMap.end() );
            double fictitiousFractureToWellTrans = fictitiousFractureToWellIt->second;

            // T^dp_mw
            double transmissibilityValue =
                RigFractureTransmissibilityEquations::effectiveMatrixToWellTransPDDHC( lumpedOriginalMatrixToFractureT,
                                                                                       fictitiousFractureToWellTrans );
            if ( transmissibilityValue > m_transmissibilityThreshold )
            {
                effectiveMatrixToWellTrans[globalMatrixCellIdx] = transmissibilityValue;
            }
        }
    }
    return effectiveMatrixToWellTrans;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigTransmissibilityCondenser::calculateCondensedTransmissibilities()
{
    if ( m_neighborTransmissibilities.empty() ) return;

    // Find all equations, and their total ordering

    union
    {
        int idxToFirstExternalEquation;
        int internalEquationCount;
    };
    idxToFirstExternalEquation = -1;
    int totalEquationCount     = -1;

    std::map<CellAddress, int> cellAddressToEqIdxMap;
    std::vector<CellAddress>   eqIdxToCellAddressMapping;
    {
        for ( const auto& adrEqIdxPair : m_neighborTransmissibilities )
        {
            cellAddressToEqIdxMap.insert( { adrEqIdxPair.first, -1 } );
            for ( const auto& adrTranspair : adrEqIdxPair.second )
            {
                cellAddressToEqIdxMap.insert( { adrTranspair.first, -1 } );
            }
        }

        int currentEqIdx = 0;

        for ( auto& adrEqIdxPair : cellAddressToEqIdxMap )
        {
            adrEqIdxPair.second = currentEqIdx;
            eqIdxToCellAddressMapping.push_back( adrEqIdxPair.first );

            if ( idxToFirstExternalEquation == -1 && adrEqIdxPair.first.m_isExternal )
            {
                idxToFirstExternalEquation = currentEqIdx;
            }
            ++currentEqIdx;
        }
        totalEquationCount = currentEqIdx;
    }

    CAF_ASSERT( idxToFirstExternalEquation != -1 );

    using namespace Eigen;

    MatrixXd totalSystem = MatrixXd::Zero( totalEquationCount, totalEquationCount );

    for ( const auto& adrToAdrTransMapPair : m_neighborTransmissibilities )
    {
        CAF_ASSERT( cellAddressToEqIdxMap.count( adrToAdrTransMapPair.first ) ); // Remove when stabilized
        int c1EquationIdx = cellAddressToEqIdxMap[adrToAdrTransMapPair.first];
        for ( const auto& adrTranspair : adrToAdrTransMapPair.second )
        {
            CAF_ASSERT( cellAddressToEqIdxMap.count( adrTranspair.first ) ); // Remove when stabilized

            int c2EquationIdx = cellAddressToEqIdxMap[adrTranspair.first];

            totalSystem( c1EquationIdx, c2EquationIdx ) += adrTranspair.second;
            totalSystem( c2EquationIdx, c1EquationIdx ) += adrTranspair.second;
            totalSystem( c1EquationIdx, c1EquationIdx ) -= adrTranspair.second;
            totalSystem( c2EquationIdx, c2EquationIdx ) -= adrTranspair.second;
        }

        ++c1EquationIdx;
    }

    // std::cout  << "T = " << std::endl <<  totalSystem << std::endl;

    int externalEquationCount = totalEquationCount - internalEquationCount;

    MatrixXd condensedSystem;
    MatrixXd Tee = totalSystem.bottomRightCorner( externalEquationCount, externalEquationCount );

    if ( internalEquationCount == 0 )
    {
        condensedSystem = Tee;
    }
    else
    {
        MatrixXd Tei = totalSystem.bottomLeftCorner( externalEquationCount, internalEquationCount );
        MatrixXd Tii = totalSystem.topLeftCorner( internalEquationCount, internalEquationCount );

        // std::ofstream outFileStream( "D:\\Data\\TestData\\TiiMatrix.txt" );
        // outFileStream << Tii;

        Eigen::FullPivLU<MatrixXd> solver( Tii );

        // outFileStream << std::endl;
        // outFileStream << "Rows x Cols: " << Tii.rows() << "x" << Tii.cols() << std::endl;
        // outFileStream << "Invertible:  " << ( solver.isInvertible() ? "True" : "False" ) << std::endl;
        // outFileStream << "Condition:   " << solver.rcond() << std::endl;
        // outFileStream << "Rank:        " << solver.rank() << std::endl;

        m_TiiInv        = solver.inverse();
        m_Tie           = totalSystem.topRightCorner( internalEquationCount, externalEquationCount );
        condensedSystem = Tee - Tei * m_TiiInv * m_Tie;
    }

    // std::cout  << "Te = " << std::endl <<  condensedSystem << std::endl << std::endl;

    for ( int exEqIdx = 0; exEqIdx < externalEquationCount; ++exEqIdx )
    {
        for ( int exColIdx = exEqIdx + 1; exColIdx < externalEquationCount; ++exColIdx )
        {
            double T = condensedSystem( exEqIdx, exColIdx );
            // if (T != 0.0)
            {
                CellAddress cell1 = eqIdxToCellAddressMapping[exEqIdx + internalEquationCount];
                CellAddress cell2 = eqIdxToCellAddressMapping[exColIdx + internalEquationCount];
                if ( cell1 < cell2 )
                    m_condensedTransmissibilities[cell1][cell2] = T;
                else
                    m_condensedTransmissibilities[cell2][cell1] = T;

                m_externalCellAddrSet.insert( cell1 );
                m_externalCellAddrSet.insert( cell2 );
            }
        }
    }
}

#include "RigFractureCell.h"
#include "RigMainGrid.h"
#include "RimStimPlanFractureTemplate.h"

void printCellAddress( std::stringstream&                        str,
                       const RigMainGrid*                        mainGrid,
                       const RigFractureGrid*                    fractureGrid,
                       RigTransmissibilityCondenser::CellAddress cellAddr )
{
    using CellAddress = RigTransmissibilityCondenser::CellAddress;

    str << ( cellAddr.m_isExternal ? "E " : "I " );

    switch ( cellAddr.m_cellIndexSpace )
    {
        case CellAddress::ECLIPSE:
        {
            if ( cellAddr.m_globalCellIdx > mainGrid->cellCount() )
            {
                str << "ECL - LGR CELL     ";
            }
            else
            {
                str << "ECL ";
                size_t i, j, k;
                mainGrid->ijkFromCellIndex( cellAddr.m_globalCellIdx, &i, &j, &k );
                str << std::setw( 5 ) << i + 1 << std::setw( 5 ) << j + 1 << std::setw( 5 ) << k + 1;
            }
        }
        break;
        case CellAddress::STIMPLAN:
        {
            str << "STP ";
            const RigFractureCell& stpCell = fractureGrid->cellFromIndex( cellAddr.m_globalCellIdx );
            str << std::setw( 5 ) << stpCell.getI() + 1 << std::setw( 5 ) << stpCell.getJ() + 1 << std::setw( 5 ) << " ";
        }
        break;

        case CellAddress::WELL:
        {
            str << "WEL ";
            str << std::setw( 5 ) << cellAddr.m_globalCellIdx << std::setw( 5 ) << " " << std::setw( 5 ) << " ";
        }
        break;
    }

    str << " ";
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::string RigTransmissibilityCondenser::neighborTransDebugOutput( const RigMainGrid*     mainGrid,
                                                                    const RigFractureGrid* fractureGrid )
{
    std::stringstream debugText;
    for ( const auto& adrEqIdxPair : m_neighborTransmissibilities )
    {
        for ( const auto& adrTransPair : adrEqIdxPair.second )
        {
            debugText << "-- ";
            printCellAddress( debugText, mainGrid, fractureGrid, adrEqIdxPair.first );
            printCellAddress( debugText, mainGrid, fractureGrid, adrTransPair.first );

            debugText << " Trans: " << std::setprecision( 10 ) << std::fixed << adrTransPair.second;
            debugText << std::endl;
        }
    }

    return debugText.str();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::string RigTransmissibilityCondenser::condensedTransDebugOutput( const RigMainGrid*     mainGrid,
                                                                     const RigFractureGrid* fractureGrid )
{
    std::stringstream debugText;
    for ( const auto& adrEqIdxPair : m_condensedTransmissibilities )
    {
        for ( const auto& adrTransPair : adrEqIdxPair.second )
        {
            debugText << "-- ";
            printCellAddress( debugText, mainGrid, fractureGrid, adrEqIdxPair.first );
            printCellAddress( debugText, mainGrid, fractureGrid, adrTransPair.first );

            debugText << " Trans: " << std::setprecision( 10 ) << std::fixed << adrTransPair.second;
            debugText << std::endl;
        }
    }

    return debugText.str();
}