/*
 *                  BioJava development code
 *
 * This code may be freely distributed and modified under the
 * terms of the GNU Lesser General Public Licence.  This should
 * be distributed with the code.  If you do not have a copy,
 * see:
 *
 *      http://www.gnu.org/copyleft/lesser.html
 *
 * Copyright for this code is held jointly by the individual
 * authors.  These should be listed in @author doc comments.
 *
 * For more information on the BioJava project and its aims,
 * or to join the biojava-l mailing list, visit the home page
 * at:
 *
 *      http://www.biojava.org/
 *
 * Created on Jan 28, 2006
 *
 */
package org.biojava.nbio.structure.align.pairwise;


import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.StrucAligParameters;
import org.biojava.nbio.structure.align.helper.AligMatEl;
import org.biojava.nbio.structure.align.helper.IndexPair;
import org.biojava.nbio.structure.align.helper.JointFragments;
import org.biojava.nbio.structure.geometry.Matrices;
import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import java.io.Serializable;
import java.text.DecimalFormat;
import java.util.List;


import javax.vecmath.Matrix4d;

/**
 * Implements a class which handles one possible (alternative) solution.
 *
 * Alternative alignments arise from different seed
 * alignments or seed FPairs. The AltAlg class contains methods
 * for refinement (Dynamic Programming based) and filtering
 * (i.e. removing probably wrongly matched APairs). In the refinement
 * phase, different seed alignments can converge to the same solution.
 *
 * @author Andreas Prlic,
 * @author Peter Lackner (original Python and C code)
 * @since 3:04:26 PM
 * @version %I% %G%
 */
public class AlternativeAlignment implements Serializable{


        private static final long serialVersionUID = -6226717654562221241L;

        private int[] idx1;
        private int[] idx2;
        private String[] pdbresnum1;
        private String[] pdbresnum2;
        //short[] alig1;
        //short[] alig2;

        private int nfrags;
        private Atom center;
        private Matrix rot;
        private Atom tr;


        // the scores...
        private int gaps0;
        private int eqr0;
        private int rms0;
        private int joined;
        private int percId;
        private int cluster;
        private float score;
        private IndexPair[] aligpath;
        private int fromia;
        private Matrix currentRotMatrix;
        private Atom currentTranMatrix;

        private double rms;

        private Matrix distanceMatrix;

        public static final Logger logger =  LoggerFactory.getLogger(AlternativeAlignment.class);


        public AlternativeAlignment() {
                super();

                nfrags = 0;
                aligpath = new IndexPair[0];
                //alig1 = new short[0];
                //alig2 = new short[0];
                idx1 = new int[0];
                idx2 = new int[0];

                center = new AtomImpl();
                rot = null;
                tr = new AtomImpl();
                eqr0 = -99;
                rms0 = 99;
                joined = 0;
                gaps0 = -99;
                fromia = 0;

                currentRotMatrix = new Matrix(0,0);
                currentTranMatrix = new AtomImpl();

                distanceMatrix = new Matrix(0,0);
        }




        /** print the idx positions of this alignment
         *
         * @return a String representation
         */
        @Override
        public String toString(){
                DecimalFormat d2 = new DecimalFormat();
                // the result can be localized. To change this and enforce UK local do...
                //(DecimalFormat)NumberFormat.getInstance(java.util.Locale.UK);
                d2.setMaximumIntegerDigits(3);
                d2.setMinimumFractionDigits(2);
                d2.setMaximumFractionDigits(2);
                StringBuffer s = new StringBuffer();
                s.append("#" + getAltAligNumber() +
                                " cluster:" + cluster +
                                " eqr:" + getEqr() +
                                " rmsd:" + d2.format(getRmsd()) +
                                " %id:" + getPercId() +
                                " gaps:" + getGaps() +
                                " score:" + d2.format(score)    );

                return s.toString();
        }


        /** get the number of the cluster this alignment belongs to
         *
         * @return an int giving the number of the cluster
         */
        public int getCluster() {
                return cluster;
        }



        /** set the number of the cluster this alignment belongs to.
         * All alignments in a cluster are quite similar.
         *
         * @param cluster the number of the cluster
         */
        public void setCluster(int cluster) {
                this.cluster = cluster;
        }




        public double getRmsd() {
                return rms;
        }

        /** the rms in the structurally equivalent regions
         *
         * @param rms
         */
        public void setRms(double rms) {
                this.rms = rms;
        }

