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}