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

#include "LBMUtils.h"
#include <cmath>
#include <limits>
#include <iostream>
#include <fstream>
#include "FileUtils.h"

#include <algorithm> 
#include <cctype>
#include <locale>

   // trim from start (in place)
static inline void ltrim(std::string& s) {
  s.erase(s.begin(), std::find_if(s.begin(), s.end(), [](int ch) {
    return !std::isspace(ch);
  }));
}

// trim from end (in place)
static inline void rtrim(std::string& s) {
  s.erase(std::find_if(s.rbegin(), s.rend(), [](int ch) {
    return !std::isspace(ch);
  }).base(), s.end());
}

// trim from both ends (in place)
static inline void trim(std::string& s) {
  ltrim(s);
  rtrim(s);
}

bool LBM::WriteCSV(std::string const& filename, const LBM::CSV& data)
{
  std::ofstream csv;
  csv.open(filename, std::ofstream::out | std::ofstream::trunc);
  if (!csv.good())
    return false;

  size_t num_rows = 0;
  size_t num_columns = data.size();
  // Write out the headers
  for (auto columns : data)
  {
    num_columns--;
    csv << columns.first;
    if(num_columns > 0)
      csv << ", ";
    num_rows = columns.second.size();
  }
  csv << "\n";

  // Write out the data
  for (size_t idx=0; idx<num_rows; idx++)
  {
    // Write out values
    num_columns = data.size();
    for (auto columns : data)
    {
      num_columns--;
      auto& v = columns.second;
      csv << v[idx];
      if (num_columns > 0)
        csv << ", ";
    }
    csv << "\n";
  }
  csv << "\n";
  csv.close();
  return true;
}

bool LBM::ReadCSV(std::string const& filename, LBM::CSV& data)
{
  data.clear();
  std::ifstream csv;
  csv.open(filename, std::ofstream::in);
  if (!csv.good())
    return false;

  int line_num = 1;
  int idx;
  bool get_names = true;
  std::string line;
  while (!csv.eof()) // To get you all the lines.
  {
    idx = 0;
    std::getline(csv, line); // Get the line.
    if (line.empty())
      break;
    std::stringstream s_stream(line); // Create stringstream from the line
    while (s_stream.good())
    {
      std::string value;
      std::getline(s_stream, value, ','); //get first string delimited by comma
      if (get_names)
      {
        trim(value);
        data.push_back(std::pair<std::string, std::vector<float>>(value, std::vector<float>()));
      }
      else
      {
        try {
          data[idx++].second.push_back((float)std::stod(value));
        } catch(...){ return false;}
      }
    }
    line_num++;
    get_names = false;
  }

  return true;
}

