/*
 * Copyright 2009-2020 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

// Standard includes
#include <locale>

// Third party includes
#include <boost/format.hpp>

// VOTCA includes
#include <votca/tools/constants.h>

// Local VOTCA includes
#include "votca/xtp/gnode.h"
#include "votca/xtp/kmccalculator.h"
#include "votca/xtp/logger.h"
#include "votca/xtp/qmstate.h"
#include "votca/xtp/rate_engine.h"
#include "votca/xtp/topology.h"

using namespace std;

namespace votca {
namespace xtp {

void KMCCalculator::ParseCommonOptions(const tools::Property& options) {
  _seed = options.get(".seed").as<Index>();

  _numberofcarriers = options.get(".numberofcarriers").as<Index>();
  _injection_name = options.get(".injectionpattern").as<std::string>();
  _maxrealtime = options.get(".maxrealtime").as<double>();
  _trajectoryfile = options.get(".trajectoryfile").as<std::string>();
  _temperature = options.get(".temperature").as<double>();

  _temperature *= (tools::conv::kB * tools::conv::ev2hrt);
  _occfile = options.get(".occfile").as<std::string>();
  _ratefile = options.get(".ratefile").as<std::string>();

  _injectionmethod = options.get(".injectionmethod").as<std::string>();
}

void KMCCalculator::LoadGraph(Topology& top) {

  std::vector<Segment>& segs = top.Segments();

  if (segs.size() < 1) {
    throw std::runtime_error("Your state file contains no segments!");
  }
  _nodes.reserve(segs.size());
  for (Segment& seg : segs) {
    bool injectable = false;
    if (tools::wildcmp(_injection_name, seg.getType())) {
      injectable = true;
    }
    _nodes.push_back(GNode(seg, _carriertype, injectable));
  }

  QMNBList& nblist = top.NBList();
  if (nblist.size() < 1) {
    throw std::runtime_error("neighborlist contains no pairs!");
  }
  if (_temperature <= 0) {
    throw std::runtime_error(
        "Your Temperature is negative or zero, please specify the temperature "
        "in Kelvin.");
  }

  Rate_Engine rate_engine(_temperature, _field);
  XTP_LOG(Log::error, _log) << "\nCalculating initial rates." << std::flush;
  XTP_LOG(Log::error, _log) << rate_engine << std::flush;
  XTP_LOG(Log::error, _log)
      << "    carriertype: " << _carriertype.ToLongString() << std::flush;

  for (const QMPair* pair : nblist) {
    Rate_Engine::PairRates rates = rate_engine.Rate(*pair, _carriertype);
    _nodes[pair->Seg1()->getId()].AddEventfromQmPair(*pair, _nodes,
                                                     rates.rate12);
    _nodes[pair->Seg2()->getId()].AddEventfromQmPair(*pair, _nodes,
                                                     rates.rate21);
  }
  _RandomVariable.setMaxInt(Index(_nodes.size()));
  XTP_LOG(Log::error, _log) << "    Rates for " << _nodes.size()
                            << " sites are computed." << std::flush;
  WriteRatestoFile(_ratefile, nblist);

  Index events = 0;
  Index max = std::numeric_limits<Index>::min();
  Index min = std::numeric_limits<Index>::max();
  double minlength = std::numeric_limits<double>::max();
  double maxlength = 0;
  for (const auto& node : _nodes) {

    Index size = Index(node.Events().size());
    for (const auto& event : node.Events()) {
      if (event.isDecayEvent()) {
        continue;
      }
      double dist = event.getDeltaR().norm();
      if (dist > maxlength) {
        maxlength = dist;
      } else if (dist < minlength) {
        minlength = dist;
      }
    }

    events += size;
    if (size == 0) {
      XTP_LOG(Log::info, _log)
          << "Node " << node.getId() << " has 0 jumps" << std::flush;
    } else if (size < min) {
      min = size;
    } else if (size > max) {
      max = size;
    }
  }
  double avg = double(events) / double(_nodes.size());
  double deviation = 0.0;
  for (const auto& node : _nodes) {
    double size = double(node.Events().size());
    deviation += (size - avg) * (size - avg);
  }
  deviation = std::sqrt(deviation / double(_nodes.size()));

  XTP_LOG(Log::error, _log)
      << "Nblist has " << nblist.size() << " pairs. Nodes contain " << events
      << " jump events\n"
      << "with avg=" << avg << " std=" << deviation << " max=" << max
      << " min=" << min << " jumps per site\n"
      << "Minimum jumpdistance =" << minlength * tools::conv::bohr2nm
      << " nm Maximum distance =" << maxlength * tools::conv::bohr2nm << " nm\n"
      << std::flush;
  double conv = std::pow(tools::conv::bohr2nm, 3);
  XTP_LOG(Log::error, _log)
      << "spatial carrier density: "
      << double(_numberofcarriers) / (top.BoxVolume() * conv) << " nm^-3"
      << std::flush;

  for (auto& node : _nodes) {
    node.InitEscapeRate();
    node.MakeHuffTree();
  }
  return;
}

void KMCCalculator::ResetForbiddenlist(
    std::vector<GNode*>& forbiddenlist) const {
  forbiddenlist.clear();
  return;
}

void KMCCalculator::AddtoForbiddenlist(
    GNode& node, std::vector<GNode*>& forbiddenlist) const {
  forbiddenlist.push_back(&node);
  return;
}

bool KMCCalculator::CheckForbidden(
    const GNode& node, const std::vector<GNode*>& forbiddenlist) const {
  bool forbidden = false;
  for (const GNode* fnode : forbiddenlist) {
    if (&node == fnode) {
      forbidden = true;
      break;
    }
  }
  return forbidden;
}

bool KMCCalculator::CheckSurrounded(
    const GNode& node, const std::vector<GNode*>& forbiddendests) const {
  bool surrounded = true;
  for (const auto& event : node.Events()) {
    bool thisevent_possible = true;
    for (const GNode* fnode : forbiddendests) {
      if (event.getDestination() == fnode) {
        thisevent_possible = false;
        break;
      }
    }
    if (thisevent_possible == true) {
      surrounded = false;
      break;
    }
  }
  return surrounded;
}

void KMCCalculator::RandomlyCreateCharges() {

  XTP_LOG(Log::error, _log) << "looking for injectable nodes..." << std::flush;
  for (Index i = 0; i < _numberofcarriers; i++) {
    Chargecarrier newCharge(i);
    RandomlyAssignCarriertoSite(newCharge);

    XTP_LOG(Log::error, _log)
        << "starting position for charge " << i << ": segment "
        << newCharge.getCurrentNodeId() << std::flush;
    _carriers.push_back(newCharge);
  }
  return;
}

void KMCCalculator::RandomlyAssignCarriertoSite(Chargecarrier& Charge) {
  Index nodeId_guess = -1;
  do {
    nodeId_guess = _RandomVariable.rand_uniform_int();
  } while (_nodes[nodeId_guess].isOccupied() ||
           _nodes[nodeId_guess].isInjectable() ==
               false);  // maybe already occupied? or maybe not injectable?
  if (Charge.hasNode()) {
    Charge.ReleaseNode();
  }
  Charge.settoNote(&_nodes[nodeId_guess]);

  return;
}

double KMCCalculator::Promotetime(double cumulated_rate) {
  double dt = 0;
  double rand_u = 1 - _RandomVariable.rand_uniform();
  dt = -1 / cumulated_rate * std::log(rand_u);
  return dt;
}

const GLink& KMCCalculator::ChooseHoppingDest(const GNode& node) {
  double u = 1 - _RandomVariable.rand_uniform();
  return *(node.findHoppingDestination(u));
}

Chargecarrier* KMCCalculator::ChooseAffectedCarrier(double cumulated_rate) {
  if (_carriers.size() == 1) {
    return &_carriers[0];
  }
  Chargecarrier* carrier = nullptr;
  double u = 1 - _RandomVariable.rand_uniform();
  for (Index i = 0; i < _numberofcarriers; i++) {
    u -= _carriers[i].getCurrentEscapeRate() / cumulated_rate;
    if (u <= 0 || i == _numberofcarriers - 1) {
      carrier = &_carriers[i];
      break;
    }
  }
  return carrier;
}
void KMCCalculator::WriteRatestoFile(std::string filename,
                                     const QMNBList& nblist) {
  XTP_LOG(Log::error, _log)
      << "\nRates are written to " << filename << std::flush;
  fstream ratefs;
  ratefs.open(filename, fstream::out);
  ratefs << "#PairID,SiteID1,SiteID2, ,rate12[1/s],rate21[1/s] at "
         << _temperature * tools::conv::hrt2ev / tools::conv::kB
         << "K for carrier:" << _carriertype.ToString() << endl;

  Rate_Engine rate_engine(_temperature, _field);
  for (const QMPair* pair : nblist) {
    Rate_Engine::PairRates rates = rate_engine.Rate(*pair, _carriertype);
    ratefs << pair->getId() << " " << pair->Seg1()->getId() << " "
           << pair->Seg2()->getId() << " " << rates.rate12 << " "
           << rates.rate21 << "\n";
  }
  ratefs << std::flush;
  ratefs.close();
}

void KMCCalculator::WriteOccupationtoFile(double simtime,
                                          std::string filename) {
  XTP_LOG(Log::error, _log)
      << "\nOccupations are written to " << filename << std::flush;
  fstream probs;
  probs.open(filename, fstream::out);
  probs << "#SiteID, Occupation prob at "
        << _temperature * tools::conv::hrt2ev / tools::conv::kB
        << "K for carrier:" << _carriertype.ToString() << endl;
  for (const GNode& node : _nodes) {
    double occupationprobability = node.OccupationTime() / simtime;
    probs << node.getId() << "\t" << occupationprobability << endl;
  }
  probs.close();
}

}  // namespace xtp
}  // namespace votca