001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 * Created on Apr 7, 2010
021 * Author: Andreas Prlic
022 *
023 */
024
025package org.biojava.nbio.structure.align.util;
026
027
028import org.biojava.nbio.structure.Atom;
029import org.biojava.nbio.structure.Calc;
030import org.biojava.nbio.structure.SVDSuperimposer;
031import org.biojava.nbio.structure.StructureException;
032import org.biojava.nbio.structure.align.model.AFPChain;
033import org.biojava.nbio.structure.jama.Matrix;
034import org.slf4j.Logger;
035import org.slf4j.LoggerFactory;
036
037public class AFPChainScorer {
038
039        private static final Logger logger = LoggerFactory.getLogger(AFPChainScorer.class);
040
041
042        public  static double getTMScore(AFPChain align, Atom[] ca1, Atom[] ca2) throws StructureException
043        {
044                if ( align.getNrEQR() == 0)
045                        return -1;
046
047
048                // Create new arrays for the subset of atoms in the alignment.
049                Atom[] ca1aligned = new Atom[align.getOptLength()];
050                Atom[] ca2aligned = new Atom[align.getOptLength()];
051                int pos=0;
052                int[] blockLens = align.getOptLen();
053                int[][][] optAln = align.getOptAln();
054                assert(align.getBlockNum() <= optAln.length);
055
056                for(int block=0;block< align.getBlockNum();block++) {
057
058                        if ( ! ( blockLens[block] <= optAln[block][0].length)) {
059                                logger.warn("AFPChainScorer getTMScore: errors reconstructing alignment block [" + block + "]. Length is " + blockLens[block] + " but should be <=" + optAln[block][0].length);
060                        }
061
062                        for(int i=0;i<blockLens[block];i++) {
063                                int pos1 = optAln[block][0][i];
064                                int pos2 = optAln[block][1][i];
065                                Atom a1 = ca1[pos1];
066                                Atom a2 = (Atom) ca2[pos2].clone();
067
068                                ca1aligned[pos] = a1;
069                                ca2aligned[pos] = a2;
070                                pos++;
071                        }
072                }
073
074                // this can happen when we load an old XML serialization which did not support modern ChemComp representation of modified residues.
075                if ( pos != align.getOptLength()){
076                        logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + align.getOptLength());
077                        // we need to resize the array, because we allocated too many atoms earlier on.
078                        ca1aligned = (Atom[]) resizeArray(ca1aligned, pos);
079                        ca2aligned = (Atom[]) resizeArray(ca2aligned, pos);
080                }
081                //Superimpose
082                SVDSuperimposer svd = new SVDSuperimposer(ca1aligned, ca2aligned);
083                Matrix matrix = svd.getRotation();
084                Atom shift = svd.getTranslation();
085
086                for(Atom a : ca2aligned) {
087                        Calc.rotate(a, matrix);
088                        Calc.shift(a, shift);
089                }
090
091                return SVDSuperimposer.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length);
092
093        }
094
095        /**
096         * Reallocates an array with a new size, and copies the contents
097         * of the old array to the new array.
098         * @param oldArray  the old array, to be reallocated.
099         * @param newSize   the new array size.
100         * @return          A new array with the same contents.
101         */
102        private static Object resizeArray (Object oldArray, int newSize) {
103                int oldSize = java.lang.reflect.Array.getLength(oldArray);
104                @SuppressWarnings("rawtypes")
105                Class elementType = oldArray.getClass().getComponentType();
106                Object newArray = java.lang.reflect.Array.newInstance(
107                                elementType,newSize);
108                int preserveLength = Math.min(oldSize,newSize);
109                if (preserveLength > 0)
110                        System.arraycopy (oldArray,0,newArray,0,preserveLength);
111                return newArray; }
112
113}