        /** the alignment score
         *
         * @return the score of this alignment
         */
        public float getScore() {
                return score;
        }

        public void setScore(float score) {
                this.score = score;
        }


        /** return the number of gaps in this alignment
         *
         * @return the number of Gaps
         */
        public int getGaps(){
                return gaps0;
        }

        /** returns the number of euqivalent residues in this alignment
         *
         * @return the number of equivalent residues
         */
        public int getEqr(){
                return eqr0 ;
        }

        /** the positions of the structure equivalent positions in atom set 1
         *
         * @return the array of the positions
         */
        public int[] getIdx1(){
                return idx1;
        }

        /** the positions of the structure equivalent atoms in atom set 2
         *
         * @return the array of the positions
         */
        public int[] getIdx2(){
                return idx2;
        }

        public int getPercId() {
                return percId;
        }

        public void setPercId(int percId) {
                this.percId = percId;
        }

        /** Set apairs according to a seed position.
         *
         * @param l
         * @param i
         * @param j
         */
        public void  apairs_from_seed(int l,int i, int j){
                aligpath = new IndexPair[l];
                idx1 = new int[l];
                idx2 = new int[l];
                for (int x = 0 ; x < l ; x++) {
                        idx1[x]=i+x;
                        idx2[x]=j+x;
                        aligpath[x] = new IndexPair((short)(i+x),(short)(j+x));
                }
        }

        /** Set apairs according to a list of (i,j) tuples.
         *
         * @param jf a JoingFragment
         */
        public void apairs_from_idxlst(JointFragments jf) {
                List<int[]> il = jf.getIdxlist();
                //System.out.println("Alt Alig apairs_from_idxlst");

                aligpath = new IndexPair[il.size()];
                idx1 = new int[il.size()];
                idx2 = new int[il.size()];
                for (int i =0 ; i < il.size();i++) {
                        int[] p = il.get(i);
                        //System.out.print(" idx1 " + p[0] + " idx2 " + p[1]);
                        idx1[i] = p[0];
                        idx2[i] = p[1];
                        aligpath[i] = new IndexPair((short)p[0],(short)p[1]);

                }
                eqr0  = idx1.length;
                //System.out.println("eqr " + eqr0);
                gaps0 = count_gaps(idx1,idx2);

        }


        /** returns the sequential number of this alternative alignment
         *
         * @return the sequential number of this alternative alignment
         */
        public int getAltAligNumber() {
                return fromia;
        }

        public void setAltAligNumber(int fromia) {
                this.fromia = fromia;
        }

        /** rotate and shift atoms with currentRotMatrix and current Tranmatrix
         *
         * @param ca
         */
        private void rotateShiftAtoms(Atom[] ca){

                for (int i  = 0 ; i < ca.length; i++){
                        Atom c = ca[i];

                        Calc.rotate(c,currentRotMatrix);
                        Calc.shift(c,currentTranMatrix);
                        //System.out.println("after " + c);
                        ca[i] = c;
                }
                //System.out.println("after " + ca[0]);
        }

        public void finish(StrucAligParameters params,Atom[]ca1,Atom[]ca2) throws StructureException{

                Atom[] ca3 = new Atom[ca2.length];
                for ( int i = 0 ; i < ca2.length;i++){
                        ca3[i] = (Atom) ca2[i].clone();

                }
                // do the inital superpos...

                super_pos_alig(ca1,ca3,idx1,idx2,true);
                rotateShiftAtoms(ca3);

                calcScores(ca1,ca2);
                logger.debug("eqr " + eqr0 + " " + gaps0 + " "  +idx1[0] + " " +idx1[1]);

                getPdbRegions(ca1,ca2);

        }

        public static Matrix getDistanceMatrix(Atom[] ca1, Atom[]ca2){

                int r = ca1.length;
                int c = ca2.length;

                Matrix out = new Matrix(r,c);

                for (int i=0; i<r; i++) {
                        Atom a1 = ca1[i];
                        for (int j=0;j<c;j++){
                                Atom b1 = ca2[j];


                                double d = Calc.getDistance(a1,b1);
                                out.set(i,j,d);

                        }
                }
                return out;
        }

