From 89ac11adfe3850865be26ba951231abe83ab80ef Mon Sep 17 00:00:00 2001 From: Spencer Bliven Date: Fri, 10 Feb 2017 11:29:19 +0100 Subject: [PATCH 1/2] Better rules for Quaternary symm self-alignments Don't allow subunits to be aligned to themselves (except for the identity alignment). --- .../symmetry/core/RotationSolver.java | 20 ++++++++----------- 1 file changed, 8 insertions(+), 12 deletions(-) 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..ccf1cf1e79 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 @@ -95,7 +95,7 @@ private void solve() { List angles = getAngles(); - for (int i = 0; i < sphereCount; i++) { + for (int i = 0; i < sphereCount; i++) { SphereSampler.getAxisAngle(i, sphereAngle); for (double angle : angles) { @@ -225,25 +225,21 @@ private boolean isValidPermutation(List 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; - } - // if this permutation is a duplicate, returns false return hashCodes.add(permutation); } + /** + * 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(); for (int i = 0; i < permutation.size(); i++) { int j = permutation.get(i); - if (seqClusterId.get(i) != seqClusterId.get(j)) { + if (i == j || seqClusterId.get(i) != seqClusterId.get(j)) { return false; } } From ac195327df3be61ea4d6567e3bd125a02a88c794 Mon Sep 17 00:00:00 2001 From: Spencer Bliven Date: Wed, 15 Feb 2017 16:23:10 +0100 Subject: [PATCH 2/2] Fix Quaternary symmetry symmetry detection bug In special cases it was possible to have pairs of axes that pass our quality thresholds but form invalid symmetry groups. Now alignments that superimpose a subunit onto itself are forbidden, which should restrict it to the chiral groups. Also add numerous notes and documentation. --- .../symmetry/core/PermutationGroup.java | 8 +- .../symmetry/core/RotationGroup.java | 8 +- .../symmetry/core/RotationSolver.java | 171 +++++++++++++----- .../symmetry/geometry/IcosahedralSampler.java | 17 +- .../symmetry/geometry/SphereSampler.java | 48 ++++- 5 files changed, 192 insertions(+), 60 deletions(-) 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 ccf1cf1e79..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(); @@ -96,8 +100,13 @@ private void solve() { List angles = getAngles(); 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 void addRotation(Rotation rot) { + evaluatedPermutations.put(rot.getPermutation(),rot); + rotations.addRotation(rot); } - private boolean evaluatePermutation(List permutation) { + /** + * 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,23 +273,31 @@ 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; + 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; } /** @@ -237,13 +308,18 @@ private boolean isValidPermutation(List permutation) { */ 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 (i == j || 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 @@ -256,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)); @@ -295,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));