//
//  Copyright (c) 2016, Guillaume GODIN
//  All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//     * Redistributions in binary form must reproduce the above
//       copyright notice, this list of conditions and the following
//       disclaimer in the documentation and/or other materials provided
//       with the distribution.
//     * Neither the name of Institue of Cancer Research.
//       nor the names of its contributors may be used to endorse or promote
//       products derived from this software without specific prior written
//       permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Adding RBF descriptors to 3D descriptors by Guillaume Godin

#include <GraphMol/RDKitBase.h>

#include "RDF.h"
#include "MolData3Ddescriptors.h"

#include <math.h>

namespace RDKit {
namespace Descriptors {
namespace {
MolData3Ddescriptors moldata3D;

std::vector<double> getG(int n) {
  std::vector<double> res(n);
  for (int i = 0; i < n; i++) {
    res[i] = 1 + i * 0.5;
  }
  return res;
}

std::vector<double> prepareIState(const ROMol& mol) {
  std::vector<double> IState = moldata3D.GetIStateDrag(
      mol);  // get the real IState value not the EState!
  return IState;
}

void getRDFDesc(double* DM, const ROMol& mol, const Conformer& conf,
                std::vector<double>& res) {
  // std::vector<double> reserror(210, 0); need to avoid the res NaNs
  // if (numAtoms < 4) return reserror;
  // if (!conf.is3D()) return reserror;

  int numAtoms = conf.getNumAtoms();

  std::vector<double> R = getG(30);
  std::vector<double> R1(30);
  std::vector<double> R2(30);
  std::vector<double> R3(30);
  std::vector<double> R4(30);
  std::vector<double> R5(30);
  std::vector<double> R6(30);
  std::vector<double> R7(30);

  std::vector<double> Mass = moldata3D.GetRelativeMW(mol);
  std::vector<double> RelativePol = moldata3D.GetRelativePol(mol);
  std::vector<double> IonPol = moldata3D.GetRelativeIonPol(mol);
  std::vector<double> RelativeElectroNeg = moldata3D.GetRelativeENeg(mol);
  std::vector<double> RelativeVdW = moldata3D.GetRelativeVdW(mol);
  std::vector<double> IState = prepareIState(mol);

  double p;
  for (size_t i = 0; i < R.size(); i++) {
    double res1 = 0.0;
    double res2 = 0.0;
    double res3 = 0.0;
    double res4 = 0.0;
    double res5 = 0.0;
    double res6 = 0.0;
    double res7 = 0.0;

    for (int j = 0; j < numAtoms - 1; j++) {
      for (int k = j + 1; k < numAtoms; k++) {
        p = exp(-100 * pow(R[i] - DM[j * numAtoms + k], 2));
        res1 += p;                                                  // "u"
        res2 += Mass[j] * Mass[k] * p;                              // "m"
        res3 += RelativeVdW[j] * RelativeVdW[k] * p;                // "v"
        res4 += RelativeElectroNeg[j] * RelativeElectroNeg[k] * p;  //"e"
        res5 += RelativePol[j] * RelativePol[k] * p;                // "p"
        res6 += IonPol[j] * IonPol[k] * p;                          // "i"
        res7 += IState[j] * IState[k] * p;                          // "s"
      }
    }
    R1[i] = round(1000 * res1) / 1000;
    R2[i] = round(1000 * res2) / 1000;
    R3[i] = round(1000 * res3) / 1000;
    R4[i] = round(1000 * res4) / 1000;
    R5[i] = round(1000 * res5) / 1000;
    R6[i] = round(1000 * res6) / 1000;
    R7[i] = round(1000 * res7) / 1000;
  }

  R1.insert(R1.end(), R2.begin(), R2.end());
  R1.insert(R1.end(), R3.begin(), R3.end());
  R1.insert(R1.end(), R4.begin(), R4.end());
  R1.insert(R1.end(), R5.begin(), R5.end());
  R1.insert(R1.end(), R6.begin(), R6.end());
  R1.insert(R1.end(), R7.begin(), R7.end());

  res = R1;
}

void getRDFDescCustom(double* DM, const ROMol& mol, const Conformer& conf,
                      std::vector<double>& res,
                      const std::string& customAtomPropName) {
  int numAtoms = conf.getNumAtoms();

  std::vector<double> R = getG(30);
  std::vector<double> R1(30);
  std::vector<double> customAtomArray =
      moldata3D.GetCustomAtomProp(mol, customAtomPropName);

  double p;
  for (size_t i = 0; i < R.size(); i++) {
    double res = 0.0;
    for (int j = 0; j < numAtoms - 1; j++) {
      for (int k = j + 1; k < numAtoms; k++) {
        p = exp(-100 * pow(R[i] - DM[j * numAtoms + k], 2));
        res += customAtomArray[j] * customAtomArray[k] * p;  // "custom"
      }
    }
    R1[i] = round(1000 * res) / 1000;
  }
  res = R1;
}

void GetRDF(double* dist3D, const ROMol& mol, const Conformer& conf,
            std::vector<double>& res) {
  getRDFDesc(dist3D, mol, conf, res);
}

void GetRDFone(double* dist3D, const ROMol& mol, const Conformer& conf,
               std::vector<double>& res, const std::string customAtomPropName) {
  getRDFDescCustom(dist3D, mol, conf, res, customAtomPropName);
}

}  // end of anonymous namespace

void RDF(const ROMol& mol, std::vector<double>& res, int confId,
         const std::string& customAtomPropName) {
  // RDF010u RDF015u RDF020u RDF025u RDF030u RDF035u RDF040u RDF045u RDF050u
  // RDF055u RDF060u RDF065u RDF070u RDF075u RDF080u RDF085u RDF090u RDF095u
  // RDF100u RDF105u RDF110u RDF115u RDF120u RDF125u RDF130u RDF135u RDF140u
  // RDF145u RDF150u RDF155u
  // RDF010m RDF015m RDF020m RDF025m RDF030m RDF035m RDF040m RDF045m RDF050m
  // RDF055m RDF060m RDF065m RDF070m RDF075m RDF080m RDF085m RDF090m RDF095m
  // RDF100m RDF105m RDF110m RDF115m RDF120m RDF125m RDF130m RDF135m RDF140m
  // RDF145m RDF150m RDF155m
  // RDF010v RDF015v RDF020v RDF025v RDF030v RDF035v RDF040v RDF045v RDF050v
  // RDF055v RDF060v RDF065v RDF070v RDF075v RDF080v RDF085v RDF090v RDF095v
  // RDF100v RDF105v RDF110v RDF115v RDF120v RDF125v RDF130v RDF135v RDF140v
  // RDF145v RDF150v RDF155v
  // RDF010e RDF015e RDF020e RDF025e RDF030e RDF035e RDF040e RDF045e RDF050e
  // RDF055e RDF060e RDF065e RDF070e RDF075e RDF080e RDF085e RDF090e RDF095e
  // RDF100e RDF105e RDF110e RDF115e RDF120e RDF125e RDF130e RDF135e RDF140e
  // RDF145e RDF150e RDF155e
  // RDF010p RDF015p RDF020p RDF025p RDF030p RDF035p RDF040p RDF045p RDF050p
  // RDF055p RDF060p RDF065p RDF070p RDF075p RDF080p RDF085p RDF090p RDF095p
  // RDF100p RDF105p RDF110p RDF115p RDF120p RDF125p RDF130p RDF135p RDF140p
  // RDF145p RDF150p RDF155p
  // RDF010i RDF015i RDF020i RDF025i RDF030i RDF035i RDF040i RDF045i RDF050i
  // RDF055i RDF060i RDF065i RDF070i RDF075i RDF080i RDF085i RDF090i RDF095i
  // RDF100i RDF105i RDF110i RDF115i RDF120i RDF125i RDF130i RDF135i RDF140i
  // RDF145i RDF150i RDF155i
  // RDF010s RDF015s RDF020s RDF025s RDF030s RDF035s RDF040s RDF045s RDF050s
  // RDF055s RDF060s RDF065s RDF070s RDF075s RDF080s RDF085s RDF090s RDF095s
  // RDF100s RDF105s RDF110s RDF115s RDF120s RDF125s RDF130s RDF135s RDF140s
  // RDF145s RDF150s RDF155s

  PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers")
  // int numAtoms = mol.getNumAtoms();
  // if (numAtoms < 4) return reserror;

  const Conformer& conf = mol.getConformer(confId);
  // if (!conf.is3D()) return reserror;

  double* dist3D =
      MolOps::get3DDistanceMat(mol, confId, false, true);  // 3D distance matrix

  if (customAtomPropName != "") {
    // std::cout << " Using CustomAtomPropertie\n";
    // do something
    res.clear();
    res.resize(30);  // 7 * 30
    GetRDFone(dist3D, mol, conf, res, customAtomPropName);
  } else {
    res.clear();
    res.resize(210);  // 7 * 30
    GetRDF(dist3D, mol, conf, res);
  }
}
}  // namespace Descriptors
}  // namespace RDKit