        private Alignable getInitalStrCompAlignment(
                        Atom[] ca1,
                        Atom[]ca2,
                        StrucAligParameters params
                        ){

                int rows = ca1.length;
                int cols = ca2.length;

                float gapOpen = params.getGapOpen();
                float gapExtension = params.getGapExtension();
                float co = params.getCreate_co();

                Alignable al = new StrCompAlignment(rows,cols);
                al.setGapExtCol(gapExtension);
                al.setGapExtRow(gapExtension);
                al.setGapOpenCol(gapOpen);
                al.setGapOpenRow(gapOpen);

                AligMatEl[][] aligmat = al.getAligMat();

                //System.out.println(rows + " " + cols );

                aligmat[0][cols]    = new AligMatEl();
                aligmat[rows][0]    = new AligMatEl();
                aligmat[rows][cols] = new AligMatEl();

                for (int i = 0; i < rows; i++) {
                        Atom a1 = ca1[i];

                        aligmat[i][0]    = new AligMatEl();
                        //aligmat[i][cols] = new AligMatEl();

                        for (int j = 0; j < cols; j++) {
                                //System.out.println("index " + i + " " + j);
                                aligmat[0][j] = new AligMatEl();

                                Atom b1 = ca2[j];
                                double d = 999;

                                d = Calc.getDistance(a1,b1);


                                AligMatEl e = new AligMatEl();
                                if (d > co) {
                                        e.setValue(0);
                                } else {
                                        double s =  2 * co / ( 1 + ( d/co) * (d/co)) - co;
                                        e.setValue((int) Math.round(Gotoh.ALIGFACTOR * s));
                                }
                                aligmat[i+1][j+1] = e;
                        }
                }

                return al;
        }

        /*private Alignable getAlignableFromSimMatrix (Matrix sim,StrucAligParameters params){
                //System.out.println("align_NPE");

                float gapOpen = params.getGapOpen();
                float gapExtension = params.getGapExtension();

                int rows = sim.getRowDimension();
                int cols = sim.getColumnDimension();

                Alignable al = new StrCompAlignment(rows,cols);
                al.setGapExtCol(gapExtension);
                al.setGapExtRow(gapExtension);
                al.setGapOpenCol(gapOpen);
                al.setGapOpenRow(gapOpen);
                //System.out.println("size of aligmat: " + rows+1 + " " + cols+1);
                //AligMatEl[][] aligmat = new AligMatEl[rows+1][cols+1];
                AligMatEl[][] aligmat = al.getAligMat();

                for (int i = 0; i < rows; i++) {
                        for (int j = 0; j < cols; j++) {

                                int e=0;
                                //if ( ( i < rows) &&
                                //        ( j < cols)) {
                                //TODO: the ALIGFACTOR calc should be hidden in Gotoh!!

                                e = (int)Math.round(Gotoh.ALIGFACTOR * sim.get(i,j));
                                //}
                                //System.out.println(e);
                                //AligMatEl am = new AligMatEl();
                                aligmat[i+1][j+1].setValue(e);
                                //.setValue(e);
                                //am.setTrack(new IndexPair((short)-99,(short)-99));
                                //igmat[i+1][j+1] = am;

                        }
                }
                //al.setAligMat(aligmat);


                return al;
        }*/
        /** Refinement procedure based on superposition and dynamic programming.
         * Performs an iterative refinement. Several methods apply such a procedure,
         * e.g. CE or ProSup. Here we additionally test for circular permutation,
         * which are in the same frame of superposition as the optimal alignment.
         * This feature may be switched off by setting permsize to -1.
         *
         * @param params the parameters
         * @param ca1 atoms of structure 1
         * @param ca2 atoms of structure 2
         * @throws StructureException
         */