bool LBM::LoadConfiguration(std::string const& config_file, LBM::Run& lbm)
{
  try
  {
    std::string line;
    std::istringstream iss;

    /////////////////
    // CONFIG FILE //
    /////////////////
    std::ifstream cFile(config_file, std::ios::in);

    float value;
    std::string name;
    while (std::getline(cFile, line))
    {
      if (line.empty() || line.find("#") != std::string::npos)// Skip a comment
        continue;
      if (line.find("ITK") != std::string::npos)// Skip a ITK configuration
        continue;
      iss.clear();
      iss.str(line);
      if (!(iss >> value >> name))
      {
        lbm.Error("Unable to parse line from: " + config_file);
        lbm.Error("  "+line);
        return false;
      }
      if (name.compare("atmospheric_pressure") == 0)
      {
        lbm.cfg.atmospheric_pressure = value; continue;
      }
      if (name.compare("density") == 0)
      {
        lbm.cfg.density = value; continue;
      }
      if (name.compare("fluid_temperature") == 0)
      {
        lbm.cfg.fluid_temperature = value; continue;
      }
      if (name.compare("initial_velocity") == 0)
      {
        lbm.cfg.initial_velocity = value; continue;
      }
      if (name.compare("speed_of_sound") == 0)
      {
        lbm.cfg.speed_of_sound = value; continue;
      }
      if (name.compare("thermal_coefficient") == 0)
      {
        lbm.cfg.thermal_coefficient = value; continue;
      }
      if (name.compare("viscosity") == 0)
      {
        lbm.cfg.viscosity = value; continue;
      }
      if (name.compare("wall_temperature") == 0)
      {
        lbm.cfg.wall_temperature = value; continue;
      }
      if (name.compare("unitless_viscosity") == 0)
      {
        lbm.cfg.unitless_viscosity = value; continue;
      }
      if (name.compare("target_flow_rate") == 0)
      {
        lbm.cfg.target_flow_rate = value; continue;
      }
      /////////////////////////
      // Boundary Conditions //
      /////////////////////////
      if (name.compare("inlet_boundary_condition") == 0)
      {
        if (value == 0)
          lbm.cfg.inlet_boundary_condition = BoundaryCondition::Pressure;
        else
          lbm.cfg.inlet_boundary_condition = BoundaryCondition::Flow;
        continue;
      }
      if (name.compare("inlet_pressure") == 0)
      {
        lbm.cfg.inlet_pressure = value; continue;
      }
      if (name.compare("inlet_volumetric_flow_rate") == 0)
      {
        lbm.cfg.inlet_volumetric_flow_rate = value; continue;
      }
      if (name.compare("outlet_boundary_condition") == 0)
      {
        if (value == 0)
          lbm.cfg.outlet_boundary_condition = BoundaryCondition::Pressure;
        else
          lbm.cfg.outlet_boundary_condition = BoundaryCondition::Flow;
        continue;
      }
      if (name.compare("outlet_pressure") == 0)
      {
        lbm.cfg.outlet_pressure = value; continue;
      }
      if (name.compare("outlet_volumetric_flow_rate") == 0)
      {
        lbm.cfg.outlet_volumetric_flow_rate = value; continue;
      }
      //////////////////////
      // Solver Execution //
      //////////////////////
      if (name.compare("max_iterations") == 0)
      {
        lbm.cfg.max_iterations = (size_t)value; continue;
      }
      if (name.compare("tolerance") == 0)
      {
        lbm.cfg.tolerance = value; continue;
      }
      ////////////////////////////
      // Geometry Specification //
      ////////////////////////////
      if (name.compare("geometry_type") == 0)
      {
        if (value == 1)
          lbm.in.geometry_type = GeometryType::Brick;
        else if (value == 2)
          lbm.in.geometry_type = GeometryType::Cylinder;
        else
        {
          lbm.Error("Unsupported geometry type" + lbm.to_string(value));
          return false;
        }
        continue;
      }
      if (name.compare("grid_spacing") == 0)
      {
        lbm.in.grid_spacing = value; continue;
      }
      if (name.compare("brick_x_length") == 0)
      {
        lbm.in.brick_x_length = value; continue;
      }
      if (name.compare("brick_y_length") == 0)
      {
        lbm.in.brick_y_length = value; continue;
      }
      if (name.compare("brick_z_length") == 0)
      {
        lbm.in.brick_z_length = value; continue;
      }
      if (name.compare("cylinder_radius") == 0)
      {
        lbm.in.cylinder_radius = value; continue;
      }
      if (name.compare("cylinder_length") == 0)
      {
        lbm.in.cylinder_length = value; continue;
      }
    }
  }
  catch (std::exception ex)
  {
    lbm.Error(ex.what());
    return false;
  }
  return true;
}

double LBM::PercentTolerance(double expected, double calculated, double epsilon)
{
  // Check for 'invalid' numbers first
  //if (Double.isNaN(expected) || Double.isNaN(calculated) || Double.isInfinite(expected) || Double.isInfinite(calculated))
  //  Log.warn("While finding percent tolerance from values 'expected' = " + expected + " and 'calculated' =" + calculated +
  //    ", invalid values (NaN or Infinity) were found.  Unexpected results may occur.");
  // Handle special cases
  if (expected == 0.0 && calculated == 0.0)
    return 0.0;
  else if (expected == 0.0 || calculated == 0.0)
  {
    if (std::fabs(expected + calculated) < epsilon)
      return 0.0;
    else
    {
      if (expected == 0.0)
        return std::numeric_limits<double>::infinity();
      else if (expected < 0.0)
        return -100.0;
      else
        return 100.0;
    }
  }
  else
    return std::fabs(calculated - expected) / expected * 100.0;
}

