//package insilico.core.descriptor.nestle;

import insilico.core.descriptor.Descriptor;
import insilico.core.descriptor.weight.QuantumNumber;
import insilico.core.molecule.InsilicoMolecule;
import org.openscience.cdk.CDKConstants;
import org.openscience.cdk.Molecule;
import org.openscience.cdk.RingSet;
import org.openscience.cdk.interfaces.IAtom;
import org.openscience.cdk.interfaces.IAtomContainer;

/**
 * Custom descriptors and fixes (to be compliant with older Dragon version).
 * 
 * @author Alberto Manganaro (a.manganaro@kode-solutions.net)
 */
public class CustomDescriptors {

    public double nHM;
    public double RBN;
    public double Hy;
    public double N_MlogP_Correction;
    
    
    public CustomDescriptors(InsilicoMolecule mol) {
        
        // Hydrophilic factor
        Hy = Hy(mol);
        
        // no. heavy atoms as in Dragon 5
        nHM = 0;
        
        // no. of rotatable bonds as in Dragon 5
        RBN = 0;
        
        // >n- correction for MLogP
        N_MlogP_Correction = 0;
        
        
        try {
            
            for (IAtom CurAt : mol.GetStructure().atoms()) {
                if (CurAt.getSymbol().equalsIgnoreCase("Si") ||
                    CurAt.getSymbol().equalsIgnoreCase("Cu") ||
                    CurAt.getSymbol().equalsIgnoreCase("Ag") ||
                    CurAt.getSymbol().equalsIgnoreCase("Cd") ||
                    CurAt.getSymbol().equalsIgnoreCase("Sn") ||
                    CurAt.getSymbol().equalsIgnoreCase("Hg") ||
                    CurAt.getSymbol().equalsIgnoreCase("Pb") ||
                    CurAt.getSymbol().equalsIgnoreCase("Bi") ||
                    CurAt.getSymbol().equalsIgnoreCase("As"))
                    nHM++;
            }
                        
        } catch (Exception e) {
            nHM = -999;
        }        
        
        try {
                        
            double[][] ConnAugMatrix = mol.GetMatrixConnectionAugmented();    
            int nSK = ConnAugMatrix.length;
            
            RingSet rings = mol.GetSSSR();
            
            for (int i=0; i<nSK; i++)
                for (int j=(i+1); j<nSK; j++) {
                    // only single bonds
                    if (ConnAugMatrix[i][j] == 1) {
                        // check non terminal
                        if (calculateVD(ConnAugMatrix, i) == 1)
                            continue;
                        if (calculateVD(ConnAugMatrix, j) == 1)
                            continue;
                        
                        // check type of atoms
                        int[] numberL = QuantumNumber.getWeights(mol.GetStructure());
                        if ( (numberL[i] < 2) || (numberL[j] < 2))
                            continue;
                        
                        // check amide
                        if (checkAmide(ConnAugMatrix, i, j))
                            continue;
                        
                        // seems that D5 discard bonds where one atom is connected to a Cl
                        if ( (checkOtherConstraints(ConnAugMatrix, i)) || (checkOtherConstraints(ConnAugMatrix, j)))
                            continue;
                        
                        // check rings
                        boolean inRing = false;
                        for (IAtomContainer curRing : rings.atomContainers())
                            if ( (curRing.contains(mol.GetStructure().getAtom(i))) &&
                                  (curRing.contains(mol.GetStructure().getAtom(j))) ) {
                                inRing = true;
                                break;
                            }
                        if (inRing) continue;
                        
                        RBN++;
                    }
                }
            
            for (int i=0; i<nSK; i++) {
                // N
                if (ConnAugMatrix[i][i] == 7) {
                    if (mol.GetStructure().getAtom(i).getFlag(CDKConstants.ISAROMATIC)) {
                        boolean nn = false;
                        int n_vd = 0;
                        for (int j=0; j<nSK; j++) {
                            if (i==j) continue;
                            if (ConnAugMatrix[i][j] > 0) {
                                n_vd++;
                                if (ConnAugMatrix[j][j] == 7) 
                                    if (mol.GetStructure().getAtom(j).getFlag(CDKConstants.ISAROMATIC))
                                        nn = true;
                            }
                        }
                        if ( (n_vd == 3) && (!nn)) {
                            N_MlogP_Correction++;
                        }
                    }
                }
            }
             
        } catch (Exception e) {
            RBN = -999;
        }
        
    }
    
    
    private int calculateVD(double[][] ConnAugMatrix, int Atom) {
        int nSK = ConnAugMatrix.length;
        int VD = 0;
        for (int i=0; i<nSK; i++) {
            if (i==Atom) continue;
            if (ConnAugMatrix[Atom][i] > 0) 
                VD++;
        }
        return VD;
    }

