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}