Source code

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 */
021
022package org.biojava.nbio.structure.symmetry.core;
023
024import org.biojava.nbio.structure.geometry.CalcPoint;
025import org.biojava.nbio.structure.geometry.UnitQuaternions;
026import org.biojava.nbio.structure.symmetry.utils.PermutationGenerator;
027
028import javax.vecmath.AxisAngle4d;
029import javax.vecmath.Matrix4d;
030import javax.vecmath.Point3d;
031import javax.vecmath.Quat4d;
032import javax.vecmath.Vector3d;
033
034import java.util.ArrayList;
035import java.util.HashSet;
036import java.util.List;
037import java.util.Set;
038
039
040/**
041 *
042 * @author Peter
043 */
044public class SystematicSolver implements QuatSymmetrySolver {
045        private QuatSymmetrySubunits subunits = null;
046        private QuatSymmetryParameters parameters = null;
047
048        private Point3d[] originalCoords = null;
049        private Point3d[] transformedCoords = null;
050        private RotationGroup rotations = new RotationGroup();
051        private Vector3d centroid = new Vector3d();
052        private Matrix4d centroidInverse = new Matrix4d();
053        private Set<List<Integer>> hashCodes = new HashSet<List<Integer>>();
054
055        public SystematicSolver(QuatSymmetrySubunits subunits, QuatSymmetryParameters parameters) {
056                if (subunits.getSubunitCount()== 2) {
057                        throw new IllegalArgumentException("SystematicSolver cannot be applied to subunits with 2 centers");
058                }
059                this.subunits = subunits;
060                this.parameters = parameters;
061        }
062
063        @Override
064        public RotationGroup getSymmetryOperations() {
065                if (rotations.getOrder() == 0) {
066                        solve();
067                        rotations.complete();
068                }
069                return rotations;
070        }
071
072        private void solve() {
073                initialize();
074                int n = subunits.getSubunitCount();
075                PermutationGenerator g = new PermutationGenerator(n);
076
077                // loop over all permutations
078                while (g.hasMore()) {
079                        int[] perm = g.getNext();
080                        List<Integer> permutation = new ArrayList<Integer>(perm.length);
081                        for (int j = 0; j < n; j++) {
082                                permutation.add(perm[j]);
083                        }
084
085                        if (! isValidPermutation(permutation)) {
086                                continue;
087                        }
088
089                        boolean newPermutation = evaluatePermutation(permutation);
090                        if (newPermutation) {
091                                completeRotationGroup();
092                        }
093
094                        if (rotations.getOrder() >= subunits.getSubunitCount()) {
095                                return;
096                        }
097                }
098        }
099
100        /**
101         * Adds translational component to rotation matrix
102         * @param rotTrans
103         * @param rotation
104         * @return
105         */
106        private void combineWithTranslation(Matrix4d rotation) {
107                rotation.setTranslation(centroid);
108                rotation.mul(rotation, centroidInverse);
109        }
110
111        private Rotation createSymmetryOperation(List<Integer> permutation, Matrix4d transformation, AxisAngle4d axisAngle, int fold, QuatSymmetryScores scores) {
112                Rotation s = new Rotation();
113                s.setPermutation(new ArrayList<Integer>(permutation));
114                s.setTransformation(new Matrix4d(transformation));
115                s.setAxisAngle(new AxisAngle4d(axisAngle));
116                s.setFold(fold);
117                s.setScores(scores);
118                return s;
119        }
120
121        private void completeRotationGroup() {
122                PermutationGroup g = new PermutationGroup();
123                for (int i = 0; i < rotations.getOrder(); i++) {
124                        Rotation s = rotations.getRotation(i);
125                        g.addPermutation(s.getPermutation());
126                }
127                g.completeGroup();
128
129//      System.out.println("Completing rotation group from: " +symmetryOperations.getSymmetryOperationCount() + " to " + g.getPermutationCount());
130
131                // the group is complete, nothing to do
132                if (g.getOrder() == rotations.getOrder()) {
133                        return;
134                }
135
136//      System.out.println("complete group: " +  rotations.getOrder() +"/" + g.getOrder());
137                // try to complete the group
138                for (int i = 0; i < g.getOrder(); i++) {
139                        List<Integer> permutation = g.getPermutation(i);
140                        if (isValidPermutation(permutation)) {
141                                  // perform permutation of subunits
142                                evaluatePermutation(permutation);
143                        }
144                }
145        }
146
147        private boolean isValidPermutation(List<Integer> permutation) {
148                if (permutation.size() == 0) {
149                        return false;
150                }
151
152                // if this permutation is a duplicate, return false
153                if (hashCodes.contains(permutation)) {
154                        return false;
155                }
156
157                // check if permutation is pseudosymmetric
158                if (! isAllowedPermuation(permutation)) {
159                        return false;
160                }
161
162                // get fold and make sure there is only one E (fold=1) permutation
163                int fold = PermutationGroup.getOrder(permutation);
164                if (rotations.getOrder() > 1 && fold == 1) {
165                        return false;
166                }
167                if (fold == 0 || subunits.getSubunitCount() % fold != 0) {
168                        return false;
169                }
170
171                // if this permutation is a duplicate, returns false
172                return hashCodes.add(permutation);
173        }
174
175        private boolean isAllowedPermuation(List<Integer> permutation) {
176                List<Integer> seqClusterId = subunits.getClusterIds();
177                for (int i = 0; i < permutation.size(); i++) {
178                        int j = permutation.get(i);
179                        if (seqClusterId.get(i) != seqClusterId.get(j)) {
180                                return false;
181                        }
182                }
183                return true;
184        }
185
186        private boolean evaluatePermutation(List<Integer> permutation) {
187                // permutate subunits
188                for (int j = 0, n = subunits.getSubunitCount(); j < n; j++) {
189                        transformedCoords[j].set(originalCoords[permutation.get(j)]);
190                }
191
192                int fold = PermutationGroup.getOrder(permutation);
193
194                // TODO implement this piece of code using at origin superposition
195                Quat4d quat = UnitQuaternions.relativeOrientation(
196                                originalCoords, transformedCoords);
197                AxisAngle4d axisAngle = new AxisAngle4d();
198                Matrix4d transformation = new Matrix4d();
199
200                transformation.set(quat);
201                axisAngle.set(quat);
202
203                Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
204                if (axis.lengthSquared() < 1.0E-6) {
205                        axisAngle.x = 0;
206                        axisAngle.y = 0;
207                        axisAngle.z = 1;
208                        axisAngle.angle = 0;
209                } else {
210                        axis.normalize();
211                        axisAngle.x = axis.x;
212                        axisAngle.y = axis.y;
213                        axisAngle.z = axis.z;
214                }
215
216                CalcPoint.transform(transformation, transformedCoords);
217
218                double subunitRmsd = CalcPoint.rmsd(transformedCoords, originalCoords);
219
220                if (subunitRmsd <parameters.getRmsdThreshold()) {
221                        // transform to original coordinate system
222                        combineWithTranslation(transformation);
223                        QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation);
224                        if (scores.getRmsd() < 0.0 || scores.getRmsd() > parameters.getRmsdThreshold()) {
225                                return false;
226                        }
227
228                        scores.setRmsdCenters(subunitRmsd);
229                        Rotation symmetryOperation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
230                        rotations.addRotation(symmetryOperation);
231                        return true;
232                }
233                return false;
234        }
235
236        private void initialize() {
237                // translation to centered coordinate system
238                centroid = new Vector3d(subunits.getCentroid());
239
240                // translation back to original coordinate system
241                Vector3d reverse = new Vector3d(centroid);
242                reverse.negate();
243                centroidInverse.set(reverse);
244                // Make sure matrix element m33 is 1.0. An old version vecmath did not set this element.
245                centroidInverse.setElement(3, 3, 1.0);
246
247                List<Point3d> centers = subunits.getCenters();
248                int n = subunits.getSubunitCount();
249
250                originalCoords = new Point3d[n];
251                transformedCoords = new Point3d[n];
252
253                for (int i = 0; i < n; i++) {
254                        originalCoords[i] = centers.get(i);
255                        transformedCoords[i] = new Point3d();
256                }
257        }
258}