        public void refine(StrucAligParameters params,Atom[]ca1,Atom[]ca2) throws StructureException{
                // System.out.println("refine Alternative Alignment #"+ getAltAligNumber()+" l1:" + ca1.length + " l2:" + ca2.length);
                //              for ( int i= 0 ; i < idx1.length;i++){
                //              System.out.println(idx1[i] + " " + idx2[i]);
                //              }

                // avoid any operations on the original Atoms ...
                Atom[] ca3 = new Atom[ca2.length];
                for ( int i = 0 ; i < ca2.length;i++){
                        ca3[i] = (Atom) ca2[i].clone();

                }
                //Atom[] ca3 = (Atom[]) ca2.clone();

                // do the inital superpos...

                super_pos_alig(ca1,ca3,idx1,idx2,false);

                //JmolDisplay jmol1 = new JmolDisplay(ca1,ca2,"after first alig " );
                //jmol1.selectEquiv(idx1,idx2);
                //int[] jmoltmp1 = idx1;
                //int[] jmoltmp2 = idx2;

                int lenalt =idx1.length;
                int lenneu = aligpath.length;
                //Matrix rmsalt = currentRotMatrix;



                int ml = Math.max(ca1.length, ca3.length);
                idx1 = new int[ml];
                idx2 = new int[ml];

                //int m1l = ca1.length;
                //int m2l = ca3.length;

                //Matrix sim = new Matrix(ca1.length,ca3.length,0);

                int maxiter = params.getMaxIter();
                for (int iter = 0 ; iter< maxiter; iter++){

                        float  subscore = 0.0f;

                        rotateShiftAtoms(ca3);

                        // do the dynamic alignment...
                        Alignable ali = getInitalStrCompAlignment(ca1,ca3, params);

                        new Gotoh(ali);

                        this.score = ali.getScore();
                        //System.out.println("score " + score);
                        IndexPair[] path = ali.getPath();

                        int pathsize = ali.getPathSize();


                        // now for a superimposable permutation. First we need to check the size and
                        // position of the quadrant left out by the optimal path
                        int firsta = path[0].getRow();
                        int firstb = path[0].getCol();
                        int lasta  = path[pathsize-1].getRow();
                        int lastb  = path[pathsize-1].getCol();

                        int quada_beg, quada_end, quadb_beg,quadb_end;

                        if ( firsta == 0){
                                quada_beg = lasta +1;
                                quada_end = ca1.length -1;
                                quadb_beg = 0;
                                quadb_end = firstb-1;
                        } else {
                                quada_beg = 0;
                                quada_end = firsta-1;
                                quadb_beg = lastb+1;
                                quadb_end = ca3.length -1;
                        }

                        // if the unaligned quadrant is larger than permsize
                        int permsize = params.getPermutationSize();

                        if ( permsize > 0 &&
                                        (quada_end - quada_beg >= permsize) &&
                                        ( quadb_end - quadb_beg >= permsize)) {
                                AligMatEl[][] aligmat = ali.getAligMat();
                                Matrix submat = new Matrix(quada_end-quada_beg, quadb_end - quadb_beg) ;
                                //then we copy the score values of the quadrant into a new matrix

                                //System.out.println(quada_beg + " " + quada_end + " " +quadb_beg+" " + quadb_end);
                                //System.out.println(sim.getRowDimension() + " " + sim.getColumnDimension());
                                Atom[] tmp1 = new Atom[quada_end - quada_beg ];
                                Atom[] tmp2 = new Atom[quadb_end - quadb_beg ];

                                for ( int s = quada_beg; s < quada_end; s++){
                                        //System.out.println(s + " " +( quada_end-s));
                                        tmp1[quada_end - s -1] = ca1[s];
                                        for ( int t = quadb_beg; t < quadb_end; t++){
                                                if ( s == quada_beg )
                                                        tmp2[quadb_end - t -1 ] = ca3[t];

                                                //System.out.println(s+" " +t);
                                                //double val = sim.get(s,t);
                                                double val = aligmat[s][t].getValue();
                                                submat.set(s-quada_beg,t-quadb_beg,val);
                                        }
                                }
                                // and perform a dp alignment again. (Note, that we fix the superposition).

                                //Alignable subali = getAlignableFromSimMatrix(submat,params);
                                Alignable subali = getInitalStrCompAlignment(tmp1, tmp2, params);
                                subscore = subali.getScore();
                                //System.out.println("join score" + score + " " + subscore);
                                this.score = score+subscore;

                                IndexPair[] subpath = subali.getPath();
                                int subpathsize = subali.getPathSize();
                                for ( int p=0;p<subpath.length;p++){
                                        IndexPair sp = subpath[p];
                                        sp.setRow((short)(sp.getRow()+quada_beg));
                                        sp.setCol((short)(sp.getCol()+quadb_beg));
                                }

                                // now we join the two solutions
                                if ( subpathsize > permsize){
                                        IndexPair[] wholepath = new IndexPair[pathsize+subpathsize];
                                        for ( int t=0; t < pathsize; t++){
                                                wholepath[t] = path[t];
                                        }
                                        for ( int t=0 ; t < subpathsize; t++){
                                                wholepath[t+pathsize] = subpath[t];
                                        }
                                        pathsize += subpathsize;
                                        path = wholepath;
                                        ali.setPath(path);
                                        ali.setPathSize(pathsize);
                                }

                        }


                        int j =0 ;
                        //System.out.println("pathsize,path " + pathsize + " " + path.length);

                        for (int i=0; i < pathsize; i++){
                                int x = path[i].getRow();
                                int y = path[i].getCol();
                                //System.out.println(x+ " " + y);
                                double d = Calc.getDistance(ca1[x], ca3[y]);
                                //double d = dismat.get(x,y);
                                // now we apply the evaluation distance cutoff
                                if ( d < params.getEvalCutoff()){
                                        //System.out.println("j:"+j+ " " + x+" " + y + " " + d );

                                        idx1[j] = x;
                                        idx2[j] = y;
                                        j+=1;
                                }
                        }

                        lenneu = j;
                        //System.out.println("lenalt,neu:" + lenalt + " " + lenneu);


                        int[] tmpidx1 = new int[j];
                        int[] tmpidx2 = new int[j];
                        //String idx1str ="idx1: ";
                        //String idx2str ="idx2: ";
                        for (int i = 0 ; i<j;i++){
                                //idx1str += idx1[i]+ " ";
                                //idx2str += idx2[i]+ " ";
                                tmpidx1[i] = idx1[i];
                                tmpidx2[i] = idx2[i];
                        }
                        //System.out.println(idx1str);
                        //System.out.println(idx2str);
                        //jmol.selectEquiv(idx1,idx2);

                        //if ( j > 0)
                        //  System.out.println("do new superimpos " + idx1.length + " " + idx2.length + " " + idx1[0] + " " + idx2[0]);
                        super_pos_alig(ca1,ca3,tmpidx1,tmpidx2,false);

                        this.aligpath = path;

                        if (lenneu == lenalt)
                                break;

                }



                idx1 = (int[]) FragmentJoiner.resizeArray(idx1,lenneu);
                idx2 = (int[]) FragmentJoiner.resizeArray(idx2,lenneu);
                //              new ... get rms...
                // now use the original atoms to get the rotmat relative to the original structure...
                super_pos_alig(ca1,ca2,idx1,idx2,true);
                eqr0 = idx1.length;
                gaps0 = count_gaps(idx1,idx2);

                getPdbRegions(ca1,ca2);
                //System.out.println("eqr " + eqr0 + " aligpath len:"+aligpath.length+ " gaps:" + gaps0 + " score " + score);
        }

