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 */ 021package org.biojava.nbio.structure.symmetry.internal; 022 023import java.util.ArrayList; 024import java.util.Arrays; 025import java.util.Iterator; 026import java.util.List; 027 028import javax.vecmath.Matrix4d; 029 030import org.biojava.nbio.structure.align.multiple.MultipleAlignment; 031import org.biojava.nbio.structure.align.util.RotationAxis; 032import org.biojava.nbio.structure.symmetry.internal.CESymmParameters.SymmetryType; 033 034/** 035 * Data Structure that stores all the symmetry axis that describe 036 * the symmetry of a structure. Generalizes to all types of symmetry, 037 * the classic ones (Cn, Dn) and any hierarchical or local symmetries. 038 * <p> 039 * Hierarchical symmetry can be visualized as a tree, where each level 040 * has a fixed branching factor. Each level of the tree is associated 041 * with a transformation operator, whose order determines the degree of 042 * nodes at that level of the tree. Leaves of the tree implicitly 043 * represent aligned repeats (indexed 0 to n-1), so care must be taken to 044 * keep external references to the repeats (e.g. rows of a 045 * {@link MultipleAlignment} in the same order implied by the tree. 046 * <p> 047 * Each node of the tree specifies an alignment between those repeats 048 * below each of its children. It is also associated with a symmetry axis, 049 * which is calculated based on the associated operator as well as any parent 050 * operators. 051 * It also stores the parts of the structure (symmetric units) involved 052 * in each axis, in addition to the way to calculate them. 053 * <p> 054 * This is intended to provide a general axis support for the multiple 055 * repeat alignment optimization and the axis display in Jmol. This 056 * object is related to a MultipleAlignment object that defines the 057 * symmetric units. 058 * 059 * @author Aleix Lafita 060 * @since 4.2.0 061 * 062 */ 063public class SymmetryAxes { 064 /* 065 * Implementation note: The tree is a nice explanation and a good image 066 * for developing algorithms, but it is not constructed explicitly. 067 * Instead, we just store one elementary axis for each level and reconstruct 068 * which operators apply to a particular leaf based on that leaf's index. 069 */ 070 071 /** 072 * Represents an axis of symmetry 073 * @author Spencer Bliven 074 * 075 */ 076 public static class Axis { 077 private Matrix4d operator; 078 private int order; 079 private SymmetryType symmType; 080 private int level; 081 //private int indexInLevel; 082 private int firstRepeat; 083 private RotationAxis rotAxis; 084 085 public Axis(Matrix4d operator, int order, SymmetryType type, int level, int firstRepeat) { 086 if (order < 2) { 087 throw new IllegalArgumentException("A symmetry axis should divide a structure in > 2 parts"); 088 } 089 if(type != SymmetryType.OPEN && type != SymmetryType.CLOSED) { 090 throw new IllegalArgumentException("Invalid symmetry type. Only OPEN and CLOSED are allowed"); 091 } 092 093 this.operator = operator; 094 this.order = order; 095 this.symmType = type; 096 setLevel(level); 097 setFirstRepeat(firstRepeat); 098 rotAxis = null; 099 } 100 /** 101 * Get the transformation operator for this axis as an homogeneous matrix 102 * @return the transformation operator 103 */ 104 public Matrix4d getOperator() { 105 return operator; 106 } 107 public void setOperator(Matrix4d op) { 108 this.operator = op; 109 } 110 /** 111 * Get the order of this axis (closed symm) or the number of repeats 112 * (open symm) 113 * @return the order 114 */ 115 public int getOrder() { 116 return order; 117 } 118 /** 119 * @return the symmType (OPEN or CLOSED only) 120 */ 121 public SymmetryType getSymmType() { 122 return symmType; 123 } 124 125 /** 126 * Get the transformation operator as a rotation axis. For open 127 * symmetry this will have a non-zero screw component. 128 * @return a RotationAxis for this Axis 129 */ 130 public RotationAxis getRotationAxis() { 131 if( rotAxis == null) { 132 rotAxis = new RotationAxis(operator); 133 } 134 return rotAxis; 135 } 136 /** 137 * @return The level of this axis within it's parent hierarchy, or -1 if unset 138 */ 139 public int getLevel() { 140 return level; 141 } 142 /** 143 * 144 * @param level The level of this axis within it's parent hierarchy. Must be positive 145 */ 146 public void setLevel(int level) { 147 if(level < 0) throw new IndexOutOfBoundsException("Level must be positive"); 148 this.level = level; 149 } 150// /** 151// * Each level can contain multiple equivalent axes. This index is 152// * used to distinguish them. 153// * @return the index of this axis relative to others at the same level 154// */ 155// public int getIndexInLevel() { 156// return indexInLevel; 157// } 158// /** 159// * 160// * @param indexInLevel the index of this axis relative to others at the same level 161// */ 162// public void setIndexInLevel(int indexInLevel) { 163// if( indexInLevel < 0 || getOrder() <= indexInLevel ) 164// throw new IndexOutOfBoundsException("Invalid index for order "+getOrder()); 165// this.indexInLevel = indexInLevel; 166// } 167 /** 168 * Get the index of the first repeat used by this axis 169 * @return the firstRepeat 170 */ 171 public int getFirstRepeat() { 172 return firstRepeat; 173 } 174 /** 175 * @param firstRepeat the index of the first repeat used by this axis 176 */ 177 public void setFirstRepeat(int firstRepeat) { 178 this.firstRepeat = firstRepeat; 179 } 180 } 181 182 /** 183 * List of all symmetry axis. They are sorted from higher to lower 184 * in the symmetry hierarchy, where higher means that they apply 185 * more globally and lower means that they apply to a local region 186 * of the higher axis division. 187 */ 188 private final List<Axis> axes; 189 190 /** 191 * Constructor. 192 * Initializes variables only. 193 */ 194 public SymmetryAxes(){ 195 axes = new ArrayList<>(); 196 } 197 198 /** 199 * Adds a new axis of symmetry to the bottom level of the tree 200 * 201 * @param axis the new axis of symmetry found 202 * @param order number of parts that this axis divides the structure in 203 * @param type indicates whether the axis has OPEN or CLOSED symmetry 204 */ 205 public void addAxis(Matrix4d axis, int order, SymmetryType type) { 206 axes.add(new Axis(axis,order,type,axes.size(),0)); 207 } 208 209 /** 210 * Return a list giving the number of times each axis must be applied 211 * to generate the given repeat. 212 * <P> 213 * For instance, for a D3 case <tt>getAxisCounts(4)</tt> would return [2,0], 214 * indicating that repeat 4 is generated by two applications of the 3-fold 215 * axis followed by 0 applications of the two-fold axis. 216 * 217 * @param repeat Index of the desired repeat 218 * @return array of the same length as axes giving the number of times 219 * to apply each axis. 220 */ 221 private int[] getAxisCounts(int repeat) { 222 int[] counts = new int[getNumLevels()]; 223 224 for(int i = counts.length-1; i >= 0; i--) { 225 int d = axes.get(i).getOrder(); 226 counts[i] = repeat % d; 227 repeat /= d; 228 } 229 assert repeat == 0 : "Invalid repeat index"; 230 return counts; 231 } 232 233// /** 234// * Inverse of {@link #getAxisCounts(int)}; Calculates the repeat for a 235// * particular number of applications of each axis 236// * @param counts Number of times to apply each axis 237// * @return Repeat index 238// */ 239// private int getRepeatIndex(int[] counts) { 240// int repeat = 0; 241// for(int i = 0; i< counts.length; i++) { 242// repeat += counts[i]*axes.get(i).getOrder(); 243// } 244// return repeat; 245// } 246 /** 247 * Updates an axis of symmetry, after the superposition changed. 248 * 249 * @param index old axis index 250 * @param newAxis 251 */ 252 public void updateAxis(Integer index, Matrix4d newAxis){ 253 axes.get(index).setOperator(newAxis); 254 } 255 256 /** 257 * Return the operator for all elementary axes of symmetry of the structure, that is, 258 * the axes stored in the List as unique and from which all the symmetry 259 * axes are constructed. 260 * 261 * @return axes elementary axes of symmetry. 262 */ 263 public List<Matrix4d> getElementaryAxes(){ 264 List<Matrix4d> ops = new ArrayList<Matrix4d>(getNumLevels()); 265 for(Axis axis : axes) { 266 ops.add(axis.getOperator()); 267 } 268 return ops; 269 } 270 271 /** 272 * Return all elementary axes of symmetry of the structure, that is, 273 * the axes stored in the List as unique and from which all the symmetry 274 * axes are constructed. 275 * 276 * @return axes elementary axes of symmetry. 277 */ 278 public List<Axis> getElementaryAxesObjects() { 279 return axes; 280 } 281 282 /** 283 * Get the indices of participating repeats in Cauchy two-line form. 284 * <p> 285 * Returns two lists of the same length. 286 * The first gives a list of all repeat indices which are aligned 287 * at the specified level of symmetry (e.g. 0 through the degree of this level). 288 * The second list gives the corresponding repeats after applying the 289 * operator once. 290 * 291 * @param level the axis index 292 * @return the double List of repeat relations, or null if the 293 * level is invalid 294 * @see #getRepeatsCyclicForm(int, int) for an equivalent specification with half the memory 295 */ 296 public List<List<Integer>> getRepeatRelation(int level){ 297 return getRepeatRelation(level,0); 298 } 299 300 public List<List<Integer>> getRepeatRelation(Axis axis){ 301 return getRepeatRelation(axis.getLevel(),axis.getFirstRepeat()); 302 } 303 304 public List<List<Integer>> getRepeatRelation(int level, int firstRepeat) { 305 Axis axis = axes.get(level); 306 int m = getNumRepeats(level+1);//size of the children 307 int d = axis.getOrder(); // degree of this node 308 int n = m*d; // number of repeats included 309 if(firstRepeat % n != 0) 310 throw new IllegalArgumentException(String.format("Repeat %d cannot start a block at level %s of this tree",firstRepeat,level)); 311 if(axis.getSymmType() == SymmetryType.OPEN) { 312 n -= m; // leave off last child for open symm 313 } 314 List<Integer> repeats = new ArrayList<>(n); 315 List<Integer> equiv = new ArrayList<>(n); 316 for(int i=0;i<n;i++) { 317 repeats.add(i+firstRepeat); 318 equiv.add( (i+m)%(m*d)+firstRepeat ); 319 } 320 return Arrays.asList(repeats,equiv); 321 322 } 323 324 /** 325 * Get the indices of participating repeats in cyclic form. 326 * <p> 327 * Each inner list gives a set of equivalent repeats and should have length 328 * equal to the order of the axis' operator. 329 * @param level 330 * @param firstRepeat 331 * @return 332 */ 333 public List<List<Integer>> getRepeatsCyclicForm(int level, int firstRepeat) { 334 Axis axis = axes.get(level); 335 int m = getNumRepeats(level+1);//size of the children 336 int d = axis.getOrder(); // degree of this node 337 int n = m*d; // number of repeats included 338 if(firstRepeat % n != 0) { 339 throw new IllegalArgumentException(String.format("Repeat %d cannot start a block at level %s of this tree",firstRepeat,level)); 340 } 341 if(axis.getSymmType() == SymmetryType.OPEN) { 342 n -= m; // leave off last child for open symm 343 } 344 345 List<List<Integer>> repeats = new ArrayList<>(m); 346 for(int i=0;i<m;i++) { 347 List<Integer> cycle = new ArrayList<>(d); 348 for(int j=0;j<d;j++) { 349 cycle.add(firstRepeat+i+j*m); 350 } 351 repeats.add(cycle); 352 } 353 return repeats; 354 } 355 public List<List<Integer>> getRepeatsCyclicForm(Axis axis) { 356 return getRepeatsCyclicForm(axis.getLevel(),axis.getFirstRepeat()); 357 } 358 public List<List<Integer>> getRepeatsCyclicForm(int level) { 359 return getRepeatsCyclicForm(level,0); 360 } 361 public String getRepeatsCyclicForm(Axis axis, List<?> repeats) { 362 if(repeats.size() != getNumRepeats()) { 363 throw new IllegalArgumentException("Mismatch in the number of repeats"); 364 } 365 return getRepeatsCyclicForm(getRepeatsCyclicForm(axis), repeats); 366 } 367 public static String getRepeatsCyclicForm(List<List<Integer>> cycleForm, List<?> repeats) { 368 StringBuilder str = new StringBuilder(); 369 for(List<Integer> cycle : cycleForm) { 370 str.append("("); 371 Iterator<Integer> cycleIt = cycle.iterator(); 372 str.append(repeats.get(cycleIt.next())); //should be at least one 373 while(cycleIt.hasNext()) { 374 str.append(";") 375 .append(repeats.get( cycleIt.next() )); 376 } 377 str.append(")"); 378 } 379 return str.toString(); 380 } 381 382 /** 383 * Return the transformation that needs to be applied to a 384 * repeat in order to superimpose onto repeat 0. 385 * 386 * @param repeat the repeat index 387 * @return transformation matrix for the repeat 388 */ 389 public Matrix4d getRepeatTransform(int repeat){ 390 391 Matrix4d transform = new Matrix4d(); 392 transform.setIdentity(); 393 394 int[] counts = getAxisCounts(repeat); 395 396 for(int t = counts.length-1; t>=0; t--) { 397 if( counts[t] == 0 ) 398 continue; 399 Matrix4d axis = new Matrix4d(axes.get(t).getOperator()); 400 for(int i=0;i<counts[t];i++) { 401 transform.mul(axis); 402 } 403 } 404 return transform; 405 } 406 407 /** 408 * Return the transformation that needs to be applied to 409 * repeat x in order to superimpose onto repeat y. 410 * 411 * @param x the first repeat index (transformed) 412 * @param y the second repeat index (fixed) 413 * @return transformation matrix for the repeat x 414 */ 415 public Matrix4d getRepeatTransform(int x, int y){ 416 417 Matrix4d transform = new Matrix4d(); 418 transform.setIdentity(); 419 420 int[] iCounts = getAxisCounts(x); 421 int[] jCounts = getAxisCounts(y); 422 423 int[] counts = new int[iCounts.length]; 424 for (int k = 0; k < iCounts.length; k++) 425 counts[k] = iCounts[k] - jCounts[k]; 426 427 for(int t = counts.length-1; t>=0; t--) { 428 if(counts[t] == 0) 429 continue; 430 if (counts[t] > 0) { 431 Matrix4d axis = new Matrix4d(axes.get(t).getOperator()); 432 for(int i=0;i<counts[t];i++) 433 transform.mul(axis); 434 } else if (counts[t] < 0) { 435 Matrix4d axis = new Matrix4d(axes.get(t).getOperator()); 436 axis.invert(); 437 for(int i=0;i<counts[t];i++) 438 transform.mul(axis); 439 } 440 } 441 return transform; 442 } 443 444 /** 445 * Return all symmetry axes of of the structure: the set of axes that 446 * describe all parts of the structure. This combines the elementary 447 * axes to generate all possible axes. The axes are returned in the repeat 448 * degrees. 449 * @return axes all symmetry axes of the structure. 450 */ 451 public List<Axis> getSymmetryAxes(){ 452 453 List<Axis> symmAxes = new ArrayList<>(); 454 455 Matrix4d prior = new Matrix4d(); 456 prior.setIdentity(); 457 458 getSymmetryAxes(symmAxes,prior,0,0); 459 460 461 return symmAxes; 462 } 463 /** 464 * Recursive helper 465 * @param symmAxes output list 466 * @param prior transformation aligning the first repeat of this axis with the first overall 467 * @param level current level 468 */ 469 private void getSymmetryAxes(List<Axis> symmAxes, Matrix4d prior, int level, int firstRepeat) { 470 if(level >= getNumLevels() ) { 471 return; 472 } 473 474 Axis elem = axes.get(level); 475 Matrix4d elemOp = elem.getOperator(); 476 477 // Current axis: 478 // elementary maps B -> A 479 // prior maps I -> A and J -> B 480 // want J -> I = J -> B -> A <- I= inv(prior) * elementary * prior 481 Matrix4d currAxisOp = new Matrix4d(prior); 482 currAxisOp.invert(); 483 currAxisOp.mul(elemOp); 484 currAxisOp.mul(prior); 485 Axis currAxis = new Axis(currAxisOp,elem.getOrder(),elem.getSymmType(),level,firstRepeat); 486 symmAxes.add(currAxis); 487 488 //Remember that all degrees are at least 2 489 getSymmetryAxes(symmAxes,prior,level+1,firstRepeat); 490 //New prior is elementary^d*prior 491 Matrix4d newPrior = new Matrix4d(elemOp); 492 newPrior.mul(prior); 493 int childSize = getNumRepeats(level+1); 494 getSymmetryAxes(symmAxes,newPrior,level+1,firstRepeat+childSize); 495 for(int d=2;d<elem.getOrder();d++) { 496 newPrior.mul(elemOp,newPrior); 497 getSymmetryAxes(symmAxes,newPrior,level+1,firstRepeat+childSize*d); 498 } 499 } 500 501 502// public Matrix4d getSymmetryAxis(int level, int axisNum) { 503// if(level == 0) { 504// if( axisNum != 0 ) 505// throw new IndexOutOfBoundsException("Axis number out of bounds"); 506// return axes.get(0); 507// } else { 508// if( axisNum >= degrees.get(level-1) ) 509// throw new IndexOutOfBoundsException("Axis number out of bounds"); 510// // Convert axisNum into a count of 511// 512// } 513 /** 514 * Get the number of repeats. This is equal to the product of all degrees. 515 * @return Number of repeats (leaves of the tree). 516 */ 517 public int getNumRepeats() { 518 return getNumRepeats(0); 519 } 520 521 /** 522 * Get the number of leaves from a node at the specified level. This is 523 * equal to the product of all degrees at or below the level. 524 * @param level level of the tree to cut at 525 * @return Number of repeats (leaves of the tree). 526 */ 527 private int getNumRepeats(int level) { 528 int size = 1; 529 // Return 1 for illegally high level 530 if(level < getNumLevels()) { 531 for(Axis axis : axes.subList(level, getNumLevels())) { 532 size *= axis.getOrder(); 533 } 534 } 535 return size; 536 } 537 538 /** 539 * Get the first repeat index of each axis of a specified level. 540 * @param level level of the tree to cut at 541 * @return List of first Repeats of each index, sorted in ascending order 542 */ 543 public List<Integer> getFirstRepeats(int level) { 544 List<Integer> firstRepeats = new ArrayList<Integer>(); 545 int m = getNumRepeats(level+1); //size of the level 546 int d = axes.get(level).getOrder(); //degree of this level 547 int n = m*d; // number of repeats included in each axis 548 for (int firstRepeat = 0; firstRepeat < getNumRepeats(); firstRepeat+=n) 549 firstRepeats.add(firstRepeat); 550 return firstRepeats; 551 } 552 553 public Axis getElementaryAxis(int level) { 554 return axes.get(level); 555 } 556 557 public int getNumLevels() { 558 return axes.size(); 559 } 560 561}