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}