        private void getPdbRegions(Atom[] ca1, Atom[] ca2){
                pdbresnum1 = new String[idx1.length];
                pdbresnum2 = new String[idx2.length];

                for (int i =0 ; i < idx1.length;i++){
                        Atom a1 = ca1[idx1[i]];
                        Atom a2 = ca2[idx2[i]];
                        Group p1 = a1.getGroup();
                        Group p2 = a2.getGroup();
                        Chain c1 = p1.getChain();
                        Chain c2 = p2.getChain();

                        String cid1 = c1.getId();
                        String cid2 = c2.getId();

                        String pdb1 = p1.getResidueNumber().toString();
                        String pdb2 = p2.getResidueNumber().toString();


                        if ( ! cid1.equals(" "))
                                pdb1 += ":" + cid1;


                        if ( ! cid2.equals(" "))
                                pdb2 += ":" + cid2;


                        pdbresnum1[i] = pdb1;
                        pdbresnum2[i] = pdb2;
                }
        }


        public String[] getPDBresnum1() {
                return pdbresnum1;
        }




        public void setPDBresnum1(String[] pdbresnum1) {
                this.pdbresnum1 = pdbresnum1;
        }




        public String[] getPDBresnum2() {
                return pdbresnum2;
        }




        public void setPDBresnum2(String[] pdbresnum2) {
                this.pdbresnum2 = pdbresnum2;
        }




        /**
         * Count the number of gaps in an alignment represented by idx1,idx2.
         *
         * @param i1
         * @param i2
         * @return the number of gaps in this alignment
         */
        private int count_gaps(int[] i1, int[] i2){

                int i0 = i1[0];
                int j0 = i2[0];
                int gaps = 0;
                for (int i =1 ; i<i1.length;i++ ){
                        if ( Math.abs(i1[i]-i0) != 1 ||
                                        ( Math.abs(i2[i]-j0) != 1)){
                                gaps +=1;
                        }
                        i0 = i1[i];
                        j0 = i2[i];
                }

                return gaps;
        }



        public void calculateSuperpositionByIdx(Atom[] ca1, Atom[] ca2)throws StructureException {

                super_pos_alig(ca1,ca2,idx1,idx2,false);

        }

