diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/PermutationGroup.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/PermutationGroup.java index 675e4facb0..19f10ba72e 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/PermutationGroup.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/PermutationGroup.java @@ -23,6 +23,7 @@ import java.util.ArrayList; import java.util.HashSet; +import java.util.Iterator; import java.util.List; import java.util.Set; @@ -30,7 +31,7 @@ * * @author Peter */ -public class PermutationGroup { +public class PermutationGroup implements Iterable> { List> permutations = new ArrayList>(); public void addPermutation(List permutation) { @@ -139,4 +140,9 @@ public String getGroupTable() { public int hashCode() { return getGroupTable().hashCode(); } + + @Override + public Iterator> iterator() { + return permutations.iterator(); + } } diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationGroup.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationGroup.java index ead65b5b20..c605aa29e7 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationGroup.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationGroup.java @@ -27,13 +27,14 @@ import java.util.ArrayList; import java.util.Collections; import java.util.Comparator; +import java.util.Iterator; import java.util.List; /** * * @author Peter */ -public class RotationGroup { +public class RotationGroup implements Iterable { private List rotations = new ArrayList(); private int principalAxisIndex = 0; private int higherOrderRotationAxis = 0; @@ -369,4 +370,9 @@ public int compare(Rotation o1, Rotation o2) { } }); } + + @Override + public Iterator iterator() { + return rotations.iterator(); + } } diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationSolver.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationSolver.java index 8802ef3bbc..bbc0d4cf9e 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationSolver.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationSolver.java @@ -21,19 +21,22 @@ package org.biojava.nbio.structure.symmetry.core; -import org.biojava.nbio.structure.symmetry.geometry.DistanceBox; -import org.biojava.nbio.structure.symmetry.geometry.MomentsOfInertia; -import org.biojava.nbio.structure.symmetry.geometry.SphereSampler; -import org.biojava.nbio.structure.symmetry.geometry.SuperPosition; +import java.util.ArrayList; +import java.util.HashMap; +import java.util.HashSet; +import java.util.List; +import java.util.Map; +import java.util.Set; import javax.vecmath.AxisAngle4d; import javax.vecmath.Matrix4d; import javax.vecmath.Point3d; import javax.vecmath.Vector3d; -import java.util.ArrayList; -import java.util.HashSet; -import java.util.List; -import java.util.Set; + +import org.biojava.nbio.structure.symmetry.geometry.DistanceBox; +import org.biojava.nbio.structure.symmetry.geometry.MomentsOfInertia; +import org.biojava.nbio.structure.symmetry.geometry.SphereSampler; +import org.biojava.nbio.structure.symmetry.geometry.SuperPosition; /** @@ -50,7 +53,8 @@ public class RotationSolver implements QuatSymmetrySolver { private Matrix4d centroidInverse = new Matrix4d(); private Point3d[] originalCoords = null; private Point3d[] transformedCoords = null; - private Set> hashCodes = new HashSet>(); + // Cache whether a permutation is invalid (null) vs has been added to rotations + private Map,Rotation> evaluatedPermutations = new HashMap<>(); private RotationGroup rotations = new RotationGroup(); @@ -95,9 +99,14 @@ private void solve() { List angles = getAngles(); - for (int i = 0; i < sphereCount; i++) { + for (int i = 0; i < sphereCount; i++) { + // Sampled orientation + //TODO The SphereSampler samples 4D orientation space. We really + // only need to sample 3D unit vectors, since we use limited + // angles. -SB SphereSampler.getAxisAngle(i, sphereAngle); + // Each valid rotation angle for (double angle : angles) { // apply rotation sphereAngle.angle = angle; @@ -113,14 +122,14 @@ private void solve() { List permutation = getPermutation(); // System.out.println("Rotation Solver: permutation: " + i + ": " + permutation); - boolean isValidPermuation = isValidPermutation(permutation); - if (! isValidPermuation) { - continue; + // check if novel + if ( evaluatedPermutations.containsKey(permutation)) { + continue; //either invalid or already added } - boolean newPermutation = evaluatePermutation(permutation); - if (newPermutation) { - completeRotationGroup(); + Rotation newPermutation = isValidPermutation(permutation); + if (newPermutation != null) { + completeRotationGroup(newPermutation); } // check if all symmetry operations have been found. @@ -131,33 +140,84 @@ private void solve() { } } - private void completeRotationGroup() { + /** + * Combine current rotations to make all possible permutations. + * If these are all valid, add them to the rotations + * @param additionalRots Additional rotations we are considering adding to this.rotations + * @return whether the rotations were valid and added + */ + private boolean completeRotationGroup(Rotation... additionalRots) { PermutationGroup g = new PermutationGroup(); - for (int i = 0; i < rotations.getOrder(); i++) { - Rotation s = rotations.getRotation(i); + for (Rotation s : rotations) { g.addPermutation(s.getPermutation()); } + for( Rotation s : additionalRots) { + g.addPermutation(s.getPermutation()); + // inputs should not have been added already + assert evaluatedPermutations.get(s.getPermutation()) == null; + } g.completeGroup(); // the group is complete, nothing to do - if (g.getOrder() == rotations.getOrder()) { - return; + if (g.getOrder() == rotations.getOrder()+additionalRots.length) { + for( Rotation s : additionalRots) { + addRotation(s); + } + return true; } // try to complete the group - for (int i = 0; i < g.getOrder(); i++) { - List permutation = g.getPermutation(i); - - boolean isValidPermutation = isValidPermutation(permutation); - if (isValidPermutation) { + List newRots = new ArrayList<>(g.getOrder()); + // First, quick check for whether they're allowed + for (List permutation : g) { + if( evaluatedPermutations.containsKey(permutation)) { + Rotation rot = evaluatedPermutations.get(permutation); + if( rot == null ) { + return false; + } + } else { + if( ! isAllowedPermutation(permutation)) { + return false; + } + } + } + // Slower check including the superpositions + for (List permutation : g) { + Rotation rot; + if( evaluatedPermutations.containsKey(permutation)) { + rot = evaluatedPermutations.get(permutation); + } else { + rot = isValidPermutation(permutation); + } - // perform permutation of subunits - evaluatePermutation(permutation); + if( rot == null ) { + // if any induced rotation is invalid, abort + return false; + } + if(!evaluatedPermutations.containsKey(permutation)){ //novel + newRots.add(rot); } } + // Add rotations + for( Rotation rot : newRots) { + addRotation(rot); + } + return true; } - private boolean evaluatePermutation(List permutation) { + private void addRotation(Rotation rot) { + evaluatedPermutations.put(rot.getPermutation(),rot); + rotations.addRotation(rot); + } + + /** + * Superimpose subunits based on the given permutation. Then check whether + * the superposition passes RMSD thresholds and create a Rotation to + * represent it if so. + * @param permutation A list specifying which subunits should be aligned by the current transformation + * @return A Rotation representing the permutation, or null if the superposition did not meet thresholds. + */ + private Rotation superimposePermutation(List permutation) { // permutate subunits for (int j = 0, n = subunits.getSubunitCount(); j < n; j++) { transformedCoords[j].set(originalCoords[permutation.get(j)]); @@ -176,17 +236,20 @@ private boolean evaluatePermutation(List permutation) { // evaluate superposition of CA traces QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation); if (scores.getRmsd() < 0.0 || scores.getRmsd() > parameters.getRmsdThreshold()) { - return false; + return null; } scores.setRmsdCenters(subunitRmsd); Rotation symmetryOperation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores); - rotations.addRotation(symmetryOperation); - return true; + return symmetryOperation; } - return false; + return null; } + /** + * Get valid rotation angles given the number of subunits + * @return The rotation angle corresponding to each fold of {@link Subunits#getFolds()} + */ private List getAngles() { int n = subunits.getSubunitCount(); // for spherical symmetric cases, n cannot be higher than 60 @@ -210,44 +273,53 @@ private boolean isSpherical() { return m.getSymmetryClass(0.05) == MomentsOfInertia.SymmetryClass.SYMMETRIC; } - private boolean isValidPermutation(List permutation) { - if (permutation.size() == 0) { - return false; - } + /** + * Checks if a particular permutation is allowed and superimposes well. + * Caches results. + * @param permutation + * @return null if invalid, or a rotation if valid + */ + private Rotation isValidPermutation(List permutation) { + if (permutation.size() == 0) { + return null; + } - // if this permutation is a duplicate, return false - if (hashCodes.contains(permutation)) { - return false; + // cached value + if (evaluatedPermutations.containsKey(permutation)) { + return evaluatedPermutations.get(permutation); } // check if permutation is allowed if (! isAllowedPermutation(permutation)) { - return false; - } - - // get fold and make sure there is only one E (fold=1) permutation - int fold = PermutationGroup.getOrder(permutation); - if (rotations.getOrder() > 1 && fold == 1) { - return false; - } - - if (fold == 0 || subunits.getSubunitCount() % fold != 0) { - return false; + evaluatedPermutations.put(permutation, null); + return null; } - // if this permutation is a duplicate, returns false - return hashCodes.add(permutation); + // check if superimposes + Rotation rot = superimposePermutation(permutation); + return rot; } + /** + * The permutation must map all subunits onto an equivalent subunit + * and no subunit onto itself + * @param permutation + * @return + */ private boolean isAllowedPermutation(List permutation) { List seqClusterId = subunits.getSequenceClusterIds(); + int selfaligned = 0; for (int i = 0; i < permutation.size(); i++) { int j = permutation.get(i); - if (seqClusterId.get(i) != seqClusterId.get(j)) { + if ( seqClusterId.get(i) != seqClusterId.get(j)) { return false; } + if(i == j ) { + selfaligned++; + } } - return true; + // either identity (all self aligned) or all unique + return selfaligned == 0 || selfaligned == permutation.size(); } /** * Adds translational component to rotation matrix @@ -260,7 +332,7 @@ private void combineWithTranslation(Matrix4d rotation) { rotation.mul(rotation, centroidInverse); } - private Rotation createSymmetryOperation(List permutation, Matrix4d transformation, AxisAngle4d axisAngle, int fold, QuatSymmetryScores scores) { + private static Rotation createSymmetryOperation(List permutation, Matrix4d transformation, AxisAngle4d axisAngle, int fold, QuatSymmetryScores scores) { Rotation s = new Rotation(); s.setPermutation(new ArrayList(permutation)); s.setTransformation(new Matrix4d(transformation)); @@ -299,6 +371,11 @@ private double calcDistanceThreshold() { return distanceThreshold; } + /** + * Compare this.transformedCoords with the original coords. For each + * subunit, return the transformed subunit with the closest position. + * @return A list mapping each subunit to the closest transformed subunit + */ private List getPermutation() { List permutation = new ArrayList(transformedCoords.length); double sum = 0.0f; diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/geometry/IcosahedralSampler.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/geometry/IcosahedralSampler.java index 08803c50a3..650f2d6c65 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/geometry/IcosahedralSampler.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/geometry/IcosahedralSampler.java @@ -24,9 +24,12 @@ import javax.vecmath.Quat4d; /** - * + * Represents an even coverage of quaternion space by 60 points. This grid is + * defined by the vertices of one half of a hexacosichoron (600-cell). It + * approximates all possible orientations to within 44.48 degrees. * @author Peter */ +// This would be better named HexacosichoronSampler, since it's sampling 4D space. -SB public final class IcosahedralSampler { private static Quat4d quat = new Quat4d(); @@ -39,7 +42,7 @@ public static int getSphereCount() { } public static Quat4d getQuat4d(int index) { - Quat4d q = new Quat4d(orientations[index]); + Quat4d q = new Quat4d(orientations[index]); //ignores 5th element return q; } @@ -48,10 +51,12 @@ public static void getAxisAngle(int index, AxisAngle4d axisAngle) { axisAngle.set(quat); } -// # Orientation set c600v, number = 60, radius = 44.48 degrees -// # $Id: c600v.quat 6102 2006-02-21 19:45:40Z ckarney $ -// # For more information, eee http://charles.karney.info/orientation/ -// format quaternion + // # Orientation set c600v, number = 60, radius = 44.48 degrees + // # $Id: c600v.quat 6102 2006-02-21 19:45:40Z ckarney $ + // # For more information, eee http://charles.karney.info/orientation/ + // format quaternion + // The fifth column gives a weighting factor. Since the 600-cell is regular, all + // orientations cover an equal fraction of orientation space and have equal weight. private static double[][] orientations = { {1.000000000f, 0.000000000f, 0.000000000f, 0.000000000f, 1.000000f}, {0.000000000f, 1.000000000f, 0.000000000f, 0.000000000f, 1.000000f}, diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/geometry/SphereSampler.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/geometry/SphereSampler.java index d6253eb4de..032ca4b957 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/geometry/SphereSampler.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/geometry/SphereSampler.java @@ -25,7 +25,11 @@ import java.util.List; -// Generate the permutations and sign changes for a Triple. +/** + * Generate the permutations and sign changes for a Triple. + * For instance, (a,0,0) becomes (±a,0,0),(0,±a,0),(0,0,±a). In the general + * case 48 tuples are returned. + */ class Permute { private List triples = new ArrayList(); @@ -38,6 +42,7 @@ class Permute { triples.add(tmp); int n = 1; + // Unpermuted ±x,±y,±z if (t.x != 0) { for (int i = 0; i < n; ++i) { Tuple3i m = triples.get(i); @@ -66,6 +71,7 @@ class Permute { return; } + // At least one differing value for (int i = 0; i < n; ++i) { Point3i m = triples.get(i); triples.add(new Point3i(m.y, m.z, m.x)); @@ -77,6 +83,7 @@ class Permute { return; } + // At least two differing values for (int i = 0; i < n; ++i) { Point3i m = triples.get(i); triples.add(new Point3i(m.y, m.x, m.z)); @@ -94,7 +101,10 @@ public Point3i get(int i) { }; /** - * + * Sample possible orientations. An orientation can be represented as a + * quaternion or as a rotation axis and angle. The orientations are somewhat + * evenly distributed over orientation space, but the current implementation + * is not suited for statistically uniform sampling. * @author Peter */ public final class SphereSampler { @@ -125,16 +135,20 @@ public final class SphereSampler { + // Approximate spacing between grid points private static final double delta = 0.15846; + // Amount of distortion towards grid corners to account for the projection back to a 4-sphere private static final double sigma = 0.00; - private static final int ntot = 7416; - private static final int ncell = 309; - private static final int nent = 18; + private static final int ntot = 7416; // total number of grid points over the 24 cells + private static final int ncell = 309; // number of grid points per cell + private static final int nent = 18; // number of distinct grid points (with k>=l>=m) + // Why not (5,5,3) and (5,5,5)? Would they overlap with other cells? -SB private static final int[] k = { 0, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5}; private static final int[] l = { 0, 1, 0, 2, 2, 1, 3, 3, 0, 2, 2, 4, 4, 4, 1, 3, 3, 5}; private static final int[] m = { 0, 1, 0, 0, 2, 1, 1, 3, 0, 0, 2, 0, 2, 4, 1, 1, 3, 1}; + // number of permutations to expect for each (k,l,m) tuple private static final int[] mult = { 1, 8, 6, 12, 8, 24, 24, 8, 6, 24, 24, 12, 24, 8, 24, 48, 24, 24}; static @@ -142,15 +156,29 @@ public final class SphereSampler { List myorientations = new ArrayList(); + // 60 equally spaced grid points covering all quaternions + // This could be removed, since the 48-cell bcc lattice has lower angle error -SB for (int i = 0; i < IcosahedralSampler.getSphereCount(); i++) { myorientations.add(IcosahedralSampler.getQuat4d(i)); } + + // SB's interpretation of this code: + // Directly form a 4D grid over the 4-sphere of orientations. Start by + // making a 3D body-centered lattice with w=1. Then use the hypercube symmetries + // to extend this grid over the whole surface. + // The permuted (k,l,m) values together make a diagonal grid out to ±(5,5,5). + // The point spacing is distorted by the pind() function so that the + // projection of the points back to the 4-sphere will be more even. + + // This is the c48u309 lattice from Karney 2006, with a max error of 10.07 + // degrees. + List grid = new ArrayList(); int ncell1 = 0; - for (int n = 0; n < nent; ++n) { + for (int n = 0; n < nent; ++n) { // for each tuple (k,l,m) above Permute p = new Permute(new Point3i(k[n], l[n], m[n])); assert (mult[n] == p.size()); - for (int i = 0; i < mult[n]; ++i) { + for (int i = 0; i < mult[n]; ++i) { // for each permutation Point3i t = p.get(i); grid.add(new Quat4d(1.0, pind(0.5 * t.x, delta, sigma), pind( 0.5 * t.y, delta, sigma), pind(0.5 * t.z, delta, sigma))); @@ -158,6 +186,7 @@ public final class SphereSampler { ncell1 += mult[n]; } assert (ncell1 == ncell); + // Apply symmetry operations to distribute grid over all values of w int nc = grid.size(); assert (nc == ncell); for (int n = 1; n < 24; ++n) { @@ -168,6 +197,7 @@ public final class SphereSampler { qs.mul(q, grid.get(i)); grid.add(qs); // s.add(times(q, s.getOrientation(i)), s.getWeight(i)); // this data set has no weights + // Actually, it does have weights but we don't care since we don't need uniform sampling. -SB } } assert (grid.size() == ntot); @@ -196,6 +226,10 @@ public static void getAxisAngle(int index, AxisAngle4f axisAngle) { axisAngle.set(orientations.get(index)); } + /** + * @param index Index between 0 and {@link #getSphereCount()}-1 + * @param axisAngle (Output) modified to contain the index-th sampled orientation + */ public static void getAxisAngle(int index, AxisAngle4d axisAngle) { axisAngle.set(orientations.get(index));