    private boolean checkAmide(double[][] ConnAugMatrix, int Atom1, int Atom2) {
        int nSK = ConnAugMatrix.length;

        if  ( ( (ConnAugMatrix[Atom1][Atom1] == 6) && (ConnAugMatrix[Atom2][Atom2] == 7) ) ||
              ( (ConnAugMatrix[Atom1][Atom1] == 7) && (ConnAugMatrix[Atom2][Atom2] == 6) ) ){
            int Carbon = ConnAugMatrix[Atom1][Atom1] == 6 ? Atom1 : Atom2;
            for (int i=0; i<nSK; i++) {
                if (i==Carbon) continue;
                if (ConnAugMatrix[Carbon][i] == 2) 
                    if ( (ConnAugMatrix[i][i] == 8) || (ConnAugMatrix[i][i] == 16) )
                        return true;
            }
        }
        
        return false;
    }
    
    
    private boolean checkOtherConstraints(double[][] ConnAugMatrix, int Atom) {
        int nSK = ConnAugMatrix.length;

        int nF=0, nCl=0, nBr=0, nI=0, VD=0;
        
        for (int i=0; i<nSK; i++) {
            if (i==Atom) continue;
            
            if (ConnAugMatrix[Atom][i] > 0) 
                VD++;
            
            // #*
            if ( (ConnAugMatrix[Atom][i] == 3) )
                return true;
            
            // Cl,F,Br,I attach
            if ( (ConnAugMatrix[Atom][i] == 1) && (ConnAugMatrix[i][i] == 17) )
                nCl++;
            if ( (ConnAugMatrix[Atom][i] == 1) && (ConnAugMatrix[i][i] == 9) )
                nF++;
            if ( (ConnAugMatrix[Atom][i] == 1) && (ConnAugMatrix[i][i] == 35) )
                nBr++;
            if ( (ConnAugMatrix[Atom][i] == 1) && (ConnAugMatrix[i][i] == 53) )
                nI++;
        }
        
        if ( (VD-1) == (nCl+nF+nBr+nI)) return true;
        
        return false;
    }    
    
    
    private double Hy(InsilicoMolecule mol) {
        
        try {
            
            double[][] ConnAugMatrix  = mol.GetMatrixConnectionAugmented();

            Molecule CurMol = mol.GetStructure();
            int nSK = CurMol.getAtomCount();

            int Nhy = 0; // hydrophilic groups
            int nC = 0; // number of carbons
            
            for (int At=0; At<nSK; At++) {
                
                int nH = CurMol.getAtom(At).getImplicitHydrogenCount();

                // hydrophilic groups (-OH, -SH, -NH)
                if ( (ConnAugMatrix[At][At] == 7) || (ConnAugMatrix[At][At] == 8) || 
                        (ConnAugMatrix[At][At] == 16) ) {

                    // WRONG COUNT! but consistent with Dragon 7 and 5
                    Nhy += nH;

                    // correct count of hydrophilic groups
                    
//                    int nSng=0, nDbl=0, nTrp=0, nAr=0;
//                    for (int j=0; j<nSK; j++) {
//                        if (j==At) continue;
//                        if (ConnAugMatrix[At][j] > 0) {
//                            if (ConnAugMatrix[At][j] == 1) nSng++;
//                            if (ConnAugMatrix[At][j] == 2) nDbl++;
//                            if (ConnAugMatrix[At][j] == 3) nTrp++;
//                            if (ConnAugMatrix[At][j] == 1.5) nAr++;
//                        }
//                    }  
//                    if ( ((nSng+nDbl+nTrp+nAr)==1) && (nSng==1) && (nH==1))
//                        Nhy++;
                }
                
                // nC
                if ( (ConnAugMatrix[At][At] == 6) ) {
                    nC++;
                }
                
            }
           
            return ( (1 + Nhy) * Log2(1 + Nhy) + nC * ( (1.0/nSK) * Log2(1.0/nSK) ) + Math.sqrt( Nhy / (nSK*nSK)) ) / Log2(1 + nSK);                    
            
        } catch (Throwable e) {
            return Descriptor.MISSING_VALUE;
        }
        
    }
 
    
    private double Log2(double val) {
        return (Math.log(val) / Math.log(2));
    }
    
}