//
// edip_multi.cpp
//
// LGPL Version 2.1 HEADER START
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
//
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library 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
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
// MA 02110-1301  USA
//
// LGPL Version 2.1 HEADER END
//

//
// Copyright (c) 2021, Regents of the University of Minnesota.
// All rights reserved.
//
// Contributors:
//    Yaser Afshar
//

#include "edip_multi.hpp"

#include <cmath>
#include <cstdio>
#include <cstring>
#include <limits>
#include <map>
#include <memory>
#include <string>
#include <vector>

#include "helper.hpp"
#include "special.hpp"
#include "utility.hpp"

namespace edip_multi {
/*!
 * \brief EDIP number of parameters per line.
 *
 * EDIP parameter file. One set of params can span multiple lines.
 * Only store the 1st entry.
 *
 * format of an entry (one or more lines) in a parameter file
 *
 *     element1   element2   element3
 *     A          B          cutoffA   cutoffC   alpha   beta   eta
 *     gamma      lambda     mu        rho       sigma   Q0
 *     u1         u2         u3        u4
 *
 * units for each parameters:
 *
 *     A, lambda -> are in eV
 *     B, cutoffA, cutoffC, gamma, sigma, -> are in Angstrom
 *     alpha, beta, eta, mu, rho, Q0, u1-u4, -> are pure numbers
 *
 */
static constexpr int kParamsPerLine{20};
}  // namespace edip_multi

int EDIPMulti::ProcessParameterFile(
    std::FILE *const parameter_file_pointer, int const max_line_size,
    std::vector<std::string> const &element_name) {
  std::unique_ptr<char> next_line_(new char[max_line_size]);
  char *next_line = next_line_.get();

  // Create a utility object
  Utility ut;

  std::size_t const nelements{element_name.size()};

  std::size_t const maxparam{nelements * nelements * nelements};

  // Resize the local arrays
  params_.resize(maxparam);

  elem3param_.clear();
  elem3param_.resize(nelements, nelements, nelements, maxparam);

  // keep track of known species
  std::map<std::string, std::size_t> species_map;
  for (std::size_t i = 0; i < element_name.size(); ++i) {
    species_map[element_name[i]] = i;
  }

  // Read each set of parameters from the EDIP parameter file.
  // One set of params can span multiple lines. Only store the 1st entry.
  char *words[edip_multi::kParamsPerLine + 1];

  std::size_t nparams{0};

  while (nparams < maxparam) {
    if (ut.GetNextLine(parameter_file_pointer, next_line, max_line_size)) {
      break;
    }

    int nwords = ut.GetNumberOfWordsInLine(next_line);

    // concatenate additional lines until have kParamsPerLine words
    while (nwords < edip_multi::kParamsPerLine) {
      std::size_t const n = std::strlen(next_line);

      if (ut.GetNextLine(parameter_file_pointer, &next_line[n],
                         max_line_size - static_cast<int>(n))) {
        std::string msg{"End of file while reading a line from the `edip` "};
        msg += "parameter file.\nIncorrect format in the `edip` parameter ";
        msg += "file.\n" + std::to_string(nwords) + " parameters are ";
        msg += "provided while the model expects 20 parameters.";
        HELPER_LOG_ERROR(msg);
        return true;
      }

      nwords = ut.GetNumberOfWordsInLine(next_line);
    }  // nwords < kParamsPerLine

    if (nwords != edip_multi::kParamsPerLine) {
      std::string msg{"Incorrect format in the `edip` parameter file.\n"};
      msg += std::to_string(nwords) + " parameters are provided while the ";
      msg += "model expects 20 parameters.";
      HELPER_LOG_ERROR(msg);
      return true;
    }

    nwords = 0;
    words[nwords++] = std::strtok(next_line, "' \t\n\r\f");
    while ((words[nwords] = std::strtok(nullptr, "' \t\n\r\f"))) {
      ++nwords;
    }

    // ielement,jelement,kelement = 1st args
    // if all 3 args are in element list, then parse this line
    // else skip to the next entry in file
    std::string ielement(words[0]);
    std::string jelement(words[1]);
    std::string kelement(words[2]);
    auto iel = species_map.find(ielement);
    if (iel == species_map.end()) {
      continue;
    }
    auto jel = species_map.find(jelement);
    if (jel == species_map.end()) {
      continue;
    }
    auto kel = species_map.find(kelement);
    if (kel == species_map.end()) {
      continue;
    }

    if (elem3param_(iel->second, jel->second, kel->second) != maxparam) {
      HELPER_LOG_ERROR("EDIP parameter file has a duplicate entry.");
      return true;
    }

    elem3param_(iel->second, jel->second, kel->second) = nparams;

    for (int i = 3; i < edip_multi::kParamsPerLine; ++i) {
      if (!ut.IsDouble(words[i])) {
        std::string msg(words[i]);
        msg += " is not a valid floating-point number.";
        HELPER_LOG_ERROR(msg);
        return true;
      }
    }  // Loop over i

    params_[nparams].A = std::atof(words[3]);
    params_[nparams].B = std::atof(words[4]);
    params_[nparams].cutoffA = std::atof(words[5]);
    params_[nparams].cutoffC = std::atof(words[6]);
    params_[nparams].alpha = std::atof(words[7]);
    params_[nparams].beta = std::atof(words[8]);
    params_[nparams].eta = std::atof(words[9]);
    params_[nparams].gamma = std::atof(words[10]);
    params_[nparams].lambda = std::atof(words[11]);
    params_[nparams].mu = std::atof(words[12]);
    params_[nparams].rho = std::atof(words[13]);
    params_[nparams].sigma = std::atof(words[14]);
    params_[nparams].Q0 = std::atof(words[15]);
    params_[nparams].u1 = std::atof(words[16]);
    params_[nparams].u2 = std::atof(words[17]);
    params_[nparams].u3 = std::atof(words[18]);
    params_[nparams].u4 = std::atof(words[19]);

    if (params_[nparams].A < 0.0 || params_[nparams].B < 0.0 ||
        params_[nparams].cutoffA < 0.0 || params_[nparams].cutoffC < 0.0 ||
        params_[nparams].alpha < 0.0 || params_[nparams].beta < 0.0 ||
        params_[nparams].eta < 0.0 || params_[nparams].gamma < 0.0 ||
        params_[nparams].lambda < 0.0 || params_[nparams].mu < 0.0 ||
        params_[nparams].rho < 0.0 || params_[nparams].sigma < 0.0) {
      HELPER_LOG_ERROR("Illegal EDIP parameter.");
      return true;
    }

    ++nparams;
  }  // while loop

  for (std::size_t i = 0; i < nelements; ++i) {
    for (std::size_t j = 0; j < nelements; ++j) {
      for (std::size_t k = 0; k < nelements; ++k) {
        if (elem3param_(i, j, k) == maxparam) {
          std::string msg = "EDIP parameter file is missing an entry for '";
          msg += element_name[i] + " " + element_name[j] + " ";
          msg += element_name[k] + "' triplet.";
          HELPER_LOG_ERROR(msg);
          return true;
        }
      }  // Loop over k
    }    // Loop over j
  }      // Loop over i

  // everything is good
  return false;
}

