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 */
021package org.biojava.nbio.structure.symmetry.internal;
022
023import java.util.ArrayList;
024import java.util.Collections;
025import java.util.List;
026import java.util.Set;
027
028import javax.vecmath.GMatrix;
029import javax.vecmath.GVector;
030
031import org.biojava.nbio.structure.Atom;
032import org.biojava.nbio.structure.StructureException;
033import org.biojava.nbio.structure.align.model.AFPChain;
034import org.biojava.nbio.structure.align.multiple.MultipleAlignment;
035import org.biojava.nbio.structure.align.util.AlignmentTools;
036import org.biojava.nbio.structure.symmetry.utils.SymmetryTools;
037import org.jgrapht.Graph;
038import org.jgrapht.alg.connectivity.ConnectivityInspector;
039import org.jgrapht.graph.DefaultEdge;
040
041/**
042 * The GraphRefiner transforms the self-alignment into a Graph and extracts its
043 * maximally connected Components. It then refines the alignment by combining
044 * the compatible Components with the following heuristic:
045 *
046 * <pre>
047 * Given a set of components and their pairwise compatibilities, iteratively
048 * add the most compatible component, which is compatible to all the components
049 * already added, to the final alignment.
050 * </pre>
051 *
052 * @author Aleix Lafita
053 *
054 */
055public class GraphComponentRefiner implements SymmetryRefiner {
056
057        @Override
058        public MultipleAlignment refine(AFPChain selfAlignment, Atom[] atoms, int order)
059                        throws StructureException, RefinerFailedException {
060
061                // Construct the alignment graph with jgrapht
062                Graph<Integer, DefaultEdge> graph = SymmetryTools
063                                .buildSymmetryGraph(selfAlignment);
064
065                // Find the maximally connected components of the graph
066                ConnectivityInspector<Integer, DefaultEdge> inspector = new ConnectivityInspector<Integer, DefaultEdge>(
067                                graph);
068                List<Set<Integer>> components = inspector.connectedSets();
069
070                // Filter components with size != order, and transform to ResidueGroups
071                List<ResidueGroup> groups = new ArrayList<ResidueGroup>();
072                for (Set<Integer> comp : components) {
073                        if (comp.size() == order) {
074                                ResidueGroup group = new ResidueGroup(comp);
075                                groups.add(group);
076                        }
077                }
078                int size = groups.size();
079                if (size == 0)
080                        throw new RefinerFailedException("Could not find any components "
081                                        + "in the alignment Graph of size != order.");
082
083                // Create a square matrix of component compatibility
084                GMatrix matrix = new GMatrix(size, size);
085                for (int i = 0; i < size; i++) {
086                        for (int j = i; j < size; j++) {
087                                // The diagonal is always 0
088                                if (i == j){
089                                        matrix.setElement(i, j, 0);
090                                        continue;
091                                }
092                                // If compatible put 1, otherwise 0
093                                ResidueGroup g1 = groups.get(i);
094                                ResidueGroup g2 = groups.get(j);
095                                if (g1.isCompatible(g2)) {
096                                        matrix.setElement(i, j, 1);
097                                        matrix.setElement(j, i, 1);
098                                } else {
099                                        matrix.setElement(i, j, 0);
100                                        matrix.setElement(j, i, 0);
101                                }
102                        }
103                }
104
105                // The compatibility score is the sum of rows of the matrix
106                List<Integer> rowScores = new ArrayList<Integer>(size);
107                for (int i = 0; i < size; i++) {
108                        GVector row = new GVector(size);
109                        matrix.getRow(i, row);
110                        // because element={0,1}, this is the sum
111                        int rowScore = (int) row.normSquared();
112                        rowScores.add(rowScore);
113                }
114
115                // Refined multiple alignment Block as a result
116                List<List<Integer>> alignRes = new ArrayList<List<Integer>>(order);
117                for (int i = 0; i < order; i++)
118                        alignRes.add(new ArrayList<Integer>());
119
120                // Iterate until no more groups left to add (all groups score 0)
121                while (true) {
122                        // Take the most compatible ResidueGroup
123                        Integer max = Collections.max(rowScores);
124                        int index = rowScores.indexOf(max);
125
126                        // Add the group to the alignment Block
127                        groups.get(index).combineWith(alignRes);
128
129                        // Zero all the scores of incompatible groups
130                        boolean allZero = true;
131                        for (int i=0; i<size; i++){
132                                if (matrix.getElement(index, i) < 1.0)
133                                        rowScores.set(i, 0);
134                                else if (rowScores.get(i) != 0)
135                                        allZero = false;
136                        }
137                        if (allZero)
138                                break;
139                }
140
141                for (int i = 0; i < order; i++)
142                        Collections.sort(alignRes.get(i));
143
144                int length = alignRes.get(0).size();
145                if (length == 0)
146                        throw new RefinerFailedException("Empty alignment");
147
148                int[][][] optAln = new int[order][2][length];
149                for (int bk = 0; bk < order; bk++) {
150                        optAln[bk] = new int[2][];
151                        optAln[bk][0] = new int[length];
152                        optAln[bk][1] = new int[length];
153                        for (int pos = 0; pos < length; pos++) {
154                                optAln[bk][0][pos] = alignRes.get(bk).get(pos);
155                                optAln[bk][1][pos] = alignRes.get((bk + 1) % order).get(pos);
156                        }
157                }
158                AFPChain afp = AlignmentTools.replaceOptAln(optAln, selfAlignment, atoms, atoms);
159                return SymmetryTools.fromAFP(afp, atoms);
160        }
161
162}