whimdescriptor.java

化学图形处理软件

     *
     * @return The parameterNames value
     */
    public String[] getParameterNames() {
        String[] pname = new String[1];
        pname[0] = "type";
        return (pname);
    }


    /**
     * Gets the parameterType attribute of the WHIMDescriptor object.
     *
     * @param name Description of the Parameter
     * @return The parameterType value
     */
    public Object getParameterType(String name) {        
        return ("");
    }


    /**
     * Calculates 11 directional and 6 non-directional WHIM descriptors for.
     * the specified weighting scheme
     *
     * @param container Parameter is the atom container.
     * @return An ArrayList containing the descriptors in the order described above.
     * @throws CDKException if the principal components decomposition fails
     */
    public DescriptorValue calculate(IAtomContainer container) throws CDKException {
        double sum = 0.0;
        Molecule ac;

        try {
            ac = (Molecule) container.clone();
        } catch (CloneNotSupportedException e) {
            throw new CDKException("Error during clone");
        }

        // do aromaticity detecttion for calculating polarizability later on
        //HueckelAromaticityDetector had = new HueckelAromaticityDetector();
        //had.detectAromaticity(ac);

        // get the coordinate matrix
        double[][] cmat = new double[ac.getAtomCount()][3];
        for (int i = 0; i < ac.getAtomCount(); i++) {
            Point3d coords = ac.getAtom(i).getPoint3d();
            if (coords == null) throw new CDKException("Molecules should have 3D coordinates");
            cmat[i][0] = coords.x;
            cmat[i][1] = coords.y;
            cmat[i][2] = coords.z;
        }

        // set up the weight vector
        Hashtable hash = null;
        double[] wt = new double[ac.getAtomCount()];

        if (this.type.equals("unity")) {
            for (int i = 0; i < ac.getAtomCount(); i++) wt[i] = 1.0;
        } else {
            if (this.type.equals("mass")) {
                hash = this.hashatwt;
            } else if (this.type.equals("volume")) {
                hash = this.hashvdw;
            } else if (this.type.equals("eneg")) {
                hash = this.hasheneg;
            } else if (this.type.equals("polar")) {
                hash = this.hashpol;
            }
            for (int i = 0; i < ac.getAtomCount(); i++) {
                String sym = ac.getAtom(i).getSymbol();
                wt[i] = ((Double) hash.get(sym)).doubleValue();
            }
        }

        PCA pcaobject = null;
        try {
            pcaobject = new PCA(cmat, wt);
        } catch (CDKException cdke) {
            logger.debug(cdke);
        }

        // directional WHIM's
        double[] lambda = pcaobject.getEigenvalues();
        double[] gamma = new double[3];
        double[] nu = new double[3];
        double[] eta = new double[3];

        for (int i = 0; i < 3; i++) sum += lambda[i];
        for (int i = 0; i < 3; i++) nu[i] = lambda[i] / sum;

        double[][] scores = pcaobject.getScores();
        for (int i = 0; i < 3; i++) {
            sum = 0.0;
            for (int j = 0; j < ac.getAtomCount(); j++)
                sum += scores[j][i] * scores[j][i] * scores[j][i] * scores[j][i];
            sum = sum / (lambda[i] * lambda[i] * ac.getAtomCount());
            eta[i] = 1.0 / sum;
        }

        // look for symmetric & asymmetric atoms for the gamma descriptor
        for (int i = 0; i < 3; i++) {
            double ns = 0.0;
            double na = 0.0;
            for (int j = 0; j < ac.getAtomCount(); j++) {
                boolean foundmatch = false;
                for (int k = 0; k < ac.getAtomCount(); k++) {
                    if (k == j) continue;
                    if (scores[j][i] == -1 * scores[k][i]) {
                        ns++;
                        foundmatch = true;
                        break;
                    }
                }
                if (!foundmatch) na++;
            }
            double n = (double) ac.getAtomCount();
            gamma[i] = -1.0 * ((ns / n) * Math.log(ns / n) / Math.log(2.0) + (na / n) * Math.log(1.0 / n) / Math.log(2.0));
            gamma[i] = 1.0 / (1.0 + gamma[i]);
            //logger.debug("ns = "+ns+" na = "+na+"  gamma = "+gamma[i]);
        }

        // non directional WHIMS's
        double t = lambda[0] + lambda[1] + lambda[2];
        double a = lambda[0] * lambda[1] + lambda[0] * lambda[2] + lambda[1] * lambda[2];
        double v = t + a + lambda[0] * lambda[1] * lambda[2];

        double k = 0.0;
        sum = 0.0;
        for (int i = 0; i < 3; i++) sum += lambda[i];
        for (int i = 0; i < 3; i++) k = (lambda[i] / sum) - (1.0 / 3.0);
        k = k / (4.0 / 3.0);

        double g = Math.pow(gamma[0] * gamma[1] * gamma[2], 1.0 / 3.0);
        double d = eta[0] + eta[1] + eta[2];

        // return all the stuff we calculated
        DoubleArrayResult retval = new DoubleArrayResult(11 + 6);
        retval.add(lambda[0]);
        retval.add(lambda[1]);
        retval.add(lambda[2]);

        retval.add(nu[0]);
        retval.add(nu[1]);

        retval.add(gamma[0]);
        retval.add(gamma[1]);
        retval.add(gamma[2]);

        retval.add(eta[0]);
        retval.add(eta[1]);
        retval.add(eta[2]);

        retval.add(t);
        retval.add(a);
        retval.add(v);
        retval.add(k);
        retval.add(g);
        retval.add(d);

        String[] names = {
                "Wlambda1", "Wlambda2", "Wlambda3",
                "Wnu1", "Wnu2",
                "Wgamma1", "Wgamma2", "Wgamma3",
                "Weta1", "Weta2", "Weta3",
                "WT", "WA", "WV", "WK", "WG", "WD"
        };
        for (int i = 0; i < names.length; i++) names[i] += "." + type;
        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval, names);
    }

    /**
     * Returns the specific type of the DescriptorResult object.
     * <p/>
     * The return value from this method really indicates what type of result will
     * be obtained from the {@link org.openscience.cdk.qsar.DescriptorValue} object. Note that the same result
     * can be achieved by interrogating the {@link org.openscience.cdk.qsar.DescriptorValue} object; this method
     * allows you to do the same thing, without actually calculating the descriptor.
     *
     * @return an object that implements the {@link org.openscience.cdk.qsar.result.IDescriptorResult} interface indicating
     *         the actual type of values returned by the descriptor in the {@link org.openscience.cdk.qsar.DescriptorValue} object
     */
    public IDescriptorResult getDescriptorResultType() {
        return new DoubleArrayResult();
    }


    class PCA {

        Matrix evec;
        Matrix t;
        double[] eval;


        public PCA(double[][] cmat, double[] wt) throws CDKException {

            int ncol = 3;
            int nrow = wt.length;

            if (cmat.length != wt.length) {
                throw new CDKException("WHIMDescriptor: number of weights should be equal to number of atoms");
            }

            // make a copy of the coordinate matrix
            double[][] d = new double[nrow][ncol];
            for (int i = 0; i < nrow; i++) {
                for (int j = 0; j < ncol; j++)
                    d[i][j] = cmat[i][j];
            }

            // do mean centering - though the first paper used
            // barymetric centering
            for (int i = 0; i < ncol; i++) {
                double mean = 0.0;
                for (int j = 0; j < nrow; j++)
                    mean += d[j][i];
                mean = mean / (double) nrow;
                for (int j = 0; j < nrow; j++)
                    d[j][i] = (d[j][i] - mean);
            }

            // get the covariance matrix
            double[][] covmat = new double[ncol][ncol];
            double sumwt = 0.;
            for (int i = 0; i < nrow; i++) sumwt += wt[i];
            for (int i = 0; i < ncol; i++) {
                double meanx = 0;
                for (int k = 0; k < nrow; k++) meanx += d[k][i];
                meanx = meanx / (double) nrow;
                for (int j = 0; j < ncol; j++) {
                    double meany = 0.0;
                    for (int k = 0; k < nrow; k++) meany += d[k][j];
                    meany = meany / (double) nrow;

                    double sum = 0.;
                    for (int k = 0; k < nrow; k++) {
                        //double dd =  wt[k] * (d[k][i] - meanx) * (d[k][j] - meany);
                        //logger.debug("("+i+","+j+") "+wts[k] + " * " + d[k][i] + "-" + meanx + " * " + d[k][j] + "-" + meany + "==" + dd);
                        sum += wt[k] * (d[k][i] - meanx) * (d[k][j] - meany);
                    }
                    //logger.debug(sum+" / "+sumwt+"="+sum/sumwt);
                    covmat[i][j] = sum / sumwt;
                }
            }

            // get the loadings (ie eigenvector matrix)
            Matrix m = new Matrix(covmat);
            EigenvalueDecomposition ed = m.eig();
            this.eval = ed.getRealEigenvalues();
            this.evec = ed.getV();
            Matrix x = new Matrix(d);
            this.t = x.times(this.evec);
        }

        double[] getEigenvalues() {
            return (this.eval);
        }

        double[][] getScores() {
            return (t.getArray());
        }
    }

}