double LBM::PercentDifference(double expected, double calculated, double epsilon)
{
  // Check for 'invalid' numbers first
  //if (Double.isNaN(expected) || Double.isNaN(calculated) || Double.isInfinite(expected) || Double.isInfinite(calculated))
  //  Log.warn("While finding percent difference from values 'expected' = " + expected + " and 'calculated' =" + calculated +
  //    ", invalid values (NaN or Infinity) were found.  Unexpected results may occur.");
  // Handle special cases
  if (expected == 0.0 && calculated == 0.0)
    return 0.0;
  else if (expected == 0.0 || calculated == 0.0)
  {
    if (std::fabs(expected + calculated) < epsilon)
      return 0.0;
    else
      return 200.0;
  }
  else
  {
    double difference = (calculated - expected);
    double average = (calculated + expected) / 2.0;

    if (average == 0.0)
      return std::numeric_limits<double>::infinity();

    return std::fabs(difference / average) * 100.0;
  }
}

bool LBM::CompareDataSet(LBM::Run const& baseline, LBM::Run const& computed, bool vtr_data_only)
{
  double diff;
  bool same = true;
  double percrent_difference = 2.0;

  ////////////////////
  // COMPARE INPUTS //
  ////////////////////

  // Compare the Points
  if (baseline.in.dimensions[0] != computed.in.dimensions[0])
  {
    computed.Error("Number of X dimension points does not match");
    same = false;
  }
  if (baseline.in.dimensions[1] != computed.in.dimensions[1])
  {
    computed.Error("Number of Y dimension points does not match");
    same = false;
  }
  if (baseline.in.dimensions[2] != computed.in.dimensions[2])
  {
    computed.Error("Number of Z dimension points does not match");
    same = false;
  }
  if (!vtr_data_only)
  {
    // TODO Percent compare?
    if (baseline.in.grid_spacing != computed.in.grid_spacing)
    {
      computed.Error("Grid spacing does not match");
      same = false;
    }
    if (baseline.in.labels.size() != computed.in.labels.size())
    {
      computed.Error("Input label size does not match");
      same = false;
    }
    else
    {
      // Check input labels
      for (unsigned int i = 0; i < baseline.in.labels.size(); i++)
      {
        if (baseline.in.labels[i] != computed.in.labels[i])
        {
          computed.Error("in.labels["+std::to_string(i)+"] does not match.");
          same = false;
        }
      }
    }
  }

  /////////////////////
  // COMPARE OUTPUTS //
  /////////////////////

  // Cell Counts
  if (baseline.out.num_cells != computed.out.num_cells)
  {
    computed.Error("Number of cells does not match");
    same = false;
  }
  if (baseline.out.num_wall_cells != computed.out.num_wall_cells)
  {
    computed.Error("Number of wall cells does not match");
    same = false;
  }
  if (baseline.out.num_inflow_cells != computed.out.num_inflow_cells)
  {
    computed.Error("Number of inflow cells does not match");
    same = false;
  }
  if (baseline.out.num_outflow_cells != computed.out.num_outflow_cells)
  {
    computed.Error("Number of outflow cells does not match");
    same = false;
  }
  if (baseline.out.num_inactive_cells != computed.out.num_inactive_cells)
  {
    computed.Error("Number of inactive cells does not match");
    same = false;
  }
  if (baseline.out.num_interior_cells != computed.out.num_interior_cells)
  {
    computed.Error("Number of interior cells does not match");
    same = false;
  }
  if (baseline.out.num_unknown_cells != computed.out.num_unknown_cells)
  {
    computed.Error("Number of unknown cells does not match");
    same = false;
  }

  // Scalar Outputs
  if (!vtr_data_only)
  {
    diff = PercentDifference(baseline.out.inlet_velocity,
      computed.out.inlet_velocity);
    if (diff > percrent_difference)
    {
      computed.Error("out.inlet_velocity does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
    diff = PercentDifference(baseline.out.inlet_area,
      computed.out.inlet_area);
    if (diff > percrent_difference)
    {
      computed.Error("out.inlet_area does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
    diff = PercentDifference(baseline.out.outlet_velocity,
      computed.out.outlet_velocity);
    if (diff > percrent_difference)
    {
      computed.Error("out.outlet_velocity does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }

    diff = PercentDifference(baseline.out.outlet_area,
      computed.out.outlet_area);
    if (diff > percrent_difference)
    {
      computed.Error("out.outlet_area does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }

    diff = PercentDifference(baseline.out.average_diameter,
      computed.out.average_diameter);
    if (diff > percrent_difference)
    {
      computed.Error("out.outlet_diameter does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }

    diff = PercentDifference(baseline.out.average_area,
      computed.out.average_area);
    if (diff > percrent_difference)
    {
      computed.Error("out.average_area does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }

    diff = PercentDifference(baseline.out.mass_flow,
      computed.out.mass_flow);
    if (diff > percrent_difference)
    {
      computed.Error("out.mass_flow does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }

    diff = PercentDifference(baseline.out.pressure_drop,
      computed.out.pressure_drop);
    if (diff > percrent_difference)
    {
      computed.Error("out.pressure_drop does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }

    diff = PercentDifference(baseline.out.volumetric_flow_rate,
      computed.out.volumetric_flow_rate);
    if (diff > percrent_difference)
    {
      computed.Error("out.volumetric_flow_rate does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
  }

  // Array Lengths
  if (baseline.out.Q_num_cell_values != computed.out.Q_num_cell_values)
  {
    computed.Error("Number of Q cell valuess does not match");
    same = false;
  }
  if (baseline.out.T_num_cell_values != computed.out.T_num_cell_values)
  {
    computed.Error("Number of T cell values does not match");
    same = false;
  }
  if (baseline.out.num_boundary_values != computed.out.num_boundary_values)
  {
    computed.Error("Number of boundary_values does not match");
    same = false;
  }
  if (baseline.out.nPDF != computed.out.nPDF)
  {
    computed.Error("Number of nPDF does not match");
    same = false;
  }

  size_t num_cells = baseline.out.num_cells;
  // Arrays
  for (size_t i = 0; i < num_cells; i++)
  {
    if (baseline.out.labels[i] != computed.out.labels[i])
    {
      computed.Error("out.labels["+std::to_string(i)+"] does not match");
      same = false;
    }
    diff = PercentDifference(baseline.out.pressure_flows_and_stress[i],
                             computed.out.pressure_flows_and_stress[i]);
    if (diff > percrent_difference)
    {
      computed.Error("out.pressure["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
    diff = PercentDifference(baseline.out.pressure_flows_and_stress[num_cells+i],
                             computed.out.pressure_flows_and_stress[num_cells+i]);
    if (diff > percrent_difference)
    {
      computed.Error("out.X_flow["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
    diff = PercentDifference(baseline.out.pressure_flows_and_stress[num_cells*2+i],
                             computed.out.pressure_flows_and_stress[num_cells*2+i]);
    if (diff > percrent_difference)
    {
      computed.Error("out.Y_flow["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
      
    diff = PercentDifference(baseline.out.pressure_flows_and_stress[num_cells*3+i],
                             computed.out.pressure_flows_and_stress[num_cells*3+i]);
    if (diff > percrent_difference)
    {
      computed.Error("out.Z_flow["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
    diff = PercentDifference(baseline.out.pressure_flows_and_stress[num_cells*4+i],
                             computed.out.pressure_flows_and_stress[num_cells*4+i]);
    if (diff > percrent_difference)
    {
      computed.Error("out.XY_stress["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
    diff = PercentDifference(baseline.out.pressure_flows_and_stress[num_cells*5+i],
                             computed.out.pressure_flows_and_stress[num_cells*5+i]);
    if (diff > percrent_difference)
    {
      computed.Error("out.XZ_stress["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
    diff = PercentDifference(baseline.out.pressure_flows_and_stress[num_cells*6+i],
                             computed.out.pressure_flows_and_stress[num_cells*6+i]);
    if (diff > percrent_difference)
    {
      computed.Error("out.YZ_stress["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
    diff = PercentDifference(baseline.out.temperature[i],
                             computed.out.temperature[i]);
    if (diff > percrent_difference)
    {
      computed.Error("out.temperature["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
      same = false;
    }
  }
  if (!vtr_data_only)
  {
    for (size_t i = 0; i < baseline.out.nPDF; i++)
    {
      diff = PercentDifference(baseline.out.fPDF[i],
                               computed.out.fPDF[i]);
      if (diff > percrent_difference)
      {
        computed.Error("out.fPDF["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
        same = false;
      }
      diff = PercentDifference(baseline.out.fnPDF[i],
                               computed.out.fnPDF[i]);
      if (diff > percrent_difference)
      {
        computed.Error("out.fnPDF["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
        same = false;
      }
    }
    for (size_t i = 0; i < baseline.out.num_boundary_values; i++)
    {
      diff = PercentDifference(baseline.out.boundary_conditions[i],
                               computed.out.boundary_conditions[i]);
      if (diff > percrent_difference)
      {
        computed.Error("out.boundary_conditions["+std::to_string(i)+"] does not match, the error is "+std::to_string(diff)+"%");
        same = false;
      }
    }
  }
  computed.Info("Finished comparison");
  return same;
}