        /** Superimposes two molecules according to residue index list idx1 and idx2.
         *  Does not change the original coordinates.
         *  as an internal result the rotation matrix and shift vectors for are set
         *
         * @param ca1 Atom set 1
         * @param ca2 Atom set 2
         * @param idx1 idx positions in set1
         * @param idx2 idx positions in set2
         * @param getRMS a flag if the RMS should be calculated
         * @throws StructureException
         */



        private void super_pos_alig(Atom[]ca1,Atom[]ca2,int[] idx1, int[] idx2, boolean getRMS) throws StructureException{

                //System.out.println("superpos alig ");
                Atom[] ca1subset = new Atom[idx1.length];
                Atom[] ca2subset = new Atom[idx2.length];

                for (int i = 0 ; i < idx1.length;i++){
                        //System.out.println("idx1 "+ idx1[i] + " " + idx2[i]);
                        int pos1 =  idx1[i];
                        int pos2 =  idx2[i];

                        ca1subset[i] =  ca1[pos1];
                        ca2subset[i] = (Atom) ca2[pos2].clone();
                }

                Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1subset),
                                Calc.atomsToPoints(ca2subset));
                this.currentRotMatrix  = Matrices.getRotationJAMA(trans);
                this.currentTranMatrix = Calc.getTranslationVector(trans);
                //currentRotMatrix.print(3,3);
                if ( getRMS) {
                        rotateShiftAtoms(ca2subset);
                        this.rms = Calc.rmsd(ca1subset,ca2subset);
                }


        }


        /** returns the rotation matrix that needs to be applied to structure 2 to rotate on structure 1
         *
         * @return the rotation Matrix
         */
        public Matrix getRotationMatrix(){
                return currentRotMatrix;
        }

        /** returns the shift vector that has to be applied on structure to to shift on structure one
         *
         * @return the shift vector
         */
        public Atom getShift(){
                return currentTranMatrix;
        }



        /** calculates  scores for this alignment ( %id )
         * @param ca1 set of Atoms for molecule 1
         * @param ca2 set of Atoms for molecule 2

         */
        public void calcScores(Atom[] ca1, Atom[] ca2){
                eqr0 = idx1.length;
                gaps0 = count_gaps(idx1,idx2);

                percId = 0;
                // calc the % id
                for (int i=0 ; i< idx1.length; i++){
                        Atom a1 = ca1[idx1[i]];
                        Atom a2 = ca2[idx2[i]];

                        Group g1 = a1.getGroup();
                        Group g2 = a2.getGroup();
                        if ( g1.getPDBName().equals(g2.getPDBName())){
                                percId++;
                        }

                }
        }



        /** create an artifical Structure object that contains the two
         * structures superimposed onto each other. Each structure is in a separate model.
         * Model 1 is structure 1 and Model 2 is structure 2.
         *
         * @param s1 the first structure. its coordinates will not be changed
         * @param s2 the second structure, it will be cloned and the cloned coordinates will be rotated according to the alignment results.
         * @return composite structure containing the 2 aligned structures as a models 1 and 2
         */
        public Structure getAlignedStructure(Structure s1, Structure s2){
                // do not change the original coords ..
                Structure s3 = s2.clone();

                currentRotMatrix.print(3,3);

                Calc.rotate(s3, currentRotMatrix);
                Calc.shift( s3, currentTranMatrix);

                Structure newpdb = new StructureImpl();
                newpdb.setPdbId(null);
                newpdb.setName("Aligned with BioJava");


                newpdb.addModel(s1.getChains(0));
                newpdb.addModel(s3.getChains(0));

                return newpdb;

        }

        /** converts the alignment to a PDB file
         * each of the structures will be represented as a model.
         *
         *
         * @param s1
         * @param s2
         * @return a PDB file as a String
         */
        public String toPDB(Structure s1, Structure s2){


                Structure newpdb = getAlignedStructure(s1, s2);

                return newpdb.toPDB();
        }



        /** The distance matrix this alignment is based on
         *
         * @return a Matrix object.
         */
        public Matrix getDistanceMatrix() {
                return distanceMatrix;
        }



        /** The distance matrix this alignment is based on
         *
         * @param distanceMatrix
         */
        public void setDistanceMatrix(Matrix distanceMatrix) {
                this.distanceMatrix = distanceMatrix;
        }




        public IndexPair[] getPath() {
                return aligpath;
        }



}