//
//  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



#include "bestfit.h"
#include "PBFRDKit.h"


using namespace RDKit;
using namespace std;

double PBFRD(ROMol& mol){
    const Conformer *conf;
    int numAtoms = mol.getNumAtoms();
    if(numAtoms<4) return 0;
    
    if(mol.getNumConformers()!=1){
        cerr << "Molecule has more than one conformation: Exiting.\n";
        exit(0);
    }
    else conf = &(mol.getConformer(0));
    if(!conf){
        cerr << "No conformation found for moelcule: Exiting.\n";
        exit(0);
    }
    if(!conf->is3D()) return 0 ;

    double *points, *plane;
    points = new double[numAtoms*3];
    plane = new double[4];    
    for(int i=0; i<numAtoms; ++i){
        const RDGeom::Point3D pos = conf->getAtomPos(i);
        points[i*3]= pos.x;
        points[i*3+1]= pos.y;
        points[i*3+2]= pos.z;

    } 
    
    getBestFitPlane(numAtoms,points,sizeof(double)*3,0,0,plane);
    
    double denom=0.0;
    for(int i=0; i<3; ++i){
        denom += plane[i]*plane[i];
    }
    denom = pow(denom,0.5);
    
    double res=0.0;
    for(int i=0; i<numAtoms; ++i){
        res+= distanceFromAPlane(points[i*3], points[i*3+1], points[i*3+2], plane, denom);
    }
    
    res = res/numAtoms;
    
    delete points;
    delete plane;
    //delete conf;
    
    return res;
    
}

double distanceFromAPlane(double x, double y, double z, double plane[], double denom){
    double numer=0.0;
    numer = abs(x*plane[0]+y*plane[1]+z*plane[2]+plane[3]);

    return numer/denom;
}