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 java.util.ArrayList;
025import java.util.HashMap;
026import java.util.HashSet;
027import java.util.List;
028import java.util.Map;
029import java.util.Set;
030
031import javax.vecmath.AxisAngle4d;
032import javax.vecmath.Matrix4d;
033import javax.vecmath.Point3d;
034import javax.vecmath.Quat4d;
035import javax.vecmath.Vector3d;
036
037import org.biojava.nbio.structure.geometry.CalcPoint;
038import org.biojava.nbio.structure.geometry.MomentsOfInertia;
039import org.biojava.nbio.structure.geometry.UnitQuaternions;
040import org.biojava.nbio.structure.symmetry.geometry.DistanceBox;
041import org.biojava.nbio.structure.symmetry.geometry.SphereSampler;
042
043
044/**
045 *
046 * @author Peter
047 */
048public class RotationSolver implements QuatSymmetrySolver {
049        private QuatSymmetrySubunits subunits = null;
050        private QuatSymmetryParameters parameters = null;
051
052        private double distanceThreshold = 0.0f;
053        private DistanceBox<Integer> box = null;
054        private Vector3d centroid = new Vector3d();
055        private Matrix4d centroidInverse = new Matrix4d();
056        private Point3d[] originalCoords = null;
057        private Point3d[] transformedCoords = null;
058        // Cache whether a permutation is invalid (null) vs has been added to rotations
059        private Map<List<Integer>,Rotation> evaluatedPermutations = new HashMap<>();
060
061        private RotationGroup rotations = new RotationGroup();
062
063        public RotationSolver(QuatSymmetrySubunits subunits, QuatSymmetryParameters parameters) {
064                if (subunits.getSubunitCount()== 2) {
065                        throw new IllegalArgumentException("RotationSolver cannot be applied to subunits with 2 centers");
066                }
067                this.subunits = subunits;
068                this.parameters = parameters;
069        }
070
071        @Override
072        public RotationGroup getSymmetryOperations() {
073                if (rotations.getOrder() == 0) {
074                        solve();
075                        completeRotationGroup();
076                        rotations.complete();
077                }
078
079                return rotations;
080        }
081
082        private void solve() {
083                initialize();
084
085                int maxSymOps = subunits.getSubunitCount();
086
087                boolean isSpherical = isSpherical();
088
089                // for cases with icosahedral symmetry n cannot be higher than 60, should check for spherical symmetry here
090                // isSpherical check added 08-04-11
091                if (maxSymOps % 60 == 0 && isSpherical) {
092                        maxSymOps = 60;
093                 }
094
095                AxisAngle4d sphereAngle = new AxisAngle4d();
096                Matrix4d transformation = new Matrix4d();
097
098                int n = subunits.getSubunitCount();
099
100                int sphereCount = SphereSampler.getSphereCount();
101
102                List<Double> angles = getAngles();
103
104                for (int i = 0; i < sphereCount; i++) {
105                        // Sampled orientation
106                        //TODO The SphereSampler samples 4D orientation space. We really
107                        // only need to sample 3D unit vectors, since we use limited
108                        // angles. -SB
109                        SphereSampler.getAxisAngle(i, sphereAngle);
110
111                        // Each valid rotation angle
112                        for (double angle : angles) {
113                                // apply rotation
114                                sphereAngle.angle = angle;
115                                transformation.set(sphereAngle);
116                                // Make sure matrix element m33 is 1.0. It's not on Linux.
117                                transformation.setElement(3, 3, 1.0);
118                                for (int j = 0; j < n; j++) {
119                                        transformedCoords[j].set(originalCoords[j]);
120                                        transformation.transform(transformedCoords[j]);
121                                }
122
123                                // get permutation of subunits and check validity/uniqueness
124                                List<Integer> permutation = getPermutation();
125        //              System.out.println("Rotation Solver: permutation: " + i + ": " + permutation);
126
127                                // check if novel
128                                if ( evaluatedPermutations.containsKey(permutation)) {
129                                        continue; //either invalid or already added
130                                }
131
132                                Rotation newPermutation = isValidPermutation(permutation);
133                                if (newPermutation != null) {
134                                        completeRotationGroup(newPermutation);
135                                }
136
137                                // check if all symmetry operations have been found.
138                                if (rotations.getOrder() >= maxSymOps) {
139                                        return;
140                                }
141                        }
142                }
143        }
144
145        /**
146         * Combine current rotations to make all possible permutations.
147         * If these are all valid, add them to the rotations
148         * @param additionalRots Additional rotations we are considering adding to this.rotations
149         * @return whether the rotations were valid and added
150         */
151        private boolean completeRotationGroup(Rotation... additionalRots) {
152                PermutationGroup g = new PermutationGroup();
153                for (Rotation s : rotations) {
154                        g.addPermutation(s.getPermutation());
155                }
156                for( Rotation s : additionalRots) {
157                        g.addPermutation(s.getPermutation());
158                        // inputs should not have been added already
159                        assert evaluatedPermutations.get(s.getPermutation()) == null;
160                }
161                g.completeGroup();
162
163                // the group is complete, nothing to do
164                if (g.getOrder() == rotations.getOrder()+additionalRots.length) {
165                        for( Rotation s : additionalRots) {
166                                addRotation(s);
167                        }
168                        return true;
169                }
170
171                // try to complete the group
172                List<Rotation> newRots = new ArrayList<>(g.getOrder());
173                // First, quick check for whether they're allowed
174                for (List<Integer> permutation : g) {
175                        if( evaluatedPermutations.containsKey(permutation)) {
176                                Rotation rot = evaluatedPermutations.get(permutation);
177                                if( rot == null ) {
178                                        return false;
179                                }
180                        } else {
181                                if( ! isAllowedPermutation(permutation)) {
182                                        return false;
183                                }
184                        }
185                }
186                // Slower check including the superpositions
187                for (List<Integer> permutation : g) {
188                        Rotation rot;
189                        if( evaluatedPermutations.containsKey(permutation)) {
190                                rot = evaluatedPermutations.get(permutation);
191                        } else {
192                                rot = isValidPermutation(permutation);
193                        }
194
195                        if( rot == null ) {
196                                // if any induced rotation is invalid, abort
197                                return false;
198                        }
199                        if(!evaluatedPermutations.containsKey(permutation)){ //novel
200                                newRots.add(rot);
201                        }
202                }
203                // Add rotations
204                for( Rotation rot : newRots) {
205                        addRotation(rot);
206                }
207                return true;
208        }
209
210        private void addRotation(Rotation rot) {
211                evaluatedPermutations.put(rot.getPermutation(),rot);
212                rotations.addRotation(rot);
213        }
214
215        /**
216         * Superimpose subunits based on the given permutation. Then check whether
217         * the superposition passes RMSD thresholds and create a Rotation to
218         * represent it if so.
219         * @param permutation A list specifying which subunits should be aligned by the current transformation
220         * @return A Rotation representing the permutation, or null if the superposition did not meet thresholds.
221         */
222        private Rotation superimposePermutation(List<Integer> permutation) {
223                // permutate subunits
224                for (int j = 0, n = subunits.getSubunitCount(); j < n; j++) {
225                        transformedCoords[j].set(originalCoords[permutation.get(j)]);
226                }
227
228                int fold = PermutationGroup.getOrder(permutation);
229
230                // get optimal transformation and axisangle by subunit superposition
231                // TODO implement this piece of code using at origin superposition
232                Quat4d quat = UnitQuaternions.relativeOrientation(
233                                originalCoords, transformedCoords);
234                AxisAngle4d axisAngle = new AxisAngle4d();
235                Matrix4d transformation = new Matrix4d();
236                
237                transformation.set(quat);
238                axisAngle.set(quat);
239                
240                Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
241                if (axis.lengthSquared() < 1.0E-6) {
242                        axisAngle.x = 0;
243                        axisAngle.y = 0;
244                        axisAngle.z = 1;
245                        axisAngle.angle = 0;
246                } else {
247                        axis.normalize();
248                        axisAngle.x = axis.x;
249                        axisAngle.y = axis.y;
250                        axisAngle.z = axis.z;
251                }
252                
253                CalcPoint.transform(transformation, transformedCoords);
254                double subunitRmsd = CalcPoint.rmsd(transformedCoords, originalCoords);
255
256                if (subunitRmsd < parameters.getRmsdThreshold()) {
257                        combineWithTranslation(transformation);
258
259                        // evaluate superposition of CA traces
260                        QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation);
261                        if (scores.getRmsd() < 0.0 || scores.getRmsd() > parameters.getRmsdThreshold()) {
262                                return null;
263                        }
264
265                        scores.setRmsdCenters(subunitRmsd);
266                        Rotation symmetryOperation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
267                        return symmetryOperation;
268                }
269                return null;
270        }
271
272        /**
273         * Get valid rotation angles given the number of subunits
274         * @return The rotation angle corresponding to each fold of {@link Subunits#getFolds()}
275         */
276        private List<Double> getAngles() {
277                int n = subunits.getSubunitCount();
278                // for spherical symmetric cases, n cannot be higher than 60
279                if (n % 60 == 0 && isSpherical()) {
280                         n = 60;
281                }
282                List<Integer> folds = subunits.getFolds();
283                List<Double> angles = new ArrayList<Double>(folds.size()-1);
284
285                // note this loop starts at 1, we do ignore 1-fold symmetry, which is the first entry
286                for (int fold: folds) {
287                        if (fold > 0 && fold <= n) {
288                                angles.add(2* Math.PI/fold);
289                        }
290                }
291                return angles;
292        }
293
294        private boolean isSpherical() {
295                MomentsOfInertia m = subunits.getMomentsOfInertia();
296                return m.getSymmetryClass(0.05) == MomentsOfInertia.SymmetryClass.SYMMETRIC;
297        }
298
299        /**
300         * Checks if a particular permutation is allowed and superimposes well.
301         * Caches results.
302         * @param permutation
303         * @return null if invalid, or a rotation if valid
304         */
305        private Rotation isValidPermutation(List<Integer> permutation) {
306                if (permutation.size() == 0) {
307                        return null;
308                }
309
310                // cached value
311                if (evaluatedPermutations.containsKey(permutation)) {
312                        return evaluatedPermutations.get(permutation);
313                }
314
315                // check if permutation is allowed
316                if (! isAllowedPermutation(permutation)) {
317                        evaluatedPermutations.put(permutation, null);
318                        return null;
319                }
320
321                // check if superimposes
322                Rotation rot = superimposePermutation(permutation);
323                return rot;
324        }
325
326        /**
327         * The permutation must map all subunits onto an equivalent subunit
328         * and no subunit onto itself
329         * @param permutation
330         * @return
331         */
332        private boolean isAllowedPermutation(List<Integer> permutation) {
333                List<Integer> seqClusterId = subunits.getClusterIds();
334                int selfaligned = 0;
335                for (int i = 0; i < permutation.size(); i++) {
336                        int j = permutation.get(i);
337                        if ( seqClusterId.get(i) != seqClusterId.get(j)) {
338                                return false;
339                        }
340                        if(i == j ) {
341                                selfaligned++;
342                        }
343                }
344                // either identity (all self aligned) or all unique
345                return selfaligned == 0 || selfaligned == permutation.size();
346        }
347        /**
348         * Adds translational component to rotation matrix
349         * @param rotation
350         * @return
351         */
352        private void combineWithTranslation(Matrix4d rotation) {
353                rotation.setTranslation(centroid);
354                rotation.mul(rotation, centroidInverse);
355        }
356
357        private static Rotation createSymmetryOperation(List<Integer> permutation, Matrix4d transformation, AxisAngle4d axisAngle, int fold, QuatSymmetryScores scores) {
358                Rotation s = new Rotation();
359                s.setPermutation(new ArrayList<Integer>(permutation));
360                s.setTransformation(new Matrix4d(transformation));
361                s.setAxisAngle(new AxisAngle4d(axisAngle));
362                s.setFold(fold);
363                s.setScores(scores);
364
365                return s;
366        }
367
368        private void setupDistanceBox() {
369                distanceThreshold = calcDistanceThreshold();
370                box = new DistanceBox<Integer>(distanceThreshold);
371
372                for (int i = 0; i < originalCoords.length; i++) {
373                        box.addPoint(originalCoords[i], i);
374                }
375        }
376
377        private double calcDistanceThreshold() {
378                double threshold = Double.MAX_VALUE;
379                int n = subunits.getSubunitCount();
380                List<Point3d> centers = subunits.getCenters();
381
382                for (int i = 0; i < n - 1; i++) {
383                        Point3d pi = centers.get(i);
384                        for (int j = i + 1; j < n; j++) {
385                                Point3d pj = centers.get(j);
386                                threshold = Math.min(threshold, pi.distanceSquared(pj));
387                        }
388                }
389                double distanceThreshold = Math.sqrt(threshold);
390
391                distanceThreshold = Math.max(distanceThreshold, parameters.getRmsdThreshold());
392
393                return distanceThreshold;
394        }
395
396        /**
397         * Compare this.transformedCoords with the original coords. For each
398         * subunit, return the transformed subunit with the closest position.
399         * @return A list mapping each subunit to the closest transformed subunit
400         */
401        private List<Integer> getPermutation() {
402                List<Integer> permutation = new ArrayList<Integer>(transformedCoords.length);
403                double sum = 0.0f;
404
405                for (Point3d t: transformedCoords) {
406                        List<Integer> neighbors = box.getNeighborsWithCache(t);
407                        int closest = -1;
408                        double minDist = Double.MAX_VALUE;
409
410                         for (int j : neighbors) {
411                                double dist = t.distanceSquared(originalCoords[j]);
412                                if (dist < minDist) {
413                                        closest = j;
414                                        minDist = dist;
415                                }
416                        }
417
418                        sum += minDist;
419                        if (closest == -1) {
420                                 break;
421                        }
422                        permutation.add(closest);
423                }
424                double rmsd = Math.sqrt(sum / transformedCoords.length);
425
426                if (rmsd > distanceThreshold || permutation.size() != transformedCoords.length) {
427                        permutation.clear();
428                        return permutation;
429                }
430
431                // check uniqueness of indices
432                Set<Integer> set = new HashSet<Integer>(permutation);
433
434                // if size mismatch, clear permutation (its invalid)
435                if (set.size() != originalCoords.length) {
436                        permutation.clear();
437                }
438
439                return permutation;
440        }
441
442        private void initialize() {
443                // translation to centered coordinate system
444                centroid = new Vector3d(subunits.getCentroid());
445                // translation back to original coordinate system
446                Vector3d reverse = new Vector3d(centroid);
447                reverse.negate();
448                centroidInverse.set(reverse);
449                // Make sure matrix element m33 is 1.0. An old version vecmath did not set this element.
450                centroidInverse.setElement(3, 3,  1.0);
451
452                List<Point3d> centers = subunits.getCenters();
453                int n = subunits.getSubunitCount();
454
455                originalCoords = new Point3d[n];
456                transformedCoords = new Point3d[n];
457
458                for (int i = 0; i < n; i++) {
459                        originalCoords[i] = centers.get(i);
460                        transformedCoords[i] = new Point3d();
461                }
462
463                setupDistanceBox();
464        }
465}