void EDIPMulti::ConvertUnit(double const convert_length_factor,
                            double const convert_energy_factor) {
  // (Angstroms in metal units).
  if (special::IsNotOne(convert_length_factor)) {
    convert_length_factor_ = convert_length_factor;

    // B, cutoffA, cutoffC, gamma, sigma, -> are in Angstrom
    for (std::size_t i = 0; i < params_.size(); ++i) {
      params_[i].B *= convert_length_factor;
      params_[i].cutoffA *= convert_length_factor;
      params_[i].cutoffC *= convert_length_factor;
      params_[i].gamma *= convert_length_factor;
      params_[i].sigma *= convert_length_factor;

    }  // Loop over i
  }

  //  (eV in the case of metal units).
  if (special::IsNotOne(convert_energy_factor)) {
    convert_energy_factor_ = convert_energy_factor;

    // A, lambda -> are in eV
    for (std::size_t i = 0; i < params_.size(); ++i) {
      params_[i].A *= convert_energy_factor;
      params_[i].lambda *= convert_energy_factor;
    }
  }
}

void EDIPMulti::CompleteSetup(double *max_cutoff) {
  for (std::size_t i = 0; i < params_.size(); ++i) {
    // set cutoff square
    params_[i].cutsq = params_[i].cutoffA * params_[i].cutoffA;
  }  // Loop over i

  cutmax_ = 0.0;
  for (std::size_t i = 0; i < params_.size(); ++i) {
    // set cutmax to max of all params
    if (cutmax_ < params_[i].cutoffA) {
      cutmax_ = params_[i].cutoffA;
    }
  }  // Loop over i

  *max_cutoff = cutmax_;
}