class BCUTDescriptor extends AbstractMolecularDescriptor implements IMolecularDescriptor {

    @Override
    public DescriptorValue calculate(IAtomContainer container) {
        int counter;
        IAtomContainer molecule;
        try {
            molecule = container.clone();
        } catch (CloneNotSupportedException e) {
            logger.debug("Error during clone");
            return getDummyDescriptorValue(new CDKException("Error occurred during clone " + e));
        }

        // add H's in case they're not present
        try {
            AtomContainerManipulator.percieveAtomTypesAndConfigureAtoms(molecule);
            CDKHydrogenAdder hAdder = CDKHydrogenAdder.getInstance(molecule.getBuilder());
            hAdder.addImplicitHydrogens(molecule);
            AtomContainerManipulator.convertImplicitToExplicitHydrogens(molecule);
        } catch (Exception e) {
            return getDummyDescriptorValue(new CDKException("Could not add hydrogens: " + e.getMessage(), e));
        }

        // do aromaticity detecttion for calculating polarizability later on
        if (this.checkAromaticity) {
            try {
                AtomContainerManipulator.percieveAtomTypesAndConfigureAtoms(molecule);
            } catch (CDKException e) {
                return getDummyDescriptorValue(new CDKException("Error in atom typing: " + e.getMessage(), e));
            }
            try {
                Aromaticity.cdkLegacy().apply(molecule);
            } catch (CDKException e) {
                return getDummyDescriptorValue(new CDKException("Error in aromaticity perception: " + e.getMessage()));
            }
        }

        // find number of heavy atoms
        int nheavy = 0;
        for (int i = 0; i < molecule.getAtomCount(); i++) {
            if (molecule.getAtom(i).getAtomicNumber() != IElement.H) nheavy++;
        }

        if (nheavy == 0) return getDummyDescriptorValue(new CDKException("No heavy atoms in the molecule"));

        double[] diagvalue = new double[nheavy];

        // get atomic mass weighted BCUT
        counter = 0;
        try {
            for (int i = 0; i < molecule.getAtomCount(); i++) {
                if (molecule.getAtom(i).getAtomicNumber() == IElement.H) continue;
                diagvalue[counter] = Isotopes.getInstance().getMajorIsotope(molecule.getAtom(i).getSymbol())
                        .getExactMass();
                counter++;
            }
        } catch (Exception e) {
            return getDummyDescriptorValue(new CDKException("Could not calculate weight: " + e.getMessage(), e));
        }

        double[][] burdenMatrix = BurdenMatrix.evalMatrix(molecule, diagvalue);
        if (hasUndefined(burdenMatrix))
            return getDummyDescriptorValue(new CDKException("Burden matrix has undefined values"));
        Matrix matrix = new Matrix(burdenMatrix);
        EigenvalueDecomposition eigenDecomposition = new EigenvalueDecomposition(matrix);
        double[] eval1 = eigenDecomposition.getRealEigenvalues();

        // get charge weighted BCUT
        LonePairElectronChecker lpcheck = new LonePairElectronChecker();
        GasteigerMarsiliPartialCharges peoe;
        try {
            lpcheck.saturate(molecule);
            double[] charges = new double[molecule.getAtomCount()];
            //            pepe = new GasteigerPEPEPartialCharges();
            //            pepe.calculateCharges(molecule);
            //            for (int i = 0; i < molecule.getAtomCount(); i++) charges[i] = molecule.getAtom(i).getCharge();
            peoe = new GasteigerMarsiliPartialCharges();
            peoe.assignGasteigerMarsiliSigmaPartialCharges(molecule, true);
            for (int i = 0; i < molecule.getAtomCount(); i++)
                charges[i] += molecule.getAtom(i).getCharge();
            for (int i = 0; i < molecule.getAtomCount(); i++) {
                molecule.getAtom(i).setCharge(charges[i]);
            }
        } catch (Exception e) {
            return getDummyDescriptorValue(new CDKException("Could not calculate partial charges: " + e.getMessage(), e));
        }
        counter = 0;
        for (int i = 0; i < molecule.getAtomCount(); i++) {
            if (molecule.getAtom(i).getAtomicNumber() == IElement.H) continue;
            diagvalue[counter] = molecule.getAtom(i).getCharge();
            counter++;
        }
        burdenMatrix = BurdenMatrix.evalMatrix(molecule, diagvalue);
        if (hasUndefined(burdenMatrix))
            return getDummyDescriptorValue(new CDKException("Burden matrix has undefined values"));
        matrix = new Matrix(burdenMatrix);
        eigenDecomposition = new EigenvalueDecomposition(matrix);
        double[] eval2 = eigenDecomposition.getRealEigenvalues();

        int[][] topoDistance = PathTools.computeFloydAPSP(AdjacencyMatrix.getMatrix(molecule));

        // get polarizability weighted BCUT
        Polarizability pol = new Polarizability();
        counter = 0;
        for (int i = 0; i < molecule.getAtomCount(); i++) {
            if (molecule.getAtom(i).getAtomicNumber() == IElement.H) continue;
            diagvalue[counter] = pol.calculateGHEffectiveAtomPolarizability(molecule, molecule.getAtom(i), false,
                    topoDistance);
            counter++;
        }
        burdenMatrix = BurdenMatrix.evalMatrix(molecule, diagvalue);
        if (hasUndefined(burdenMatrix))
            return getDummyDescriptorValue(new CDKException("Burden matrix has undefined values"));
        matrix = new Matrix(burdenMatrix);
        eigenDecomposition = new EigenvalueDecomposition(matrix);
        double[] eval3 = eigenDecomposition.getRealEigenvalues();

        // return only the n highest & lowest eigenvalues
        int lnlow, lnhigh, enlow, enhigh;
        if (nlow > nheavy) {
            lnlow = nheavy;
            enlow = nlow - nheavy;
        } else {
            lnlow = nlow;
            enlow = 0;
        }

        if (nhigh > nheavy) {
            lnhigh = nheavy;
            enhigh = nhigh - nheavy;
        } else {
            lnhigh = nhigh;
            enhigh = 0;
        }

        DoubleArrayResult retval = new DoubleArrayResult((lnlow + enlow + lnhigh + enhigh) * 3);

        for (int i = 0; i < lnlow; i++)
            retval.add(eval1[i]);
        for (int i = 0; i < enlow; i++)
            retval.add(Double.NaN);
        for (int i = 0; i < lnhigh; i++)
            retval.add(eval1[eval1.length - i - 1]);
        for (int i = 0; i < enhigh; i++)
            retval.add(Double.NaN);

        for (int i = 0; i < lnlow; i++)
            retval.add(eval2[i]);
        for (int i = 0; i < enlow; i++)
            retval.add(Double.NaN);
        for (int i = 0; i < lnhigh; i++)
            retval.add(eval2[eval2.length - i - 1]);
        for (int i = 0; i < enhigh; i++)
            retval.add(Double.NaN);

        for (int i = 0; i < lnlow; i++)
            retval.add(eval3[i]);
        for (int i = 0; i < enlow; i++)
            retval.add(Double.NaN);
        for (int i = 0; i < lnhigh; i++)
            retval.add(eval3[eval3.length - i - 1]);
        for (int i = 0; i < enhigh; i++)
            retval.add(Double.NaN);

        return new DescriptorValue(getSpecification(), getParameterNames(),
                                   getParameters(), retval,
                                   getDescriptorNames());
    }

    public  BCUTDescriptor();

}