//
//  Copyright (C) 2020 Brian P. Kelley
//
//   @@ All Rights Reserved @@
//  This file is part of the RDKit.
//  The contents are covered by the terms of the BSD license
//  which is included in the file license.txt, found at the root
//  of the RDKit source tree.
//
#ifdef RDK_HAS_EIGEN3
#include "BCUT.h"
#include "Crippen.h"
#include <Eigen/Dense>
#include <GraphMol/RDKitBase.h>
#include "GraphMol/PartialCharges/GasteigerCharges.h"
#include "GraphMol/PartialCharges/GasteigerParams.h"
#include <RDGeneral/types.h>

namespace RDKit {
namespace Descriptors {
  // diagonal elements are a property (atomic num, charge, etc)
  // off diagonal are 1/sqrt(bond_order)
  //  Original burden matrix was .1, .2, .3, .15 for single,double,triple or aromatic
  //  all other elements are .001
namespace {
std::unique_ptr<Eigen::MatrixXd>  make_burden(const ROMol &m) {
  auto num_atoms = m.getNumAtoms();
  std::unique_ptr<Eigen::MatrixXd> burden(
	 new Eigen::MatrixXd(num_atoms, num_atoms));
    
  for(unsigned int i=0;i<num_atoms; ++i) {
    for(unsigned int j=0;j<num_atoms; ++j) {
      (*burden)(i,j) = (*burden)(j,i) = 0.001;
    }
  }
  
  for(auto &bond : m.bonds()) {
    unsigned int i = bond->getBeginAtomIdx();
    unsigned int j = bond->getEndAtomIdx();
    double score = 0.0;
    switch(bond->getBondType()) {
    case Bond::AROMATIC:
      // score = 0.15; orig burden
      score = 0.8164965809277261; // 1/sqrt(1.5)
      break;
    case Bond::SINGLE:
      // score = 0.1;
      score = 1.0; // 1/sqrt(1.0)
      break;
    case Bond::DOUBLE:
      // score = 0.2;
      score = 0.7071067811865475; // 1/sqrt(2.0)
      break;
    case Bond::TRIPLE:
      // score = 0.3;	
      score = 0.5773502691896258; // 1/sqrt(3);
      break;
    default:
      CHECK_INVARIANT(0, "Bond order must be Single, Double, Triple or Aromatic");
    }
    (*burden)(i,j) = (*burden)(j,i) = score;
  }    
  return burden;
}
  
std::pair<double,double> BCUT2D(std::unique_ptr<Eigen::MatrixXd> &burden,
				const std::vector<double> &atom_props) {
  for(unsigned int i=0; i<atom_props.size(); ++i) {
    (*burden)(i,i) = atom_props[i];
  }
  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(*burden);
  auto eivals = es.eigenvalues();
  double lowest = eivals(0);
  double highest = eivals(atom_props.size()-1);
  return std::pair<double,double>(highest,lowest);
}
}
    
std::pair<double,double> BCUT2D(const ROMol &m, const std::vector<double> &atom_props) {
  unsigned int num_atoms = m.getNumAtoms();
  PRECONDITION(atom_props.size() == num_atoms, "Number of atom props not equal to number of atoms");
  
  if (num_atoms == 0) {
    return std::pair<double,double>(0,0);
  }
  auto burden = make_burden(m);
  return BCUT2D(burden, atom_props);
}
  
std::pair<double,double> BCUT2D(const ROMol &m, const std::string &atom_double_prop) {
  std::vector<double> props;
  props.reserve(m.getNumAtoms());
  for(auto &atom : m.atoms()) {
    props.push_back(atom->getProp<double>(atom_double_prop));
  }
  return BCUT2D(m, props);
}

std::vector<double> BCUT2D(const ROMol &m) {
  std::unique_ptr<ROMol> mol(MolOps::removeHs(m));
  std::vector<double> masses;
  std::vector<double> charges;
  unsigned int num_atoms = mol->getNumAtoms();
  masses.reserve(num_atoms);
  charges.reserve(num_atoms);

  RDKit::computeGasteigerCharges(*mol, 12, true);
  for(auto &atom: mol->atoms()) {
    masses.push_back(atom->getMass());
    charges.push_back(atom->getProp<double>(common_properties::_GasteigerCharge));
  }
  
  std::vector<double> slogp(num_atoms, 0.0);
  std::vector<double> cmr(num_atoms, 0.0);
  getCrippenAtomContribs(*mol, slogp, cmr);
  
  // polarizability? - need model
  // slogp?  sasa?
  auto burden = make_burden(m);
  auto atom_bcut = BCUT2D(burden, masses);
  auto gasteiger = BCUT2D(burden, charges);
  auto logp = BCUT2D(burden, slogp);
  auto mr = BCUT2D(burden, cmr);
  std::vector<double> res = {atom_bcut.first, atom_bcut.second,
			     gasteiger.first, gasteiger.second,
			     logp.first, logp.second,
			     mr.first, mr.second
  };
  return res;
}
}
}

#endif