//
//  Copyright (c) 2012, Institue of Cancer Research.
//  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.
//
// For more information on the Plane of Best Fit please see http://pubs.acs.org/doi/abs/10.1021/ci300293f
//
//  If this code has been useful to you, please include the reference
//  in any work which has made use of it:

//  Plane of Best Fit: A Novel Method to Characterize the Three-Dimensionality of Molecules, Nicholas C. Firth, Nathan Brown, and Julian Blagg, Journal of Chemical Information and Modeling 2012 52 (10), 2516-2525

//
//
// Created by Nicholas Firth, November 2011
// Modified by Greg Landrum for inclusion in the RDKit distribution November 2012
//

#include "PBFRDKit.h"
#include <Numerics/Matrix.h>
#include <Numerics/SquareMatrix.h>
#include <Numerics/SymmMatrix.h>
#include <boost/foreach.hpp>

#include <Eigen/Dense>
using namespace RDKit;

void getSmallestEigenVector(double fSumXX,double fSumXY,double fSumXZ,
                            double fSumYY,double fSumYZ,double fSumZZ,
                            double &x,double &y, double &z);

double distanceFromAPlane(const RDGeom::Point3D &pt,const std::vector<double> &plane, double denom){
  double numer=0.0;
  numer = std::abs(pt.x*plane[0]+pt.y*plane[1]+pt.z*plane[2]+plane[3]);

  return numer/denom;
}

bool getBestFitPlane(const std::vector<RDGeom::Point3D> &points,
                     std::vector<double> &plane,
                     const std::vector<double> *weights) {
  PRECONDITION((!weights || weights->size()>=points.size()),"bad weights vector");
  RDGeom::Point3D origin(0,0,0);
  double wSum=0.0;

  for(unsigned int i=0;i<points.size();++i){
    if(weights){
      double w=(*weights)[i];
      wSum+=w;
      origin+=points[i]*w;
    } else {
      wSum+=1;
      origin+=points[i];
    }
  }
  origin /= wSum;

  double sumXX=0,sumXY=0,sumXZ=0,sumYY=0,sumYZ=0,sumZZ=0;
  for(unsigned int i=0;i<points.size();++i){
    RDGeom::Point3D delta=points[i]-origin;
    if(weights){
      double w=(*weights)[i];
      delta *= w;
    }
    sumXX += delta.x*delta.x;
    sumXY += delta.x*delta.y;
    sumXZ += delta.x*delta.z;
    sumYY += delta.y*delta.y;
    sumYZ += delta.y*delta.z;
    sumZZ += delta.z*delta.z;
  }
  sumXX/=wSum;
  sumXY/=wSum;
  sumXZ/=wSum;
  sumYY/=wSum;
  sumYZ/=wSum;
  sumZZ/=wSum;

  Eigen::Matrix3d mat;
  mat << sumXX, sumXY, sumXZ,
    sumXY, sumYY, sumYZ,
    sumXZ, sumYZ, sumZZ;
  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(mat);
  if(eigensolver.info()!=Eigen::Success){
    BOOST_LOG(rdErrorLog)<<"eigenvalue calculation did not converge"<<std::endl;
    return 0.0;
  }
  RDGeom::Point3D normal;
  normal.x=eigensolver.eigenvectors()(0,0);
  normal.y=eigensolver.eigenvectors()(1,0);
  normal.z=eigensolver.eigenvectors()(2,0);

  plane[0] = normal.x;
  plane[1] = normal.y;
  plane[2] = normal.z;
  plane[3] = -1*normal.dotProduct(origin);
  
}

double PBFRD(ROMol& mol,int confId){
  PRECONDITION(mol.getNumConformers()>=1,"molecule has no conformers")
  int numAtoms = mol.getNumAtoms();
  if(numAtoms<4) return 0;

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

  std::vector<RDGeom::Point3D> points;
  points.reserve(numAtoms);
  for(unsigned int i=0; i<numAtoms; ++i){
    points.push_back(conf.getAtomPos(i));
  } 
    
  std::vector<double> plane(4);
  getBestFitPlane(points,plane,0);

  double denom=0.0;
  for(unsigned int i=0; i<3; ++i){
    denom += plane[i]*plane[i];
  }
  denom = pow(denom,0.5);
    
  double res=0.0;
  for(unsigned int i=0; i<numAtoms; ++i){
    res+= distanceFromAPlane(points[i], plane, denom);
  }
  res /= numAtoms;

  return res;
}