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