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. 200 * The repeats that participate in this axis and their superposition 201 * relation should also be indicated. 202 * 203 * @param axis the new axis of symmetry found 204 * @param superposition repeats participating and superposition relation 205 * @param repeats number of times the transformation is applied to every 206 * repeat. index1=repeat, index2=times. 207 * @param division number of parts that this axis divides the structure in 208 * 209 * @throws IllegalArgumentException if the repeat relation is in the 210 * wrong format: should be double List of equal sizes. 211 * @deprecated Use {@link #addAxis(Matrix4d, int, SymmetryType)} instead. 212 * Repeats and Superposition are now inferred automatically. 213 */ 214 @Deprecated 215 public void addAxis(Matrix4d axis, List<List<Integer>> superposition, 216 List<Integer> repeats, Integer division) { 217 218 //Check correct format of repeat relations 219 if (superposition.size() != 2){ 220 throw new IllegalArgumentException( 221 "Wrong superposition format: should be double List."); 222 } else if (superposition.get(0).size() != superposition.get(1).size()){ 223 throw new IllegalArgumentException( 224 "Wrong superposition format: not equal List sizes."); 225 } 226 // Now ignores superposition & repeats except to guess symmetry type 227 SymmetryType type; 228 // Closed if superposition has a circular permutation 229 List<Integer> superPos1 = superposition.get(1); 230 if(superPos1.get(0) > superPos1.get(superPos1.size()-1)) { 231 type = SymmetryType.CLOSED; 232 } else { 233 type = SymmetryType.OPEN; 234 } 235 this.addAxis(axis,division,type); 236 } 237 /** 238 * Adds a new axis of symmetry to the bottom level of the tree 239 * 240 * @param axis the new axis of symmetry found 241 * @param order number of parts that this axis divides the structure in 242 * @param type indicates whether the axis has OPEN or CLOSED symmetry 243 */ 244 public void addAxis(Matrix4d axis, int order, SymmetryType type) { 245 axes.add(new Axis(axis,order,type,axes.size(),0)); 246 } 247 248 /** 249 * Return a list giving the number of times each axis must be applied 250 * to generate the given repeat. 251 * <P> 252 * For instance, for a D3 case <tt>getAxisCounts(4)</tt> would return [2,0], 253 * indicating that repeat 4 is generated by two applications of the 3-fold 254 * axis followed by 0 applications of the two-fold axis. 255 * 256 * @param repeat Index of the desired repeat 257 * @return array of the same length as axes giving the number of times 258 * to apply each axis. 259 */ 260 private int[] getAxisCounts(int repeat) { 261 int[] counts = new int[getNumLevels()]; 262 263 for(int i = counts.length-1; i >= 0; i--) { 264 int d = axes.get(i).getOrder(); 265 counts[i] = repeat % d; 266 repeat /= d; 267 } 268 assert repeat == 0 : "Invalid repeat index"; 269 return counts; 270 } 271 272// /** 273// * Inverse of {@link #getAxisCounts(int)}; Calculates the repeat for a 274// * particular number of applications of each axis 275// * @param counts Number of times to apply each axis 276// * @return Repeat index 277// */ 278// private int getRepeatIndex(int[] counts) { 279// int repeat = 0; 280// for(int i = 0; i< counts.length; i++) { 281// repeat += counts[i]*axes.get(i).getOrder(); 282// } 283// return repeat; 284// } 285 /** 286 * Updates an axis of symmetry, after the superposition changed. 287 * 288 * @param index old axis index 289 * @param newAxis 290 */ 291 public void updateAxis(Integer index, Matrix4d newAxis){ 292 axes.get(index).setOperator(newAxis); 293 } 294 295 /** 296 * Return the operator for all elementary axes of symmetry of the structure, that is, 297 * the axes stored in the List as unique and from which all the symmetry 298 * axes are constructed. 299 * 300 * @return axes elementary axes of symmetry. 301 */ 302 public List<Matrix4d> getElementaryAxes(){ 303 List<Matrix4d> ops = new ArrayList<Matrix4d>(getNumLevels()); 304 for(Axis axis : axes) { 305 ops.add(axis.getOperator()); 306 } 307 return ops; 308 } 309 310 /** 311 * Return all elementary axes of symmetry of the structure, that is, 312 * the axes stored in the List as unique and from which all the symmetry 313 * axes are constructed. 314 * 315 * @return axes elementary axes of symmetry. 316 */ 317 public List<Axis> getElementaryAxesObjects() { 318 return axes; 319 } 320 321 /** 322 * Get the indices of participating repeats in Cauchy two-line form. 323 * <p> 324 * Returns two lists of the same length. 325 * The first gives a list of all repeat indices which are aligned 326 * at the specified level of symmetry (e.g. 0 through the degree of this level). 327 * The second list gives the corresponding repeats after applying the 328 * operator once. 329 * 330 * @param level the axis index 331 * @return the double List of repeat relations, or null if the 332 * level is invalid 333 * @see #getRepeatsCyclicForm(int, int) for an equivalent specification with half the memory 334 */ 335 public List<List<Integer>> getRepeatRelation(int level){ 336 return getRepeatRelation(level,0); 337 } 338 339 public List<List<Integer>> getRepeatRelation(Axis axis){ 340 return getRepeatRelation(axis.getLevel(),axis.getFirstRepeat()); 341 } 342 343 public List<List<Integer>> getRepeatRelation(int level, int firstRepeat) { 344 Axis axis = axes.get(level); 345 int m = getNumRepeats(level+1);//size of the children 346 int d = axis.getOrder(); // degree of this node 347 int n = m*d; // number of repeats included 348 if(firstRepeat % n != 0) 349 throw new IllegalArgumentException(String.format("Repeat %d cannot start a block at level %s of this tree",firstRepeat,level)); 350 if(axis.getSymmType() == SymmetryType.OPEN) { 351 n -= m; // leave off last child for open symm 352 } 353 List<Integer> repeats = new ArrayList<>(n); 354 List<Integer> equiv = new ArrayList<>(n); 355 for(int i=0;i<n;i++) { 356 repeats.add(i+firstRepeat); 357 equiv.add( (i+m)%(m*d)+firstRepeat ); 358 } 359 return Arrays.asList(repeats,equiv); 360 361 } 362 363 /** 364 * Get the indices of participating repeats in cyclic form. 365 * <p> 366 * Each inner list gives a set of equivalent repeats and should have length 367 * equal to the order of the axis' operator. 368 * @param level 369 * @param firstRepeat 370 * @return 371 */ 372 public List<List<Integer>> getRepeatsCyclicForm(int level, int firstRepeat) { 373 Axis axis = axes.get(level); 374 int m = getNumRepeats(level+1);//size of the children 375 int d = axis.getOrder(); // degree of this node 376 int n = m*d; // number of repeats included 377 if(firstRepeat % n != 0) { 378 throw new IllegalArgumentException(String.format("Repeat %d cannot start a block at level %s of this tree",firstRepeat,level)); 379 } 380 if(axis.getSymmType() == SymmetryType.OPEN) { 381 n -= m; // leave off last child for open symm 382 } 383 384 List<List<Integer>> repeats = new ArrayList<>(m); 385 for(int i=0;i<m;i++) { 386 List<Integer> cycle = new ArrayList<>(d); 387 for(int j=0;j<d;j++) { 388 cycle.add(firstRepeat+i+j*m); 389 } 390 repeats.add(cycle); 391 } 392 return repeats; 393 } 394 public List<List<Integer>> getRepeatsCyclicForm(Axis axis) { 395 return getRepeatsCyclicForm(axis.getLevel(),axis.getFirstRepeat()); 396 } 397 public List<List<Integer>> getRepeatsCyclicForm(int level) { 398 return getRepeatsCyclicForm(level,0); 399 } 400 public String getRepeatsCyclicForm(Axis axis, List<?> repeats) { 401 if(repeats.size() != getNumRepeats()) { 402 throw new IllegalArgumentException("Mismatch in the number of repeats"); 403 } 404 return getRepeatsCyclicForm(getRepeatsCyclicForm(axis), repeats); 405 } 406 public static String getRepeatsCyclicForm(List<List<Integer>> cycleForm, List<?> repeats) { 407 StringBuilder str = new StringBuilder(); 408 for(List<Integer> cycle : cycleForm) { 409 str.append("("); 410 Iterator<Integer> cycleIt = cycle.iterator(); 411 str.append(repeats.get(cycleIt.next())); //should be at least one 412 while(cycleIt.hasNext()) { 413 str.append(";") 414 .append(repeats.get( cycleIt.next() )); 415 } 416 str.append(")"); 417 } 418 return str.toString(); 419 } 420 421 /** 422 * Return the transformation that needs to be applied to a 423 * repeat in order to superimpose onto repeat 0. 424 * 425 * @param repeat the repeat index 426 * @return transformation matrix for the repeat 427 */ 428 public Matrix4d getRepeatTransform(int repeat){ 429 430 Matrix4d transform = new Matrix4d(); 431 transform.setIdentity(); 432 433 int[] counts = getAxisCounts(repeat); 434 435 for(int t = counts.length-1; t>=0; t--) { 436 if( counts[t] == 0 ) 437 continue; 438 Matrix4d axis = new Matrix4d(axes.get(t).getOperator()); 439 for(int i=0;i<counts[t];i++) { 440 transform.mul(axis); 441 } 442 } 443 return transform; 444 } 445 446 /** 447 * Return the transformation that needs to be applied to 448 * repeat x in order to superimpose onto repeat y. 449 * 450 * @param x the first repeat index (transformed) 451 * @param y the second repeat index (fixed) 452 * @return transformation matrix for the repeat x 453 */ 454 public Matrix4d getRepeatTransform(int x, int y){ 455 456 Matrix4d transform = new Matrix4d(); 457 transform.setIdentity(); 458 459 int[] iCounts = getAxisCounts(x); 460 int[] jCounts = getAxisCounts(y); 461 462 int[] counts = new int[iCounts.length]; 463 for (int k = 0; k < iCounts.length; k++) 464 counts[k] = iCounts[k] - jCounts[k]; 465 466 for(int t = counts.length-1; t>=0; t--) { 467 if(counts[t] == 0) 468 continue; 469 if (counts[t] > 0) { 470 Matrix4d axis = new Matrix4d(axes.get(t).getOperator()); 471 for(int i=0;i<counts[t];i++) 472 transform.mul(axis); 473 } else if (counts[t] < 0) { 474 Matrix4d axis = new Matrix4d(axes.get(t).getOperator()); 475 axis.invert(); 476 for(int i=0;i<counts[t];i++) 477 transform.mul(axis); 478 } 479 } 480 return transform; 481 } 482 483 /** 484 * Return all symmetry axes of of the structure: the set of axes that 485 * describe all parts of the structure. This combines the elementary 486 * axes to generate all possible axes. The axes are returned in the repeat 487 * degrees. 488 * @return axes all symmetry axes of the structure. 489 */ 490 public List<Axis> getSymmetryAxes(){ 491 492 List<Axis> symmAxes = new ArrayList<>(); 493 494 Matrix4d prior = new Matrix4d(); 495 prior.setIdentity(); 496 497 getSymmetryAxes(symmAxes,prior,0,0); 498 499 500 return symmAxes; 501 } 502 /** 503 * Recursive helper 504 * @param symmAxes output list 505 * @param prior transformation aligning the first repeat of this axis with the first overall 506 * @param level current level 507 */ 508 private void getSymmetryAxes(List<Axis> symmAxes, Matrix4d prior, int level, int firstRepeat) { 509 if(level >= getNumLevels() ) { 510 return; 511 } 512 513 Axis elem = axes.get(level); 514 Matrix4d elemOp = elem.getOperator(); 515 516 // Current axis: 517 // elementary maps B -> A 518 // prior maps I -> A and J -> B 519 // want J -> I = J -> B -> A <- I= inv(prior) * elementary * prior 520 Matrix4d currAxisOp = new Matrix4d(prior); 521 currAxisOp.invert(); 522 currAxisOp.mul(elemOp); 523 currAxisOp.mul(prior); 524 Axis currAxis = new Axis(currAxisOp,elem.getOrder(),elem.getSymmType(),level,firstRepeat); 525 symmAxes.add(currAxis); 526 527 //Remember that all degrees are at least 2 528 getSymmetryAxes(symmAxes,prior,level+1,firstRepeat); 529 //New prior is elementary^d*prior 530 Matrix4d newPrior = new Matrix4d(elemOp); 531 newPrior.mul(prior); 532 int childSize = getNumRepeats(level+1); 533 getSymmetryAxes(symmAxes,newPrior,level+1,firstRepeat+childSize); 534 for(int d=2;d<elem.getOrder();d++) { 535 newPrior.mul(elemOp,newPrior); 536 getSymmetryAxes(symmAxes,newPrior,level+1,firstRepeat+childSize*d); 537 } 538 } 539 540 541// public Matrix4d getSymmetryAxis(int level, int axisNum) { 542// if(level == 0) { 543// if( axisNum != 0 ) 544// throw new IndexOutOfBoundsException("Axis number out of bounds"); 545// return axes.get(0); 546// } else { 547// if( axisNum >= degrees.get(level-1) ) 548// throw new IndexOutOfBoundsException("Axis number out of bounds"); 549// // Convert axisNum into a count of 550// 551// } 552 /** 553 * Get the number of repeats. This is equal to the product of all degrees. 554 * @return Number of repeats (leaves of the tree). 555 */ 556 public int getNumRepeats() { 557 return getNumRepeats(0); 558 } 559 560 /** 561 * Get the number of leaves from a node at the specified level. This is 562 * equal to the product of all degrees at or below the level. 563 * @param level level of the tree to cut at 564 * @return Number of repeats (leaves of the tree). 565 */ 566 private int getNumRepeats(int level) { 567 int size = 1; 568 // Return 1 for illegally high level 569 if(level < getNumLevels()) { 570 for(Axis axis : axes.subList(level, getNumLevels())) { 571 size *= axis.getOrder(); 572 } 573 } 574 return size; 575 } 576 577 /** 578 * Get the first repeat index of each axis of a specified level. 579 * @param level level of the tree to cut at 580 * @return List of first Repeats of each index, sorted in ascending order 581 */ 582 public List<Integer> getFirstRepeats(int level) { 583 List<Integer> firstRepeats = new ArrayList<Integer>(); 584 int m = getNumRepeats(level+1); //size of the level 585 int d = axes.get(level).getOrder(); //degree of this level 586 int n = m*d; // number of repeats included in each axis 587 for (int firstRepeat = 0; firstRepeat < getNumRepeats(); firstRepeat+=n) 588 firstRepeats.add(firstRepeat); 589 return firstRepeats; 590 } 591 592 public Axis getElementaryAxis(int level) { 593 return axes.get(level); 594 } 595 596 public int getNumLevels() { 597 return axes.size(); 598 } 599 600}