/* Distributed under the Apache License, Version 2.0.
   See accompanying NOTICE file for details.*/

#include "stdafx.h"
#include "substance/SESubstanceTransport.h"
#include "Eigen/Core"
#include "properties/SEScalarVolume.h"

//#define VERBOSE
#define ZERO_APPROX 1e-10

template <typename GraphType, typename FluxUnit, typename QuantityUnit, typename ExtensiveUnit, typename IntensiveUnit>
SESubstanceTransporter<GraphType, FluxUnit, QuantityUnit, ExtensiveUnit, IntensiveUnit>::SESubstanceTransporter(const FluxUnit& fUnit,
                                                                                                                const QuantityUnit& qUnit, 
                                                                                                                const ExtensiveUnit& eUnit, 
                                                                                                                const IntensiveUnit& iUnit, 
                                                                                                                Logger* logger) : Loggable(logger),
                                                                                                                m_FluxUnit(fUnit),
                                                                                                                m_QuantityUnit(qUnit),
                                                                                                                m_ExtensiveUnit(eUnit),
                                                                                                                m_IntensiveUnit(iUnit)
{

}


//--------------------------------------------------------------------------------------------------
/// \brief
/// Calculates the intensive and extensive properties for this Graph
///
/// \param  timestep time-step of the system
///
/// \details
/// Calculates the intensive (liquid concentration or gas volume fraction) and extensive (mass in liquid or volume of gas)
/// properties for all vertexs in the graph.
/// This uses an Eigen solver.
/// An INTENSIVE property is a bulk property, meaning that it is a physical property of a system that does not depend on the system size or the amount of material in the system.
/// For us, this includes Concentration and VolumeFraction
/// An EXTENSIVE property is additive for independent, non-interacting subsystems. The property is proportional to the amount of material in the system.
/// For us, this includes Mass and Volume
//--------------------------------------------------------------------------------------------------
template <typename GraphType, typename FluxUnit, typename QuantityUnit, typename ExtensiveUnit, typename IntensiveUnit>
void SESubstanceTransporter<GraphType, FluxUnit, QuantityUnit, ExtensiveUnit, IntensiveUnit>::Transport(GraphType& graph, double timeStep_s)
{
  //We can assume the flows and volumes are already solved for this time-step, but no substances moved

  auto& verticies = graph.GetVerticies();
  size_t numVerts = verticies.size();

  // We assume every vertex has the same amount of substances and are in the same order
  size_t numSubstances = verticies[0]->GetTransportSubstances().size();
  if (numSubstances == 0)
    return;

  // The A matrix is a square matrix for all nodes
  //All zeros - we'll modify the values later
  //Ax=b
  Eigen::MatrixXd AMatrix = Eigen::MatrixXd::Zero(numVerts, numVerts);
  Eigen::MatrixXd AMatrixInverse = Eigen::MatrixXd::Zero(numVerts, numVerts);

  //1) Determine A matrix - this will be the same for all substances -------------------------------------------------
  for (size_t i = 0; i < numVerts; i++)
  {
    auto* v = verticies[i];
    auto& quantity = v->GetQuantity();
    //Put the current quantity in the matrix
    //If it doesn't have a quantity, just leave it zero
    if (quantity.IsInfinity())
    {
      //Infinite quantity means intensive property is constant
      AMatrix(i, i) = 1.0;
      continue;
    }
    else if (!quantity.IsValid() || quantity.IsZero())
    {
      //If there's no quantity and therefore no intensive or extensive properties, we need to mimic a super small volume
      //Setting it to zero and adding extensive balance from paths in and out will not solve
      //Also setting it to a constant value of zero and no adding the in and out path will not solve
      //Basically divide by zero / un-invertable matrix issues
      //The solution is to use an approximately zero quantity
      AMatrix(i, i) += 1e-20;// ZERO_APPROX;
    }
    else
    {
      if (quantity.IsNegative())
      {
        ///\error Fatal: Can't transport with a negative quantity included.
        std::stringstream ss;
        ss << "Can't transport with a negative volume included. Node = " << v->GetName()
          << ". Volume = " << v->GetQuantity().GetValue(VolumeUnit::mL) << " mL";
        ///\error Fatal: Can't transport with a negative volume included.
        Fatal(ss);
      }

      AMatrix(i, i) += quantity.GetValue(m_QuantityUnit);
    }

    auto* sEdges = graph.GetSourceEdges(*v);
    if (sEdges != nullptr)
    {
      for (auto* srcEdge : *sEdges)
      {
        if (srcEdge->HasFlux())
        {
          double flux = srcEdge->GetFlux().GetValue(m_FluxUnit);

          if (std::abs(flux) < ZERO_APPROX)
          {
            //This is basically zero and not worth dealing with
            continue;
          }

          if (flux < 0.0)
          {
            //Flow in
            AMatrix(i, graph.GetVertexIndex(srcEdge->GetTargetVertex())) -= std::abs(flux) * timeStep_s;
          }
          else
          {
            //Flow out
            AMatrix(i, i) += std::abs(flux) * timeStep_s;
          }
        }
      }
    }

    auto* tEdges = graph.GetTargetEdges(*v);
    if (tEdges != nullptr)
    {
      for (auto* tgtEdge : *tEdges)
      {
        if (tgtEdge->HasFlux())
        {
          double flux = tgtEdge->GetFlux().GetValue(m_FluxUnit);

          if (std::abs(flux) < ZERO_APPROX)
          {
            //This is basically zero and not worth dealing with
            continue;
          }

          if (flux < 0.0)
          {
            //Flow out
            AMatrix(i, i) += std::abs(flux) * timeStep_s;
          }
          else
          {
            //Flow in          
            AMatrix(i, graph.GetVertexIndex(tgtEdge->GetSourceVertex())) -= std::abs(flux) * timeStep_s;
          }
        }
      }
    }
  }

  //Since A doesn't change, we'll just calculate it's inverse and reuse it to solve for x for each substance
  AMatrixInverse = AMatrix.inverse();

  //2) Loop over the extensive/intensive properties per substances -------------------------------------------------  
  for (size_t s=0; s<numSubstances; s++)
  {
    //Ax = b
    Eigen::VectorXd xVector(numVerts); //x
    Eigen::VectorXd bVector(numVerts); //b
    bVector = Eigen::VectorXd::Zero(numVerts);

    //2a) Put the extensive property in the b vector for each node -------------------------------------------------
    for (size_t i = 0; i<numVerts; i++)
    {
      auto* v = verticies[i];
      if (!v->HasQuantity())
        continue;
      auto* q = v->GetTransportSubstances()[s];
      if (q->HasIntensive())
      {
        if (v->GetQuantity().IsInfinity())
        {
          //Concentration doesn't change - this happens at reference nodes a lot (e.g. the Environment)
          bVector(i) = q->GetIntensive().GetValue(m_IntensiveUnit);
        }
        else
        {
          bVector(i) += q->GetExtensive().GetValue(m_ExtensiveUnit);
        }
      }
    }

    //2b) Solve for the intensive property x vector -------------------------------------------------
  xVector = AMatrixInverse * bVector;

    //2c) Parse the intensive property and calculate the new intensive property -------------------------------------------------
    for (size_t i = 0; i<numVerts; i++)
    {
      auto* v = verticies[i];
      auto* q = v->GetTransportSubstances()[s];
      if (q->HasExtensive())
      {
        if (!v->HasQuantity())
        {
          q->Invalidate();
          continue;
        }
        double newIntensiveVlaue = xVector(i);
        if (std::abs(newIntensiveVlaue) < ZERO_APPROX)
        {
          newIntensiveVlaue = 0.0;
        }
        //We force this value becuase it can get above 1 when there is a disconinuity that causes numerical error.
        //That generally only occurs when an extreme action or condition is applied.
        //The balance call after this loop will adjust everything to sum nicely to 1.0.
        q->GetIntensive().ForceValue(newIntensiveVlaue, m_IntensiveUnit);        
      }
    }
  }
  graph.BalanceByIntensive();
}
