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