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 * Created on Sep 25, 2009
021 * Author: Andreas Prlic
022 *
023 */
024
025package org.biojava.nbio.structure.align.ce;
026
027import org.biojava.nbio.core.alignment.matrices.ScaledSubstitutionMatrix;
028import org.biojava.nbio.core.alignment.template.SubstitutionMatrix;
029import org.biojava.nbio.structure.Atom;
030import org.biojava.nbio.structure.Calc;
031import org.biojava.nbio.structure.Group;
032import org.biojava.nbio.structure.StructureException;
033import org.biojava.nbio.structure.StructureTools;
034import org.biojava.nbio.structure.align.model.AFP;
035import org.biojava.nbio.structure.align.model.AFPChain;
036import org.biojava.nbio.structure.align.util.AFPAlignmentDisplay;
037import org.biojava.nbio.structure.align.util.AFPChainScorer;
038import org.biojava.nbio.structure.geometry.Matrices;
039import org.biojava.nbio.structure.geometry.SuperPositions;
040import org.biojava.nbio.structure.jama.Matrix;
041import org.biojava.nbio.core.sequence.compound.AminoAcidCompound;
042import org.biojava.nbio.core.sequence.compound.AminoAcidCompoundSet;
043
044import java.util.ArrayList;
045import java.util.List;
046
047import javax.vecmath.Matrix4d;
048
049
050
051/** This is based on the original Combinatorial Extension (CE) source code from 2003 or 2004 (CE version 2.3),
052 * as has been originally developed by I. Shindyalov and P.Bourne (1998).
053 * The original CE paper is available from here: <a href="http://peds.oxfordjournals.org/cgi/content/short/11/9/739">http://peds.oxfordjournals.org/cgi/content/short/11/9/739</a>.
054 *
055 * This class is a pretty much exact 1:1 port from C, where I cared about exact reproduce of the CE results
056 * and not about Java style.
057 *
058 * @author Andreas Prlic
059
060 *
061 */
062public class CECalculator {
063
064        protected static final boolean isPrint = false;
065        private static final boolean debug = false;
066
067        int[] f1;
068        int[] f2;
069        double[][]dist1;
070        double[][]dist2;
071        protected double[][]mat;
072        protected int[] bestTrace1;
073        protected int[] bestTrace2;
074        protected int[][] bestTraces1;
075        protected int[][] bestTraces2;
076        protected int nBestTrace;
077        protected int nBestTraces;
078        double[] d_ = new double[20];
079        protected int[] bestTracesN;
080        protected double bestTraceScore;
081        protected int nTrace;
082        protected double[] bestTracesScores;
083        protected int[] trace1;
084        protected int[] trace2;
085
086        protected static final  double  zThr=-0.1;
087
088        long timeStart ;
089        long timeEnd;
090        private int nAtom;
091
092        // the equivalent positions in the alignment...
093        private int[] align_se1;
094        private int[] align_se2;
095
096
097        private int lcmp;
098        private int[] bestTraceLen;
099        private Matrix r;
100        private Atom t;
101        protected int nTraces;
102
103        private double z;
104        private int[] traceIndexContainer;
105
106        protected CeParameters params;
107        // SHOULD these fields be PARAMETERS?
108
109        protected static final int nIter = 1;
110        private static final boolean distAll = false;
111
112        List<MatrixListener> matrixListeners;
113
114
115
116        public CECalculator(CeParameters params){
117                timeStart = System.currentTimeMillis();
118                dist1= new double[0][0];
119                dist2= new double[0][0];
120                this.params = params;
121                matrixListeners = new ArrayList<MatrixListener>();
122
123        }
124
125        /**
126         *
127         * @param afpChain A new AFPChain, which will be filled in by this function
128         * @param ca1
129         * @param ca2
130         * @return afpChain
131         * @throws StructureException
132         */
133        public AFPChain extractFragments(AFPChain afpChain,
134                        Atom[] ca1, Atom[] ca2) throws StructureException{
135
136                int nse1 = ca1.length;
137                int nse2 = ca2.length;
138
139                afpChain.setCa1Length(nse1);
140                afpChain.setCa2Length(nse2);
141
142                int traceMaxSize=nse1<nse2?nse1:nse2;
143
144                f1 = new int[nse1];
145                f2 = new int[nse2];
146
147                dist1 = initIntraDistmatrix(ca1, nse1);
148                dist2 = initIntraDistmatrix(ca2, nse2);
149
150
151                if ( debug )
152                        System.out.println("parameters: " + params);
153
154                if ( params.getScoringStrategy() == CeParameters.ScoringStrategy.SEQUENCE_CONSERVATION){
155                        if ( params.getSeqWeight() < 1)
156                                params.setSeqWeight(2);
157                }
158
159                int winSize = params.getWinSize();
160
161                int winSizeComb1 = (winSize-1)*(winSize-2)/2;
162
163                traceIndexContainer = new int[traceMaxSize];
164
165                // CE: unused code. distAll is always false and both loops do the same???
166                // CE v2.3 calls this Weight factors for trace extension
167                if(distAll ) {
168                        for(int i=0; i<traceMaxSize; i++)
169                                traceIndexContainer[i]=(i+1)*i*winSize*winSize/2+(i+1)*winSizeComb1;
170                } else {
171                        for(int i=0; i<traceMaxSize; i++) {
172                                traceIndexContainer[i]=(i+1)*i*winSize/2+(i+1)*winSizeComb1;
173
174
175                        }
176                }
177
178                // verified: a[] is set correctly.
179
180                mat = initSumOfDistances(nse1, nse2, winSize, winSizeComb1, ca1, ca2);
181
182
183
184//              try {
185//                      Matrix m2 = new Matrix(mat).copy();
186//                      JPanel panel = GuiWrapper.getScaleableMatrixPanel(m2);
187//                      JFrame frame = new JFrame();
188//                      frame.addWindowListener(new WindowAdapter(){
189//                              public void windowClosing(WindowEvent e){
190//                                      JFrame f = (JFrame) e.getSource();
191//                                      f.setVisible(false);
192//                                      f.dispose();
193//                              }
194//                      });
195//
196//
197//                      frame.getContentPane().add(panel);
198//
199//                      frame.pack();
200//                      frame.setVisible(true);
201//              } catch (Exception e) {
202//                      e.printStackTrace();
203//              }
204
205
206                // Set the distance matrix
207                //afpChain.setDistanceMatrix(new Matrix(mat.clone()));
208
209
210                //
211                //                         double rmsdThr = params.getRmsdThr();
212                //                         StringBuffer buf = new StringBuffer("  ");
213                //                         for(int i=0; i<nse2; i++)
214                //                            buf.append(String.format("%c", i%10==0?(i%100)/10+48:32));
215                //                         buf.append("\n");
216                //                         for(int i=0; i<nse1; i++) {
217                //                            buf.append(String.format("%c ", i%10==0?(i%100)/10+48:32));
218                //                            for(int j=0; j<nse2; j++)
219                //                               buf.append(String.format("%c", (mat[i][j])<rmsdThr?'+':'X'));
220                //                            //printf("%c", ((int)*(mat[i]+j)/40)>9?'*':((int)*(mat[i]+j)/40)+48);
221                //                            buf.append("\n");
222                //                         }
223                //                         buf.append("\n");
224                //
225                //                         System.out.println(buf.toString());
226                //
227
228
229                return afpChain;
230        }
231
232        /**
233         * Evaluates the distance between two atoms
234         * Several scoring functions are implemented and can be changed by calling
235         * {@link CeParameters#setScoringStrategy(Integer) setScoringStrategy()}
236         * on {@link CeParameters parameter} object this CECalculator was created with.
237         * <p>
238         * Scoring Strategies:<dl>
239         * <dt>DEFAULT_SCORING_STRATEGY</dt>
240         * <dd>Strategy of the original CE publication; CA-CA distance</dd>
241         *
242         * <dt>SIDE_CHAIN_SCORING</dt>
243         * <dd>CB-CB distance. This performs better for sheets and helices than CA.</dd>
244         *
245         * <dt>SIDE_CHAIN_ANGLE_SCORING</dt>
246         * <dd>Use the dot product (eg the cosine) of the two CA-CB vectors.</dd>
247         *
248         * <dt>CA_AND_SIDE_CHAIN_ANGLE_SCORING</dt>
249         * <dd>Equivalent to DEFAULT_SCORING_STRATEGY + SIDE_CHAIN_ANGLE_SCORING</dd>
250         * </dl>
251         *
252         *  <dt>SEQUENCE_CONSERVATION</dt>
253         * <dd>A mix between the DEFAULT_SCORING_STRATEGY and a scoring function that favors the alignment of sequence conserved positions in the alignment</dd>
254         * </dl>
255         *
256         *
257         *
258         * @param ca1 The CA of the first residue
259         * @param ca2 The CA of the second residue
260         * @return The distance between the two fragments, according to the selected
261         * scoring strategy. Lower distances are better alignments.
262         * @throws StructureException
263         */
264        private double getDistanceWithSidechain(Atom ca1, Atom ca2) throws StructureException {
265
266                if ( params.getScoringStrategy() == CeParameters.ScoringStrategy.CA_SCORING) {
267
268                        return Calc.getDistance(ca1,ca2);
269
270                }
271
272                double dist;
273                Group g1 = ca1.getGroup();
274                Atom cb1 = null;
275                if ( g1.hasAtom(StructureTools.CB_ATOM_NAME)) {
276                        cb1 = g1.getAtom(StructureTools.CB_ATOM_NAME);
277                }
278                //
279                Group g2 = ca2.getGroup();
280                Atom cb2 = null;
281                if ( g2.hasAtom(StructureTools.CB_ATOM_NAME)) {
282                        cb2 = g2.getAtom(StructureTools.CB_ATOM_NAME);
283                }
284
285
286                if ( params.getScoringStrategy() == CeParameters.ScoringStrategy.SIDE_CHAIN_SCORING) {
287
288
289                        // here we are using side chain orientation for scoring...
290
291
292                        // score type 1    consider side chain distances
293                        if ( cb1 != null && cb2 != null) {
294                                // CB distance
295                                dist = Calc.getDistance(cb1,cb2);
296                                //dist = dist / 2.;
297                        } else {
298                                dist = Calc.getDistance(ca1,ca2);
299                        }
300
301                        return dist;
302                }
303
304                else if ( params.getScoringStrategy() == CeParameters.ScoringStrategy.SIDE_CHAIN_ANGLE_SCORING){
305
306                        // score type 2 add angle info
307
308
309                        if ( cb1 != null && cb2 != null) {
310                                // If the CA were overlaid, what is the distance between the CB?
311                                // Recall c^2 = a^2 + b^2 -2ab*cos(theta), so this is a function of angle
312                                Atom c1 = Calc.subtract(cb1, ca1);
313                                Atom c2 = Calc.subtract(cb2, ca2);
314                                Atom newA = Calc.subtract(c2, c1);
315                                dist = Calc.amount(newA);
316                        }  else {
317                                //dist += Calc.getDistance(ca1,ca2);
318                                dist = 0;
319                        }
320
321                        return dist;
322
323                }
324                else if ( params.getScoringStrategy() == CeParameters.ScoringStrategy.CA_AND_SIDE_CHAIN_ANGLE_SCORING){
325
326                        // score type 3
327                        // CA distance + cos(angle)
328                        dist = 0;
329                        if ( cb1 != null && cb2 != null) {
330                                Atom cacb1 = Calc.subtract(cb1, ca1);
331                                Atom cacb2 = Calc.subtract(cb2, ca2);
332                                Atom newA = Calc.subtract(cacb2, cacb1);
333                                //System.out.format("CACB 1: %s\nCACB 2: %s\ndiff: %s\nd: %f\n",cacb1.toString(),cacb2.toString(),newA.toString(),Calc.amount(newA));
334                                dist += Calc.amount(newA);
335                        }
336                        dist += Calc.getDistance(ca1,ca2);
337
338                        return dist;
339
340                } else if ( params.getScoringStrategy() == CeParameters.ScoringStrategy.SEQUENCE_CONSERVATION){
341                        if ( cb1 != null && cb2 != null) {
342                                // CB distance
343                                dist = Calc.getDistance(cb1,cb2);
344                                //dist = dist / 2.;
345                        } else {
346                                dist = Calc.getDistance(ca1,ca2);
347                        }
348                        return dist;
349
350
351                }
352                else {
353                        // unsupported scoring scheme
354                        return Calc.getDistance(ca1,ca2);
355                }
356        }
357
358        /** build up intramolecular distance matrix dist1 & dist2
359         *
360         * @param ca
361         * @param nse
362         * @return
363         * @throws StructureException
364         */
365        private double[][] initIntraDistmatrix(Atom[] ca, int nse) throws StructureException
366        {
367
368
369                double[][] intraDist = new double[nse][nse];
370
371                //
372                for(int ise1=0; ise1<nse; ise1++)  {
373
374                        for(int ise2=0; ise2<nse; ise2++)  {
375                                intraDist[ise1][ise2] = getDistanceWithSidechain(ca[ise1], ca[ise2]);
376
377                        }
378                }
379                return intraDist;
380        }
381
382
383        public double[][] initSumOfDistances(int nse1, int nse2, int winSize, int  winSizeComb1, Atom[] ca1, Atom[] ca2) {
384
385                double d;
386
387                double[][] mat   = new double[nse1][nse2];
388
389                // init the initial mat[] array.
390                // at this stage mat contains the sum of the distances of fragments of the matrices dist1, dist
391                for(int ise1=0; ise1<nse1; ise1++) {
392                        for(int ise2=0; ise2<nse2; ise2++) {
393
394                                mat[ise1][ise2]=-1.0;
395
396                                if(ise1>nse1-winSize || ise2>nse2-winSize) continue;
397
398                                d=0.0;
399                                // this sums up over the distances of the fragments
400                                for(int is1=0; is1<winSize-2; is1++)
401                                        for(int is2=is1+2; is2<winSize; is2++) {
402                                                //System.out.println("pos1 :" +  (ise1+is1) + " " + (ise1+is2) +  " " + (ise2+is1) + " " + (ise2+is2));
403                                                // is this abs or floor? check!
404                                                d+=Math.abs(dist1[ise1+is1][ise1+is2]-dist2[ise2+is1][ise2+is2]);
405                                        }
406                                mat[ise1][ise2]=d/winSizeComb1;
407
408                                //System.out.println("mat ["+ise1+"]["+ise2+"]="+mat[ise1][ise2]);
409                        }
410
411                }
412
413                // verified: mat[][] probably ok.
414
415                return mat;
416        }
417
418
419
420
421
422
423        @SuppressWarnings("unused")
424        public void traceFragmentMatrix( AFPChain afpChain,
425                        Atom[] ca1, Atom[] ca2) {
426
427                double rmsdThr = params.getRmsdThr();
428
429
430                double oldBestTraceScore=10000.0;
431                bestTraceScore = 100.0;
432                nBestTrace=0;
433                int nBestTrace0 = 0;
434                int winSize = params.getWinSize();
435                int winSizeComb1=(winSize-1)*(winSize-2)/2;
436                boolean distAll = false;
437
438                int winSizeComb2=distAll?winSize*winSize:winSize;
439                double rmsdThrJoin = params.getRmsdThrJoin();
440
441                double z0;
442
443                //double bestTraceZScore=-1.0;
444
445                int nse1 = ca1.length;
446                int nse2 = ca2.length;
447
448                //System.out.println("nse1 :" +nse1 + " nse2: " + nse2);
449
450                int traceMaxSize=nse1<nse2?nse1:nse2;
451
452                bestTrace1 = new int [traceMaxSize];
453                bestTrace2 = new int [traceMaxSize];
454                trace1     = new int [traceMaxSize];
455                trace2     = new int [traceMaxSize];
456
457                int[] traceIndex     = new int [traceMaxSize];
458                int[] traceIterLevel = new int [traceMaxSize];
459
460                int ise11;
461                int ise12;
462                int ise21;
463                int ise22;
464
465                int ise1;
466                int ise2;
467
468                int gapMax=params.getMaxGapSize();
469
470                int iterDepth;
471                if ( gapMax > 0){
472                        iterDepth       =gapMax*2+1;
473                } else {
474                        iterDepth = traceMaxSize;
475                }
476                double[][] traceScore = new double[traceMaxSize][iterDepth];
477
478                nTraces =0;
479                long tracesLimit=(long)5e7;
480                double score =-1;
481                double score0 = -1;
482                double score1 = -1 ;
483                double score2 = -1;
484
485                int mse1;
486                int mse2;
487                int jgap;
488                int jdir;
489                int jse1=0;
490                int jse2=0;
491
492                int bestTracesMax=30;
493                bestTraces1 = new int[bestTracesMax][traceMaxSize];
494                bestTraces2 = new int[bestTracesMax][ traceMaxSize];
495                bestTracesN=new int [bestTracesMax];
496                bestTracesScores = new double [bestTracesMax];
497                for(int it=0; it<bestTracesMax; it++) {
498                        bestTracesN[it]=0;
499                        bestTracesScores[it]=100;
500                }
501
502                nBestTraces=0;
503                int newBestTrace=0;
504
505                double traceTotalScore=0;
506                double traceScoreMax =0;
507                double userRMSDMax = params.getMaxOptRMSD();
508                int kse1;
509                int kse2;
510
511                iterLoop:
512                        for(int iter=0; iter<nIter; iter++) {
513
514                                if(iter>2) {
515                                        if(oldBestTraceScore<=bestTraceScore) break;
516                                }
517                                oldBestTraceScore=bestTraceScore;
518
519                                if(iter==1) {
520                                        z0=zStrAlign(winSize, nBestTrace, bestTraceScore,
521                                                        bestTrace1[nBestTrace]+winSize-bestTrace1[0]+
522                                                        bestTrace2[nBestTrace]+winSize-bestTrace2[0]-
523                                                        nBestTrace*2*winSize);
524                                        if(z0<zThr) break;
525                                        nBestTrace0=nBestTrace;
526                                        nBestTrace=0;
527                                        bestTraceScore=100.0;
528
529                                        nTraces=0;
530                                }
531
532
533                                if(iter==0) {
534                                        ise11=0; ise12=nse1;
535                                        ise21=0; ise22=nse2;
536
537                                }
538                                else {
539                                        if(iter==1) {
540                                                ise11=bestTrace1[0]; ise12=bestTrace1[0]+1;
541                                                ise21=bestTrace2[0]; ise22=bestTrace2[0]+1;
542                                        }
543                                        else {
544                                                ise11=bestTrace1[0]-1; ise12=bestTrace1[0]+2;
545                                                ise21=bestTrace2[0]-1; ise22=bestTrace2[0]+2;
546                                        }
547                                        if(ise11<0) ise11=0;
548                                        if(ise12>nse1) ise12=nse1;
549                                        if(ise21<0) ise21=0;
550                                        if(ise22>nse2) ise22=nse2;
551                                }
552
553                                //System.out.println("ise1Loop: " + ise11 + " " + ise12 + " " + ise21 + " " + ise22);
554                                ise1Loop:
555                                        for(int ise1_=ise11; ise1_<ise12; ise1_++) {
556                                                ise2Loop:
557                                                        for(int ise2_=ise21; ise2_<ise22; ise2_++) {
558
559                                                                ise1=ise1_;
560                                                                ise2=ise2_;
561                                                                if(iter>1 && ise1==ise11+1 && ise2==ise21+1) continue ise2Loop;
562
563                                                                //if(ise2==ise21)       System.out.println(String.format("(%d, %d)",ise1, nTraces));
564
565
566                                                                if(iter==0 && (ise1>nse1-winSize*(nBestTrace-1) ||
567                                                                                ise2>nse2-winSize*(nBestTrace-1))) continue ise2Loop;
568
569                                                                if(mat[ise1][ise2]<0.0) continue ise2Loop;
570                                                                if(mat[ise1][ise2]>rmsdThr) continue ise2Loop;
571                                                                if (mat[ise1][ise2]>userRMSDMax) continue ise2Loop;
572                                                                nTrace=0;
573                                                                trace1[nTrace]=ise1;
574                                                                trace2[nTrace]=ise2;
575                                                                traceIndex[nTrace]=0;
576                                                                traceIterLevel[nTrace]=0;
577
578                                                                score0=mat[ise1][ise2];
579
580
581                                                                nTrace++;
582                                                                boolean isTraceUp=true;
583                                                                int traceIndex_=0;
584
585                                                                traceLoop:
586                                                                        while(nTrace>0) {
587
588                                                                                kse1=trace1[nTrace-1]+winSize;
589                                                                                kse2=trace2[nTrace-1]+winSize;
590
591                                                                                //System.out.println("isTraceUp " + isTraceUp + " " + nTrace + " " + kse1 + " " + kse2);
592
593                                                                                while(true) {
594                                                                                        if(kse1>nse1-winSize-1) break;
595                                                                                        if(kse2>nse2-winSize-1) break;
596                                                                                        if(mat[kse1][kse2]>=0.0) break;
597                                                                                        kse1++;
598                                                                                        kse2++;
599                                                                                }
600
601
602                                                                                traceIndex_=-1;
603
604                                                                                if(isTraceUp) {
605
606                                                                                        int nBestExtTrace=nTrace;
607                                                                                        double bestExtScore=100.0;
608
609
610                                                                                        // extension of the alignment path
611                                                                                        // condition 4, 5
612                                                                                        itLoop:
613                                                                                                for(int it=0; it<iterDepth; it++) {
614
615                                                                                                        jgap=(it+1)/2;
616                                                                                                        jdir=(it+1)%2;
617
618                                                                                                        if(jdir==0) {
619                                                                                                                mse1=kse1+jgap;
620                                                                                                                mse2=kse2;
621                                                                                                        }
622                                                                                                        else {
623                                                                                                                mse1=kse1;
624                                                                                                                mse2=kse2+jgap;
625                                                                                                        }
626
627                                                                                                        if(mse1>nse1-winSize-1) continue itLoop;
628                                                                                                        if(mse2>nse2-winSize-1) continue itLoop;
629
630                                                                                                        if(mat[mse1][mse2]<0.0)     continue itLoop;
631                                                                                                        if(mat[mse1][mse2]>rmsdThr) continue itLoop;
632                                                                                                        if(mat[mse1][mse2]>userRMSDMax) continue itLoop;
633
634                                                                                                        nTraces++;
635                                                                                                        if(nTraces>tracesLimit) {
636
637                                                                                                                return;
638                                                                                                        }
639
640                                                                                                        score=0.0;
641
642                                                                                                        //                                                                                                      if(!distAll) {
643                                                                                                                //System.out.println("getting score " + mse1 + " " + mse2 + " " + winSize + " " + jgap + " " + jdir + " " + it + " " + kse1 + " " + kse2);
644                                                                                                                score = getScoreFromDistanceMatrices(mse1,mse2,winSize);
645                                                                                                                //System.out.println("got score: " + score);
646                                                                                                                score1=score/(nTrace*winSize);
647
648                                                                                                                //                                                                                                      } else {
649                                                                                                                        //                                                                                                              // all dist
650                                                                                                                        //                                                                                                              for(int itrace=0; itrace<nTrace; itrace++) {
651                                                                                                                //                                                                                                                      for(int is1=0; is1<winSize; is1++)
652                                                                                                                //                                                                                                                              for(int is2=0; is2<winSize; is2++)
653                                                                                                                //                                                                                                                                      score+=Math.abs(dist1[trace1[itrace]+is1][mse1+is2]-
654                                                                                                                //                                                                                                                                                      dist2[trace2[itrace]+is1][mse2+is2]);
655                                                                                                                //                                                                                                              }
656                                                                                                                //                                                                                                              score1=score/(nTrace*winSize*winSize);
657                                                                                                                //                                                                                                      }
658
659
660                                                                                                                //System.out.println("up: " + nTrace + " "  + score + " " + score0 + " " + score1 + " " + winSize + " " + traceIndex_ + " " + it + " ");
661                                                                                                                if(score1>rmsdThrJoin)
662                                                                                                                        continue itLoop;
663                                                                                                                if(score1>userRMSDMax)
664                                                                                                                   continue itLoop;
665                                                                                                                score2=score1;
666
667                                                                                                                // this just got checked, no need to check again..
668                                                                                                                //if(score2>rmsdThrJoin)
669                                                                                                                //      continue itLoop;
670
671                                                                                                                if(nTrace>nBestExtTrace || (nTrace==nBestExtTrace &&
672                                                                                                                                score2<bestExtScore)) {
673                                                                                                                        //System.out.println("setting traceindex to " + it + " " + score2);
674                                                                                                                        bestExtScore=score2;
675                                                                                                                        nBestExtTrace=nTrace;
676                                                                                                                        traceIndex_=it;
677                                                                                                                        traceScore[nTrace-1][traceIndex_]=score1;
678                                                                                                                }
679
680                                                                                                }
681                                                                                }
682
683                                                                                if(traceIndex_!=-1) {
684                                                                                        jgap=(traceIndex_+1)/2;
685                                                                                        jdir=(traceIndex_+1)%2;
686                                                                                        if(jdir==0) {
687                                                                                                jse1=kse1+jgap;
688                                                                                                jse2=kse2;
689                                                                                        }
690                                                                                        else {
691                                                                                                jse1=kse1;
692                                                                                                jse2=kse2+jgap;
693                                                                                        }
694
695                                                                                        if(iter==0){
696
697                                                                                                score1=(traceScore[nTrace-1][traceIndex_]*winSizeComb2*nTrace+
698                                                                                                                mat[jse1][jse2]*winSizeComb1)/(winSizeComb2*nTrace+
699                                                                                                                                winSizeComb1);
700
701                                                                                                score2 = getScore2(jse1, jse2, traceScore, traceIndex_, traceIndex, winSizeComb1, winSizeComb2, score0, score1);
702
703                                                                                                if(score2>rmsdThrJoin)
704                                                                                                        traceIndex_=-1;
705                                                                                                else if ( score2 > userRMSDMax)
706                                                                                                   traceIndex_=-1;
707                                                                                                else {
708                                                                                                        traceScore[nTrace-1][traceIndex_]=score2;
709
710                                                                                                        traceTotalScore=score2;
711                                                                                                }
712
713                                                                                        }
714                                                                                        else {
715                                                                                                if(traceScoreMax>rmsdThrJoin && nBestTrace>=nBestTrace0)
716                                                                                                        traceIndex_=-1;
717                                                                                                traceTotalScore=traceScoreMax;
718                                                                                        }
719                                                                                }
720
721                                                                                //System.out.println("middle: " + nTrace + " " + score + " " + score0 + " " + score1 + "  " + score2  + " " + traceIndex_);
722
723                                                                                if(traceIndex_==-1) {
724                                                                                        //System.out.println("continue traceLoop " + nTrace);
725                                                                                        //if(iterLevel==1) break;
726                                                                                        nTrace--;
727                                                                                        isTraceUp=false;
728                                                                                        continue traceLoop;
729                                                                                }
730                                                                                else {
731                                                                                        traceIterLevel[nTrace-1]++;
732                                                                                        trace1[nTrace]=jse1;
733                                                                                        trace2[nTrace]=jse2;
734                                                                                        traceIndex[nTrace]=traceIndex_;
735                                                                                        traceIterLevel[nTrace]=0;
736                                                                                        nTrace++;
737                                                                                        isTraceUp=true;
738
739                                                                                        if(nTrace>nBestTrace ||
740                                                                                                        (nTrace==nBestTrace  &&
741                                                                                                                        bestTraceScore>traceTotalScore)) {
742
743                                                                                                for(int itrace=0; itrace<nTrace; itrace++) {
744                                                                                                        bestTrace1[itrace]=trace1[itrace];
745                                                                                                        bestTrace2[itrace]=trace2[itrace];
746                                                                                                }
747                                                                                                bestTraceScore=traceTotalScore;
748                                                                                                nBestTrace=nTrace;
749                                                                                        }
750
751                                                                                        if(iter==0) {
752                                                                                                //System.out.println("doing iter0 " + newBestTrace + " " + traceTotalScore + " " + bestTracesMax);
753                                                                                                newBestTrace = doIter0(newBestTrace,traceTotalScore, bestTracesMax);
754
755
756                                                                                        }
757                                                                                }
758                                                                        }
759                                                        }
760                                        }
761
762                                if ( isPrint) {
763                                        System.out.println("fragment length: " + params.getWinSize());
764                                        System.out.println("ntraces : " + nTraces );
765                                }
766
767
768
769                        }
770
771//              try {
772//                      Matrix m2 = new Matrix(traceScore).copy();
773//                      JPanel panel = GuiWrapper.getScaleableMatrixPanel(m2);
774//                      JFrame frame = new JFrame();
775//                      frame.addWindowListener(new WindowAdapter(){
776//                              public void windowClosing(WindowEvent e){
777//                                      JFrame f = (JFrame) e.getSource();
778//                                      f.setVisible(false);
779//                                      f.dispose();
780//                              }
781//                      });
782//
783//
784//                      frame.getContentPane().add(panel);
785//
786//                      frame.pack();
787//                      frame.setVisible(true);
788//              } catch (Exception e) {
789//                      e.printStackTrace();
790//              }
791
792
793                if ( params.isShowAFPRanges()){
794                        System.out.println("fragment length: " + params.getWinSize());
795                        System.out.println("ntraces : " + nTraces );
796
797                }
798
799        }
800
801        protected double getScore2(int jse1, int jse2, double[][] traceScore, int traceIndex_,int[] traceIndex,int winSizeComb1, int winSizeComb2, double score0, double score1 ) {
802
803
804
805                /*double score2=
806                        ((nTrace>1?traceScore[nTrace-2][traceIndex[nTrace-1]]:score0)
807                 *a[nTrace-1]+score1*(a[nTrace]-a[nTrace-1]))/a[nTrace];
808                 */
809                double val = 0;
810                if ( nTrace>1)
811                        val =traceScore[nTrace-2][traceIndex[nTrace-1]];
812                else
813                        val = score0;
814
815                double score2 =  (val * traceIndexContainer[nTrace-1]+score1*(traceIndexContainer[nTrace]-traceIndexContainer[nTrace-1]))/traceIndexContainer[nTrace];
816
817                //System.out.println("check: score0 " + score0 + " score 1 " + score1 + " sc2: " + score2 + " val: " + val + " nTrace:" + nTrace+ " " +  traceIndexContainer[nTrace-1] + " " +  traceIndexContainer[nTrace-1] + " " + traceIndexContainer[nTrace] );
818
819                return score2;
820
821
822        }
823
824        protected int doIter0(int newBestTrace, double traceTotalScore, double bestTracesMax) {
825
826
827                // only do the whole method if this criteria is fulfilled...
828                if(nTrace>bestTracesN[newBestTrace] ||
829                                (nTrace==bestTracesN[newBestTrace] &&
830                                                bestTracesScores[newBestTrace]>traceTotalScore)) {
831
832
833                        for(int itrace=0; itrace<nTrace; itrace++) {
834                                bestTraces1[newBestTrace][itrace]=trace1[itrace];
835                                bestTraces2[newBestTrace][itrace]=trace2[itrace];
836                                bestTracesN[newBestTrace]=nTrace;
837                                bestTracesScores[newBestTrace]=traceTotalScore;
838                                //System.out.println("bestTracesScrores ["+newBestTrace+"]=" +traceTotalScore);
839                        }
840
841                        if(nTrace>nBestTrace) nBestTrace=nTrace;
842
843                        if(nBestTraces<bestTracesMax) {
844                                nBestTraces++;
845                                newBestTrace++;
846                        }
847
848                        if(nBestTraces==bestTracesMax) {
849                                //System.out.println("nBestTraces == bestTracesmax " + nBestTraces);
850                                newBestTrace=0;
851                                double scoreTmp=bestTracesScores[0];
852                                int nTraceTmp=bestTracesN[0];
853                                for(int ir=1; ir<nBestTraces; ir++) {
854                                        if(bestTracesN[ir]<nTraceTmp ||
855                                                        (bestTracesN[ir]==nTraceTmp &&
856                                                                        scoreTmp<bestTracesScores[ir])) {
857                                                nTraceTmp=bestTracesN[ir];
858                                                scoreTmp=bestTracesScores[ir];
859                                                newBestTrace=ir;
860                                                //System.out.println("setting new bestTracesScore to " + ir + " " + scoreTmp);
861                                        }
862                                }
863                        }
864                }
865
866                //System.out.println("iter0 : " + newBestTrace + " " + bestTracesN[newBestTrace] + " " + traceTotalScore + " " + nTrace);
867
868
869
870                /*
871z=zStrAlign(winSize, nTrace, traceTotalScore,
872trace1[nTrace-1]-trace1[0]+trace2[nTrace-1]-trace2[0]-
8732*(nTrace-1)*winSize);
874if(z>bestTraceZScore) {
875for(int itrace=0; itrace<nTrace; itrace++) {
876bestTrace1[itrace]=trace1[itrace];
877bestTrace2[itrace]=trace2[itrace];
878}
879bestTraceZScore=z;
880bestTraceScore=*(traceScore[nTrace-2]+traceIndex_);
881nBestTrace=nTrace;
882}
883                 */
884                return newBestTrace;
885
886        }
887
888
889        protected double getScoreFromDistanceMatrices(int mse1, int mse2,int winSize) {
890
891                double score = 0;
892                // (winSize) "best" dist
893
894                // reduce sign. values to C code.. 6 digits..
895
896                for(int itrace=0; itrace<nTrace; itrace++) {
897                        score+=  Math.abs(dist1[trace1[itrace]][mse1]-
898                                        dist2[trace2[itrace]][mse2]);
899
900                        score+=  Math.abs(dist1[trace1[itrace]+winSize-1]
901                                                [mse1+winSize-1]-
902                                                dist2[trace2[itrace]+winSize-1][mse2+winSize-1]);
903
904                        for(int id=1; id<winSize-1; id++)
905                                score+=  Math.abs(dist1[trace1[itrace]+id][mse1+winSize-1-id]-
906                                                dist2[trace2[itrace]+id][mse2+winSize-1-id]);
907
908                }
909
910                return score;
911        }
912
913        public void nextStep( AFPChain afpChain,
914                        Atom[] ca1, Atom[] ca2) throws StructureException{
915
916
917                if(nBestTrace>0) {
918                        checkBestTraces(afpChain,ca1,ca2);
919                } else {
920                        noBestTrace();
921                }
922
923                convertAfpChain(afpChain, ca1, ca2);
924                AFPAlignmentDisplay.getAlign(afpChain, ca1, ca2);
925                double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2,false);
926                afpChain.setTMScore(tmScore);
927        }
928
929
930
931        // this part is modified from the original CeCalculator
932        @SuppressWarnings("unused")
933        private void checkBestTraces( AFPChain afpChain,
934                        Atom[] ca1, Atom[] ca2) throws StructureException{
935
936                z=0.0;
937
938
939                int nGaps;
940                int winSize = params.getWinSize();
941                int nse1 = ca1.length;
942                int nse2 = ca2.length;
943                int traceMaxSize=nse1<nse2?nse1:nse2;
944                int idir;
945
946
947                align_se1=new int [nse1+nse2];
948                align_se2=new int [nse1+nse2];
949                lcmp = 0;
950
951                // we now support alignment using any particular atoms..
952
953                Atom[] strBuf1 = new Atom[traceMaxSize];
954                Atom[] strBuf2 = new Atom[traceMaxSize];
955
956                double rmsdNew;
957
958
959
960                // removing some loops that are run in orig CE
961                // and which did not do anything
962                if ( debug ){
963                        checkPrintRmsdNew(traceMaxSize, winSize, ca1, ca2);
964                }
965
966                double rmsd=100.0;
967
968                int iBestTrace=0;
969
970                for(int ir=0; ir<nBestTraces; ir++) {
971                        if(bestTracesN[ir]!=nBestTrace) continue;
972
973                        rmsdNew = getRMSDForBestTrace(ir, strBuf1, strBuf2, bestTracesN,bestTraces1, bestTrace2,winSize,ca1,ca2);
974                        if ( isPrint)
975                                System.out.println(String.format("%d %d %d %.2f", ir, bestTracesN[ir], nBestTrace, rmsdNew));
976
977                        if(rmsd>rmsdNew) {
978                                iBestTrace=ir;
979                                rmsd=rmsdNew;
980                                //System.out.println(" iBestTrace:" + iBestTrace + " new rmsd = " + rmsd);
981                        }
982                }
983                for(int it=0; it<bestTracesN[iBestTrace]; it++) {
984                        bestTrace1[it]=bestTraces1[iBestTrace][it];
985                        bestTrace2[it]=bestTraces2[iBestTrace][it];
986                }
987
988                //System.out.println("iBestTrace: "+iBestTrace+" = bestTracesScores " + bestTracesScores[iBestTrace]);
989
990                nBestTrace=bestTracesN[iBestTrace];
991
992                bestTraceScore=bestTracesScores[iBestTrace];
993
994
995                //printf("\nOptimizing gaps...\n");
996
997                int[] traceLen=new int [traceMaxSize];
998                bestTraceLen=new int [traceMaxSize];
999
1000
1001                int strLen=0;
1002
1003                int jt;
1004                strLen=0;
1005                nGaps=0;
1006                nTrace=nBestTrace;
1007
1008                for(jt=0; jt<nBestTrace; jt++) {
1009                        trace1[jt]=bestTrace1[jt];
1010                        trace2[jt]=bestTrace2[jt];
1011                        traceLen[jt]=winSize;
1012
1013                        if(jt<nBestTrace-1) {
1014                                nGaps+=bestTrace1[jt+1]-bestTrace1[jt]-winSize+
1015                                bestTrace2[jt+1]-bestTrace2[jt]-winSize;
1016                        }
1017                }
1018                nBestTrace=0;
1019                for(int it=0; it<nTrace; ) {
1020                        int cSize=traceLen[it];
1021                        for(jt=it+1; jt<nTrace; jt++) {
1022                                if(trace1[jt]-trace1[jt-1]-traceLen[jt-1]!=0 ||
1023                                                trace2[jt]-trace2[jt-1]-traceLen[jt-1]!=0) break;
1024                                cSize+=traceLen[jt];
1025                        }
1026                        bestTrace1[nBestTrace]=trace1[it];
1027                        bestTrace2[nBestTrace]=trace2[it];
1028                        bestTraceLen[nBestTrace]=cSize;
1029                        nBestTrace++;
1030                        strLen+=cSize;
1031                        it=jt;
1032                }
1033
1034
1035                int is=0;
1036                for(jt=0; jt<nBestTrace; jt++) {
1037                        for(int i=0; i<bestTraceLen[jt]; i++) {
1038                                setStrBuf(strBuf1,is+i,ca1,bestTrace1[jt]+i );
1039                                setStrBuf(strBuf2,is+i,ca2,bestTrace2[jt]+i);
1040                        }
1041                        is+=bestTraceLen[jt];
1042                }
1043                //sup_str(strBuf1, strBuf2, strLen, d_);
1044
1045                rmsd=calc_rmsd(strBuf1, strBuf2, strLen,true);
1046
1047                if ( isPrint)
1048                        System.out.println("got first rmsd: " + rmsd);
1049                boolean isCopied=false;
1050
1051                outer_loop:
1052                        for(int it=1; it<nBestTrace; it++) {
1053
1054                                /* not needed...
1055                        int igap;
1056                        if(bestTrace1[it]-bestTrace1[it-1]-bestTraceLen[it-1]>0) igap=0;
1057                        if(bestTrace2[it]-bestTrace2[it-1]-bestTraceLen[it-1]>0) igap=1;
1058                                 */
1059
1060
1061                                boolean wasBest=false;
1062                                main_loop:
1063                                        for(idir=-1; idir<=1; idir+=2) {
1064                                                if(wasBest) break;
1065
1066                                                inner_loop:
1067                                                        for(int idep=1; idep<=winSize/2; idep++) {
1068
1069                                                                // isCopied indicates that bestTrace has changed and needs to be re-copied
1070                                                                if(!isCopied)
1071                                                                        for(jt=0; jt<nBestTrace; jt++) {
1072                                                                                trace1[jt]=bestTrace1[jt];
1073                                                                                trace2[jt]=bestTrace2[jt];
1074                                                                                traceLen[jt]=bestTraceLen[jt];
1075                                                                        }
1076                                                                isCopied=false;
1077
1078                                                                // Move an atom from the previous trace to the current on, or vice versa
1079                                                                traceLen[it-1]+=idir;
1080                                                                traceLen[it]-=idir;
1081                                                                trace1[it]+=idir;
1082                                                                trace2[it]+=idir;
1083
1084                                                                // Copy atoms from the current trace into strBuf
1085                                                                is=0;
1086                                                                for(jt=0; jt<nBestTrace; jt++) {
1087                                                                        for(int i=0; i<traceLen[jt]; i++) {
1088                                                                                if(ca1[trace1[jt]+i].getX()>1e10 || ca2[trace2[jt]+i].getX()>1e10)
1089                                                                                        continue main_loop;
1090                                                                                strBuf1[is+i]=ca1[trace1[jt]+i];
1091                                                                                strBuf2[is+i]=ca2[trace2[jt]+i];
1092                                                                        }
1093                                                                        is+=traceLen[jt];
1094                                                                }
1095                                                                // Check new RMSD
1096                                                                //sup_str(strBuf1, strBuf2, strLen, d_);
1097                                                                rmsdNew=calc_rmsd(strBuf1, strBuf2, strLen, true);
1098                                                                //System.out.println(String.format("step %d %d %d %.2f old: %.2f", it, idir, idep, rmsdNew, rmsd));
1099
1100                                                                // Update best trace if RMSD improved
1101                                                                if(rmsdNew<rmsd) {
1102
1103                                                                        for(jt=0; jt<nBestTrace; jt++) {
1104                                                                                bestTrace1[jt]  = trace1[jt];
1105                                                                                bestTrace2[jt]  = trace2[jt];
1106                                                                                bestTraceLen[jt]= traceLen[jt];
1107                                                                        }
1108                                                                        isCopied=true;
1109                                                                        wasBest=true;
1110                                                                        rmsd=rmsdNew;
1111                                                                        continue inner_loop;
1112                                                                }
1113                                                                // AP
1114                                                                //bad_ca: break;
1115                                                                continue main_loop;
1116                                                        }
1117                                        }
1118                        }
1119                rmsdNew=calc_rmsd(strBuf1, strBuf2, strLen,true);
1120                if ( isPrint)
1121                        System.out.println("rmsdNew: " + rmsdNew + " rmsd " + rmsd);
1122                afpChain.setTotalRmsdIni(rmsdNew);
1123                afpChain.setTotalLenIni(strBuf1.length);
1124
1125
1126                nAtom = strLen;
1127
1128                //System.out.println("zStrAlign: " + winSize + " strLen " + strLen  + " s/w " + (strLen/winSize) + " " + bestTraceScore + " " + nGaps);
1129                z=zStrAlign(winSize, strLen/winSize, bestTraceScore, nGaps);
1130
1131                if(params.isShowAFPRanges()) {
1132                        System.out.println("win size: " + winSize + " strLen/winSize: " + strLen/winSize + " best trace score: " + String.format("%.2f",bestTraceScore) + " nr gaps: " + nGaps + " nr residues: " + nAtom);
1133
1134                        System.out.println(String.format("size=%d rmsd=%.2f z=%.1f gaps=%d(%.1f%%) comb=%d",
1135                                        nAtom, rmsd, z, nGaps, nGaps*100.0/nAtom,
1136                                        nTraces));
1137
1138                        System.out.println("Best Trace, before optimization");
1139                        for(int k=0; k<nBestTrace; k++)
1140                                System.out.println(String.format("(%d,%d,%d) ", bestTrace1[k]+1, bestTrace2[k]+1,
1141                                                bestTraceLen[k]));
1142
1143                }
1144
1145                // start to convert CE internal datastructure to generic AFPChain one...
1146                List<AFP> afpSet = new ArrayList<AFP>();
1147                for (int afp=0;afp<nBestTrace;afp++){
1148                        // fill in data from nBestTrace into AFP
1149
1150                        AFP afpI = new AFP();
1151
1152                        afpI.setFragLen(bestTraceLen[afp]);
1153                        afpI.setP1(bestTrace1[afp]+1);
1154                        afpI.setP2(bestTrace2[afp]+1);
1155
1156                        afpSet.add( afpI);
1157                }
1158
1159                afpChain.setAfpSet(afpSet);
1160
1161                //System.out.println("z:"+z + " zThr" + zThr+ " bestTraceScore " + bestTraceScore + " " + nGaps );
1162                if(params.isOptimizeAlignment() && z>=zThr) {
1163                        nGaps = optimizeSuperposition(afpChain,nse1, nse2, strLen, rmsd, ca1, ca2,nGaps,strBuf1,strBuf2);
1164                        //            if(isPrint) {
1165                        //              /*
1166                        //              FILE *f=fopen("homologies", "a");
1167                        //              fprintf(f, "%s(%d) %s(%d) %3d %4.1f %4.1f %d(%d) ",
1168                        //                      name1, nse1, name2, nse2, nAtom, rmsd, z,
1169                        //                      nGaps, nGaps*100/nAtom);
1170                        //              for(int k=0; k<nBestTrace; k++)
1171                        //                fprintf(f, "(%d,%d,%d) ", bestTrace1[k]+1, bestTrace2[k]+1,
1172                        //                        bestTraceLen[k]);
1173                        //              fprintf(f, "\n");
1174                        //              fclose(f);
1175                        //              */
1176                        //            }
1177                }
1178                else {
1179                        int lali_x_ = 0;
1180                        for(int k=0; k<nBestTrace; k++) {
1181                                for(int l=0; l<bestTraceLen[k]; l++) {
1182                                        align_se1[lcmp+l]=bestTrace1[k]+l;
1183                                        align_se2[lcmp+l]=bestTrace2[k]+l;
1184                                }
1185                                lali_x_+=bestTraceLen[k];
1186                                lcmp+=bestTraceLen[k];
1187                                if(k<nBestTrace-1) {
1188                                        if(bestTrace1[k]+bestTraceLen[k]!=bestTrace1[k+1])
1189                                                for(int l=bestTrace1[k]+bestTraceLen[k]; l<bestTrace1[k+1]; l++) {
1190                                                        align_se1[lcmp]=l;
1191                                                        align_se2[lcmp]=-1;
1192                                                        lcmp++;
1193                                                }
1194                                        if(bestTrace2[k]+bestTraceLen[k]!=bestTrace2[k+1])
1195                                                for(int l=bestTrace2[k]+bestTraceLen[k]; l<bestTrace2[k+1]; l++) {
1196                                                        align_se1[lcmp]=-1;
1197                                                        align_se2[lcmp]=l;
1198                                                        lcmp++;
1199                                                }
1200                                }
1201                        }
1202                        nAtom=lali_x_;
1203                        afpChain.setTotalRmsdOpt(afpChain.getTotalRmsdIni());
1204                }
1205
1206                timeEnd = System.currentTimeMillis();
1207                long time_q=(timeEnd-timeStart);
1208
1209                double gapsP = ( nGaps*100.0/nAtom) ;
1210                if(isPrint) {
1211                        String msg = String.format("Alignment length = %d Rmsd = %.2fA Z-Score = %.1f Gaps = %d(%.1f%%)",nAtom,rmsd,z,nGaps, gapsP);
1212                        System.out.println(msg + " CPU = " + time_q);
1213                }
1214
1215                //      if ( params.isShowAFPRanges()){
1216
1217                // this is actually the final alignment...
1218                //         System.out.println("Best Trace: (index1,index2,len)");
1219                //         for(int k=0; k<nBestTrace; k++)
1220                //            System.out.println(
1221                //                  String.format("(%d,%d,%d) ", bestTrace1[k]+1, bestTrace2[k]+1, bestTraceLen[k]));
1222                //
1223                //
1224                //
1225                //      }
1226
1227                afpChain.setCalculationTime(time_q);
1228                afpChain.setGapLen(nGaps);
1229
1230                int[] optLen = new int[]{nAtom};
1231                afpChain.setOptLen(optLen);
1232                afpChain.setOptLength(nAtom);
1233                afpChain.setAlnLength(lcmp);
1234
1235                afpChain.setProbability(z);
1236
1237
1238        }
1239
1240        /** set the Atoms for a particular residue position.
1241         * Requires that atom.getParent returns the correct group!
1242         * take care during cloning of atoms. Best to use StructureTools.cloneCaAtoms();
1243         *
1244         * @param strBuf
1245         * @param i
1246         * @param ca
1247         * @param j
1248         */
1249        private void setStrBuf(Atom[] strBuf, int i, Atom[] ca, int j) {
1250                // TODO Auto-generated method stub
1251                //TODO
1252                Group parent = ca[j].getGroup();
1253                int pos = 0;
1254                String atomName = ca[j].getName();
1255
1256                Atom a = null;
1257
1258                        a= parent.getAtom(atomName);
1259                        if ( a != null){
1260                                strBuf[i]=a;
1261                        } else {
1262                                //      probably a GLY and no CB was found...
1263                                //e.printStackTrace();
1264                        }
1265                strBuf[i+pos] = a;
1266                pos++;
1267
1268
1269
1270        }
1271
1272        // TODO:  consider all requested Atoms?
1273        private double getRMSDForBestTrace(int ir, Atom[] strBuf1, Atom[] strBuf2,
1274                        int[] bestTracesN2, int[][] bestTraces12, int[] bestTrace22,
1275                        int winSize,Atom[] ca1, Atom[] ca2 ) throws StructureException {
1276                int is=0;
1277                for(int jt=0; jt<bestTracesN[ir]; jt++) {
1278                        for(int i=0; i<winSize; i++) {
1279
1280                                setStrBuf(strBuf1, is+i, ca1, bestTraces1[ir][jt]+i);
1281                                setStrBuf(strBuf2, is+i, ca2, bestTraces2[ir][jt]+i);
1282                        }
1283                        is+=winSize;
1284                }
1285                //sup_str(strBuf1, strBuf2, bestTracesN[ir]*winSize, d_);
1286                double rmsdNew=calc_rmsd(strBuf1, strBuf2, bestTracesN[ir]*winSize, true);
1287                return rmsdNew;
1288
1289        }
1290
1291
1292
1293
1294
1295
1296        /** calc initial RMSD for bestTrace1 in debug only
1297         *
1298         */
1299        private void checkPrintRmsdNew(int traceMaxSize, int winSize, Atom[] ca1, Atom[] ca2) throws StructureException{
1300
1301                int is = 0;
1302                // calc initial RMSD for bestTrace1
1303                Atom[] strBuf1 = new Atom[traceMaxSize];
1304                Atom[] strBuf2 = new Atom[traceMaxSize];
1305                for(int jt=0; jt<nBestTrace; jt++) {
1306                        for(int i=0; i<winSize; i++) {
1307                                setStrBuf(strBuf1, is+i, ca1, bestTrace1[jt]+i );
1308                                setStrBuf(strBuf2, is+i, ca2, bestTrace2[jt]+i );
1309                        }
1310                        is+=winSize;
1311                }
1312
1313                //sup_str(strBuf1, strBuf2, nBestTrace*winSize, d_);
1314                double rmsdNew=calc_rmsd(strBuf1, strBuf2, nBestTrace*winSize, true);
1315                //afpChain.setTotalRmsdIni(rmsdNew);
1316
1317
1318                if ( isPrint){
1319                        System.out.println("rmsdNew after trace: " +rmsdNew);
1320
1321                        for(int k=0; k<nBestTrace; k++)
1322                                System.out.println(String.format("(%d,%d,%d) ", bestTrace1[k]+1, bestTrace2[k]+1,8));
1323                }
1324
1325                if ( isPrint){
1326                        System.out.println("best traces: " + nBestTraces);
1327                }
1328
1329
1330        }
1331
1332
1333
1334
1335
1336
1337        private static char getOneLetter(Group g){
1338
1339                if (g==null) return StructureTools.UNKNOWN_GROUP_LABEL;
1340
1341                return StructureTools.get1LetterCode(g.getPDBName());
1342        }
1343
1344
1345
1346        private int optimizeSuperposition(AFPChain afpChain, int nse1, int nse2, int strLen, double rmsd, Atom[] ca1, Atom[] ca2,int nGaps,
1347                        Atom[] strBuf1, Atom[] strBuf2 ) throws StructureException {
1348
1349                //System.out.println("optimizing Superimposition...");
1350
1351                //nAtom=strLen;
1352                // optimization on superposition
1353                Atom[] ca3=new Atom[nse2];
1354
1355
1356                double rmsdLen  = 0.0;
1357
1358                // this flag tests if the RMSDLen has been assigned.
1359                // this is to enforce that the alignment ends up
1360                // smaller than 95% of the original alignment.
1361                // +/-
1362                boolean isRmsdLenAssigned=false;
1363                int nAtomPrev=-1;
1364
1365                double oRmsdThr = params.getORmsdThr();
1366
1367                double distanceIncrement = params.getDistanceIncrement();
1368                double maxUserRMSD = params.getMaxOptRMSD();
1369                nAtom=0;
1370                int counter = -1;
1371
1372                int maxNrIterations = params.getMaxNrIterationsForOptimization();
1373                while((nAtom<strLen*0.95 ||
1374                                (isRmsdLenAssigned && rmsd<rmsdLen*1.1 && nAtomPrev!=nAtom)) && ( counter< maxNrIterations)) {
1375
1376                        counter++;
1377                        if ( debug)
1378                           System.out.println("nAtom: " + nAtom + " " + nAtomPrev + " " + rmsdLen + " " + isRmsdLenAssigned + " strLen:" + strLen);
1379                        nAtomPrev=nAtom;
1380                        oRmsdThr += distanceIncrement;
1381
1382                        rot_mol(ca2, ca3, nse2, r,t);
1383
1384                        for(int ise1=0; ise1<nse1; ise1++) {
1385                                for(int ise2=0; ise2<nse2; ise2++) {
1386
1387                                        //mat[ise1][ise2]=-0.001;
1388
1389                                        // this needs to be a parameter...
1390
1391
1392                                        double dist = getDistanceWithSidechain(ca1[ise1], ca3[ise2]);
1393                                        mat[ise1][ise2] = oRmsdThr - dist;
1394
1395                                        //double distold = Calc.getDistance(ca1[ise1],ca3[ise2]);
1396                                        //double scoreOld  = oRmsdThr - distold ;
1397                                        //mat[ise1][ise2] = scoreOld;
1398                                        //mat[ise1][ise2] = oRmsdThr - Calc.getDistance(ca1[ise1],ca3[ise2]);
1399
1400                                        //if ( counter == 0 &&  ise1 == ise2) {
1401
1402                                        // System.out.println("mat[" + ise1 + "][" + ise2 + "] " + mat[ise1][ise2] + " scoreOld:" + scoreOld + " oRmsdThr: " + oRmsdThr +" dist: " + dist + " distold:" + distold );
1403                                        // }
1404
1405
1406                                }
1407                        }
1408
1409                        if ( params.getScoringStrategy() == CeParameters.ScoringStrategy.SEQUENCE_CONSERVATION){
1410                                mat = updateMatrixWithSequenceConservation(mat,ca1,ca2, params);
1411                        }
1412
1413                        mat = notifyMatrixListener(mat);
1414
1415                        double gapOpen = params.getGapOpen();
1416                        double gapExtension = params.getGapExtension();
1417
1418                        double score = dpAlign( nse1, nse2, gapOpen , gapExtension , false, false);
1419
1420                        if (debug)
1421                                System.out.println("iter: "+ counter + "  score:"  + score + " " + " nAtomPrev: " + nAtomPrev + " nAtom:" + nAtom + " oRmsdThr: " + oRmsdThr);
1422
1423                        afpChain.setAlignScore(score);
1424
1425
1426                        nAtom=0; nGaps=0;
1427                        for(int ia=0; ia<lcmp; ia++) {
1428                                if(align_se1[ia]!=-1 && align_se2[ia]!=-1) {
1429
1430                                        strBuf1[nAtom]=ca1[align_se1[ia]];
1431                                        strBuf2[nAtom]=ca2[align_se2[ia]];
1432
1433                                        nAtom++;
1434
1435                                }
1436                                else {
1437                                        nGaps++;
1438                                }
1439                        }
1440
1441
1442
1443                        if(nAtom<4) continue;
1444
1445                        //sup_str(strBuf1, strBuf2, nAtom, _d);
1446                        // here we don't store the rotation matrix for the user!
1447                        rmsd= calc_rmsd(strBuf1, strBuf2, nAtom,false);
1448                        if ( isPrint )
1449                                System.out.println("iter: " + counter + " nAtom " + nAtom + " rmsd: " + rmsd);
1450                        //afpChain.setTotalRmsdOpt(rmsd);
1451                        //System.out.println("rmsd: " + rmsd);
1452
1453                        if(!(nAtom<strLen*0.95) && (!isRmsdLenAssigned)) {
1454                                rmsdLen=rmsd;
1455                                isRmsdLenAssigned=true;
1456                        }
1457                        //System.out.println(String.format("nAtom %d %d rmsd %.1f", nAtom, nAtomPrev, rmsd));
1458
1459
1460                        afpChain.setBlockRmsd(new double[]{rmsd});
1461                        afpChain.setOptRmsd(new double[]{rmsd});
1462                        afpChain.setTotalRmsdOpt(rmsd);
1463                        afpChain.setChainRmsd(rmsd);
1464
1465                        if ( rmsd >= maxUserRMSD) {
1466                      break;
1467                        }
1468
1469                }
1470
1471
1472
1473                //System.out.println("done optimizing");
1474                /*
1475                nAtom=0; nGaps=0;
1476                for(int ia=0; ia<lcmp; ia++)
1477                if(align_se1[ia]!=-1 && align_se2[ia]!=-1) {
1478                if(ca1[align_se1[ia]].X<1e10 && ca2[align_se2[ia]].X<1e10) {
1479                strBuf1[nAtom]=ca1[align_se1[ia]];
1480                strBuf2[nAtom]=ca2[align_se2[ia]];
1481                nAtom++;
1482                }
1483                }
1484                else {
1485                nGaps++;
1486                }
1487
1488                sup_str(strBuf1, strBuf2, nAtom, _d);
1489                rmsd=calc_rmsd(strBuf1, strBuf2, nAtom, _d);
1490                 */
1491                nBestTrace=0;
1492                boolean newBestTrace=true;
1493                for(int ia=0; ia<lcmp; ia++) {
1494                        if(align_se1[ia]!=-1 && align_se2[ia]!=-1) {
1495                                //System.out.println(" " +align_se1[ia] + " " + align_se2[ia]);
1496
1497                                if(newBestTrace) {
1498                                        bestTrace1[nBestTrace]=align_se1[ia];
1499                                        bestTrace2[nBestTrace]=align_se2[ia];
1500                                        bestTraceLen[nBestTrace]=0;
1501                                        newBestTrace=false;
1502                                        nBestTrace++;
1503                                }
1504                                bestTraceLen[nBestTrace-1]++;
1505
1506                        }
1507                        else {
1508                                newBestTrace=true;
1509                        }
1510                }
1511
1512                return nGaps;
1513
1514                // end of optimization on superposition
1515
1516        }
1517
1518        /** Modifies an alignment matrix by favoring the alignment of similar and identical amino acids and penalizing the alignment of unrelated ones.
1519         *
1520         * @param max alignment matrix
1521         * @param ca1 Atoms for protein 1
1522         * @param ca2 Atoms for Protein 2
1523         * @param params alignment parameters
1524         * @return modified alignment matrix
1525         */
1526        public static double[][] updateMatrixWithSequenceConservation(double[][] max, Atom[] ca1, Atom[] ca2, CeParameters params) {
1527                Matrix origM = new Matrix(max);
1528
1529                SubstitutionMatrix<AminoAcidCompound> substMatrix =
1530                        params.getSubstitutionMatrix();
1531
1532                int internalScale = 1;
1533                if ( substMatrix instanceof ScaledSubstitutionMatrix) {
1534                        ScaledSubstitutionMatrix scaledMatrix = (ScaledSubstitutionMatrix) substMatrix;
1535                        internalScale = scaledMatrix.getScale();
1536                }
1537
1538
1539                AminoAcidCompoundSet set = AminoAcidCompoundSet.getAminoAcidCompoundSet();
1540
1541                for (int i = 0 ; i < origM.getRowDimension() ; i++){
1542                        for ( int j =0; j < origM.getColumnDimension() ; j ++ ) {
1543                                double val = origM.get(i,j);
1544                                Atom a1 = ca1[i];
1545                                Atom a2 = ca2[j];
1546
1547                                AminoAcidCompound ac1 =
1548                                        set.getCompoundForString(a1.getGroup().getChemComp().getOne_letter_code());
1549                                AminoAcidCompound ac2 =
1550                                        set.getCompoundForString(a2.getGroup().getChemComp().getOne_letter_code());
1551
1552
1553                                if ( ac1 == null || ac2 == null)
1554                                        continue;
1555
1556                                short aaScore = substMatrix.getValue(ac1,ac2);
1557
1558                                double weightedScore = (aaScore / internalScale) * params.getSeqWeight();
1559
1560
1561                                val += weightedScore;
1562                                origM.set(i,j,val);
1563
1564                        }
1565                }
1566                max = origM.getArray();
1567
1568                //SymmetryTools.showMatrix((Matrix)origM.clone(), "in optimizer "  + loopCount  );
1569                //SymmetryTools.showMatrix(origM, "iteration  matrix " + loopCount + " after");
1570
1571                return max;
1572        }
1573
1574        private double[][] notifyMatrixListener(double[][] mat2) {
1575                for (MatrixListener li : matrixListeners) {
1576                        mat2 = li.matrixInOptimizer(mat2);
1577                }
1578                return mat2;
1579        }
1580
1581        private boolean[][] notifyBreakFlagListener(boolean[][] brkFlag){
1582                for (MatrixListener li : matrixListeners) {
1583                        brkFlag = li.initializeBreakFlag(brkFlag);
1584                }
1585                return brkFlag;
1586        }
1587
1588        public void addMatrixListener(MatrixListener li){
1589                matrixListeners.add(li);
1590        }
1591
1592
1593        /**
1594         * On input, mat[i][j] should give the score for aligning positions i and j.
1595         * On output, mat[i][j] gives the maximum score possible for aligning 1..i
1596         * of protein 1 with 1..j of protein 2.
1597         *
1598         * @param nSeq1 The length of protein 1 (mat.length)
1599         * @param nSeq2 The length of protein 2 (mat[0].length)
1600         * @param gapI gap initiation penalty
1601         * @param gapE gap extension penalty
1602         * @param isGlobal1 The alignment is global for protein 1
1603         * @param isGlobal2 The alignment is global for protein 2
1604         * @return The maximum score
1605         */
1606        private double dpAlign(int nSeq1, int nSeq2, double gapI, double gapE,
1607                        boolean isGlobal1, boolean isGlobal2) {
1608
1609                // isGlobal1,isGlobal2 are always false...
1610
1611                int i, j, is, js, iMax, jMax, k;
1612                boolean ge=(gapE!=0.0?true:false);
1613                double sum, sum_ret, sum_brk;
1614
1615                boolean[][] brk_flg=new boolean [nSeq1][nSeq2];
1616                for(i=0; i<nSeq1; i++) brk_flg[i]=new boolean [nSeq2];
1617
1618                brk_flg = notifyBreakFlagListener(brk_flg);
1619
1620                // ge = true here...
1621                /*
1622                  for(i=0; i<nSeq1; i++)
1623                   {
1624                     printf("\n");
1625                     for(j=0; j<nSeq2; j++)
1626                       {
1627                         printf("%4d", (int)(*(mat[i]+j)*10));
1628                       }
1629                   }
1630                 printf("\n\n\n");
1631                 */
1632                if(!ge)
1633                {
1634                        for(i=nSeq1-1; i>=0; i--)
1635                                for(j=nSeq2-1; j>=0; j--)
1636                                {
1637                                        brk_flg[i][j]=false;
1638                                        if(j<nSeq2-1 && i<nSeq1-1)
1639                                        {
1640                                                sum=mat[i+1][j+1];
1641                                        }
1642                                        else
1643                                        {
1644                                                sum=0.0;
1645                                                if((isGlobal1 && i!=nSeq1-1) || (isGlobal2 && j!=nSeq2-1))
1646                                                        sum=-gapI;
1647                                        }
1648                                        if(j+1<nSeq2)
1649                                                for(k=i+2; k<nSeq1; k++)
1650                                                {
1651                                                        if(mat[k][j+1]-gapI>sum)
1652                                                                sum=mat[k][j+1]-gapI;
1653                                                }
1654                                        if(i+1<nSeq1)
1655                                                for(k=j+2; k<nSeq2; k++)
1656                                                {
1657                                                        if(mat[i+1][k]-gapI>sum)
1658                                                                sum=mat[i+1][k]-gapI;
1659                                                }
1660                                        sum+=mat[i][j];
1661                                        sum_brk=(isGlobal1?-gapI:0.0)+(isGlobal2?-gapI:0.0);
1662                                        if(sum<sum_brk)
1663                                        {
1664                                                sum=sum_brk;
1665                                                brk_flg[i][j]=true;
1666                                                //System.out.println("break at: " + i + " " + j);
1667                                        }
1668                                        mat[i][j]=sum;
1669                                }
1670                }
1671                else
1672                {
1673                        for(i=nSeq1-1; i>=0; i--)
1674                                for(j=nSeq2-1; j>=0; j--)
1675                                {
1676                                        brk_flg[i][j]=false;
1677                                        if(j<nSeq2-1 && i<nSeq1-1)
1678                                        {
1679                                                sum=mat[i+1][j+1];
1680                                        }
1681                                        else
1682                                        {
1683                                                sum=0.0;
1684                                                if(isGlobal1 && i!=nSeq1-1) sum=-gapI-gapE*(nSeq1-i-1);
1685                                                if(isGlobal2 && j!=nSeq2-1) sum=-gapI-gapE*(nSeq2-j-1);
1686                                        }
1687                                        if(j+1<nSeq2)
1688                                                for(k=i+2; k<nSeq1; k++)
1689                                                        if(mat[k][j+1]-gapI-gapE*(k-i-1)>sum)
1690                                                                sum=mat[k][j+1]-gapI-gapE*(k-i-1);
1691                                        if(i+1<nSeq1)
1692                                                for(k=j+2; k<nSeq2; k++)
1693                                                        if(mat[i+1][k]-gapI-gapE*(k-j-1)>sum)
1694                                                                sum=mat[i+1][k]-gapI-gapE*(k-j-1);
1695                                        sum+=mat[i][j];
1696                                        sum_brk=(isGlobal1?(-gapI-gapE*(nSeq1-1-i)):0.0)+(isGlobal2?(-gapI-gapE*(nSeq2-1-j)):0.0);
1697                                        if(sum<sum_brk)
1698                                        {
1699                                                sum=sum_brk;
1700                                                brk_flg[i][j]=true;
1701                                        }
1702                                        mat[i][j]=sum;
1703                                }
1704                }
1705
1706                //              if (debug ){
1707                //                      ScaleableMatrixPanel smp = new ScaleableMatrixPanel();
1708                //                      JFrame frame = new JFrame("CE alignment matrix in dpAlign " );
1709                //                      frame.addWindowListener(new WindowAdapter(){
1710                //                              public void windowClosing(WindowEvent e){
1711                //                                      JFrame f = (JFrame) e.getSource();
1712                //                                      f.setVisible(false);
1713                //                                      f.dispose();
1714                //                              }
1715                //
1716                //                      });
1717                //
1718                //                      smp.getMatrixPanel().setScalevalue(100);
1719                //                      Matrix mx = (Matrix) new Matrix(mat).clone();
1720                //                      smp.setMatrix(mx);
1721                //
1722                //                      frame.getContentPane().add(smp);
1723                //
1724                //                      frame.pack();
1725                //                      frame.setVisible(true);
1726                //              }
1727
1728                /*
1729                 for(i=0; i<nSeq1; i++)
1730                   {
1731                     printf("\n");
1732                     for(j=0; j<nSeq2; j++)
1733                       {
1734                         printf("%4d", (int)(*(mat[i]+j)*10));
1735                       }
1736                   }
1737                 printf("\n\n\n");
1738                 for(i=0; i<nSeq1; i++)
1739                   {
1740                     printf("\n");
1741                     for(j=0; j<nSeq2; j++)
1742                       {
1743                         printf("%4d", (int)(*(brk_flg[i]+j)));
1744                       }
1745                   }
1746                 // exit(0);
1747                 */
1748
1749                is=0; js=0; lcmp=0;
1750                // no nc-end penalty - begin
1751                sum_ret=mat[0][0];
1752
1753                // look for the highest score in mat[i][j]
1754                for(i=0; i<nSeq1; i++)
1755                        for(j=0; j<nSeq2; j++)
1756                        {
1757                                if(i==0 && j==0) continue;
1758                                sum=mat[i][j];
1759                                if(isGlobal1) sum+=-gapI-gapE*i;
1760                                if(isGlobal2) sum+=-gapI-gapE*j;
1761                                if(sum>sum_ret)
1762                                {
1763                                        sum_ret=sum;
1764                                        is=i; js=j;
1765                                }
1766                        }
1767
1768                //System.out.println("start at " + is + "  " + js);
1769                //for(k=0; k<is; k++) align1[k]=-1;
1770                //for(k=0; k<js; k++) align2[k]=-1;
1771                // no nc-end penalty - end
1772
1773                for(i=is, j=js; i<nSeq1 && j<nSeq2; i++, j++)
1774                {
1775                        iMax=i; jMax=j;
1776                        sum=mat[i][j];
1777                        if(!ge)
1778                        {
1779                                for(k=i+1; k<nSeq1; k++)
1780                                        if(mat[k][j]-gapI>sum)
1781                                        {
1782                                                iMax=k; jMax=j;
1783                                                sum=mat[k][j]-gapI;
1784                                        }
1785
1786                                for(k=j+1; k<nSeq2; k++)
1787                                        if(mat[i][k]-gapI>sum)
1788                                        {
1789                                                iMax=i; jMax=k;
1790                                                sum=mat[i][k]-gapI;
1791                                        }
1792                        }
1793                        else
1794                        {
1795                                for(k=i+1; k<nSeq1; k++)
1796                                        if(mat[k][j]-gapI-gapE*(k-i)>sum)
1797                                        {
1798                                                //System.out.println("gap1 " + k + " " + j + " " + sum + "<" +(mat[k][j]-gapI-gapE*(k-i)));
1799                                                iMax=k; jMax=j;
1800                                                sum=mat[k][j]-gapI-gapE*(k-i);
1801                                        }
1802
1803                                for(k=j+1; k<nSeq2; k++)
1804                                        if(mat[i][k]-gapI-gapE*(k-j)>sum)
1805                                        {
1806                                                //System.out.println("gap2 " + i + " " + k + " " + sum + "<"+ (mat[i][k]-gapI-gapE*(k-j)));
1807                                                iMax=i; jMax=k;
1808                                                sum=mat[i][k]-gapI-gapE*(k-j);
1809                                        }
1810                        }
1811
1812                        //if ( i != iMax || j != jMax )
1813                        //      System.out.println("FOUND GAP AT: " + i+ " " + iMax + " " + j + " " + jMax);
1814
1815                        //System.out.println(" iMax " + iMax + " jMax " +  jMax);
1816                        // set the gap positions:
1817                        //printf("%d %d\n", iMax, jMax);
1818
1819
1820                        for(k=i; k<iMax; k++, i++) {
1821                                align_se1[lcmp]=k;
1822                                align_se2[lcmp]=-1;
1823
1824
1825                                lcmp++;
1826                        }
1827
1828                        for(k=j; k<jMax; k++, j++) {
1829                                align_se1[lcmp]=-1;
1830                                align_se2[lcmp]=k;
1831
1832
1833                                lcmp++;
1834                        }
1835
1836
1837                        align_se1[lcmp]=iMax;
1838                        align_se2[lcmp]=jMax;
1839                        lcmp++;
1840
1841                        if(brk_flg[i][j]) {
1842                                //System.out.println("hit break flag at: " + i + "  " + j + " sum " + sum_ret + " lcmp " + lcmp);
1843                                break;
1844
1845                        }
1846                }
1847
1848
1849                return sum_ret;
1850        }
1851
1852
1853        private void rot_mol(Atom[] caA, Atom[] caB, int nse2, Matrix m , Atom shift) throws StructureException{
1854
1855                for(int l=0; l<nse2; l++) {
1856                        Atom a = caA[l];
1857                        Group g = (Group)a.getGroup().clone();
1858                        //Group g = (Group)a.getParent();
1859
1860                        Calc.rotate( g, m);
1861                        Calc.shift(  g, shift);
1862                        caB[l] = g.getAtom(a.getName());
1863                }
1864        }
1865
1866        //void rot_mol(XYZ *molA, XYZ *molB, int nAtom, double d_[20] ) {
1867        //                        double dx, dy, dz;
1868        //                        for(int l=0; l<nAtom; l++) {
1869        //                          if(molA[l].X<1e10) {
1870        //                            dx=molA[l].X;
1871        //                            dy=molA[l].Y;
1872        //                            dz=molA[l].Z;
1873        //                            molB[l].X=dx*d_[0]+dy*d_[1]+dz*d_[2]+d_[9];
1874        //                            molB[l].Y=dx*d_[3]+dy*d_[4]+dz*d_[5]+d_[10];
1875        //                            molB[l].Z=dx*d_[6]+dy*d_[7]+dz*d_[8]+d_[11];
1876        //                          }
1877        //                          else {
1878        //                            molB[l]=molA[l];
1879        //                          }
1880        //                        }
1881        //                      }
1882        //
1883
1884        /** superimpose and get rmsd
1885         *
1886         * @param pro1
1887         * @param pro2
1888         * @param strLen Number of atoms from pro1 and pro2 to use
1889         * @param storeTransform Store rotation and shift matrices locally
1890         * @return RMSD
1891         * @throws StructureException
1892         */
1893        public double calc_rmsd(Atom[] pro1, Atom[] pro2, int strLen, 
1894                        boolean storeTransform) throws StructureException {
1895
1896                Atom[] cod1 = getAtoms(pro1,  strLen,false);
1897                Atom[] cod2 = getAtoms(pro2,  strLen,true);
1898
1899                Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(cod1), 
1900                                Calc.atomsToPoints(cod2));
1901
1902                Matrix matrix = Matrices.getRotationJAMA(trans);
1903                Atom shift = Calc.getTranslationVector(trans);
1904
1905                if ( storeTransform) {
1906                        r = matrix;
1907                        t = shift;
1908                }
1909                for (Atom a : cod2)
1910                        Calc.transform(a.getGroup(), trans);
1911                        
1912                return Calc.rmsd(cod1, cod2);
1913
1914        }
1915
1916        /**
1917         * Copies the first length atoms from the input array
1918         * @param ca The array to copy
1919         * @param length the number of atoms to copy
1920         * @param clone If true, preform a deep copy, cloning the underlying Groups
1921         * @return An array with the first length items of ca, possibly cloning the Atoms.
1922         * @throws StructureException
1923         */
1924        private Atom[] getAtoms(Atom[] ca,  int length, boolean clone) throws StructureException{
1925
1926                List<Atom> atoms = new ArrayList<Atom>();
1927                for ( int i = 0 ; i < length ; i++){
1928
1929                        Atom a;
1930                        if ( clone ){
1931                                Group g = (Group)ca[i].getGroup().clone();
1932                                a = g.getAtom(ca[i].getName());
1933                        }
1934                        else {
1935                                a = ca[i];
1936                        }
1937                        atoms.add(a);
1938                }
1939                return atoms.toArray(new Atom[atoms.size()]);
1940        }
1941
1942
1943
1944        private void noBestTrace(){
1945
1946                if(isPrint) {
1947                        timeEnd = System.currentTimeMillis();
1948                        long time_q=(timeEnd-timeStart);
1949
1950                        String msg = String.format("size=0 time=%d comb=%d\n", (int)(time_q), nTraces);
1951                        System.out.println(msg);
1952                }
1953        }
1954
1955
1956
1957        private double zToP(double z) {
1958                int ind=(int)(z/0.1);
1959                if(ind<0) ind=0;
1960                if(ind>149) ind=149;
1961                return(tableZtoP[ind]);
1962        }
1963        ///////////////////////////////////////////////////////////////////////////
1964        private         double pToZ(double p) {
1965                int ind=(int)(-Math.log10(p)*3.0);
1966                if(ind<0) ind=0;
1967                if(ind>149) ind=149;
1968                return(tablePtoZ[ind]);
1969        }
1970        ///////////////////////////////////////////////////////////////////////////
1971        private double zByZ(double z1, double z2) {
1972                double p1=zToP(z1);
1973                double p2=zToP(z2);
1974                return(pToZ(p1*p2));
1975        }
1976
1977        protected double zStrAlign(int winSize, int nTrace, double score, int nGaps) {
1978                double z1=zScore(winSize, nTrace, score);
1979                double z2=zGaps(winSize, nTrace, nGaps);
1980                return(zByZ(z1, z2));
1981        }
1982
1983        double zScore(int winSize, int nTrace, double score) {
1984
1985                if(winSize==8) {
1986
1987                        if(nTrace<1) return(0.0);
1988
1989                        double scoreAv_, scoreSd_;
1990                        if(nTrace<21) {
1991                                scoreAv_=scoreAv8[nTrace-1];
1992                                scoreSd_=scoreSd8[nTrace-1];
1993                        }
1994                        else {
1995                                scoreAv_=0.209874*nTrace+2.944714;
1996                                scoreSd_=0.039487*nTrace+0.675735;
1997                        }
1998                        if(score>scoreAv_) return(0.0);
1999                        return((scoreAv_-score)/scoreSd_);
2000                }
2001
2002                if(winSize==6) {
2003
2004                        if(nTrace<1) return(0.0);
2005
2006                        double scoreAv_, scoreSd_;
2007                        if(nTrace<21) {
2008                                scoreAv_=scoreAv6[nTrace-1];
2009                                scoreSd_=scoreSd6[nTrace-1];
2010                        }
2011                        else {
2012                                scoreAv_=0.198534*nTrace+2.636477;
2013                                scoreSd_=0.040922*nTrace+0.715636;
2014                        }
2015                        if(score>scoreAv_) return(0.0);
2016                        return((scoreAv_-score)/scoreSd_);
2017                }
2018
2019                return(0.0);
2020
2021        }
2022        ///////////////////////////////////////////////////////////////////////////
2023        double zGaps(int winSize, int nTrace, int nGaps) {
2024
2025                if(nTrace<2) return(0.0);
2026                double scoreAv_, scoreSd_;
2027
2028                if(winSize==8) {
2029                        if(nTrace<21) {
2030                                scoreAv_=gapsAv8[nTrace-1];
2031                                scoreSd_=gapsSd8[nTrace-1];
2032                        }
2033                        else {
2034                                scoreAv_=14.949173*nTrace-14.581193;
2035                                scoreSd_=2.045067*nTrace+13.191095;
2036                        }
2037                        if(nGaps>scoreAv_) return(0.0);
2038                        return((scoreAv_-nGaps)/scoreSd_);
2039                }
2040
2041                if(winSize==6) {
2042                        if(nTrace<21) {
2043                                scoreAv_=gapsAv6[nTrace-1];
2044                                scoreSd_=gapsSd6[nTrace-1];
2045                        }
2046                        else {
2047                                scoreAv_=13.574490*nTrace-13.977223;
2048                                scoreSd_=1.719977*nTrace+19.615014;
2049                        }
2050                        if(nGaps>scoreAv_) return(0.0);
2051                        return((scoreAv_-nGaps)/scoreSd_);
2052                }
2053
2054                return(0.0);
2055        }
2056
2057        private static final double[] scoreAv8 ={2.54, 2.51, 2.72, 3.01, 3.31, 3.61, 3.90, 4.19, 4.47, 4.74,
2058                4.99, 5.22, 5.46, 5.70, 5.94, 6.13, 6.36, 6.52, 6.68, 6.91};
2059        private static final double[] scoreSd8 ={1.33, 0.88, 0.73, 0.71, 0.74, 0.80, 0.86, 0.92, 0.98, 1.04,
2060                1.08, 1.10, 1.15, 1.19, 1.23, 1.25, 1.32, 1.34, 1.36, 1.45};
2061        private static final double[] gapsAv8 ={0.00, 11.50, 23.32, 35.95, 49.02, 62.44, 76.28, 90.26,
2062                104.86, 119.97, 134.86, 150.54, 164.86, 179.57, 194.39,
2063                209.38, 224.74, 238.96, 253.72, 270.79};
2064        private static final double[] gapsSd8 ={0.00, 9.88, 14.34, 17.99, 21.10, 23.89, 26.55, 29.00, 31.11,
2065                33.10, 35.02, 36.03, 37.19, 38.82, 41.04, 43.35, 45.45,
2066                48.41, 50.87, 52.27};
2067        private static final double[] scoreAv6 ={1.98, 1.97, 2.22, 2.54, 2.87, 3.18, 3.48, 3.77, 4.05, 4.31,
2068                4.57, 4.82, 5.03, 5.24, 5.43, 5.64, 5.82, 6.02, 6.21, 6.42};
2069        private static final double[] scoreSd6 ={1.15, 0.73, 0.63, 0.64, 0.71, 0.80, 0.87, 0.95, 1.01, 1.07,
2070                1.13, 1.19, 1.22, 1.25, 1.28, 1.32, 1.35, 1.39, 1.45, 1.50};
2071        private static final double[] gapsAv6 ={0.00, 10.12, 20.25, 31.29, 42.95, 55.20, 67.53, 80.15,
2072                93.30, 106.47, 120.52, 134.38, 148.59, 162.58, 176.64,
2073                191.23, 204.12, 218.64, 231.82, 243.43};
2074        private static final double[] gapsSd6 ={0.00, 9.80, 14.44, 18.14, 21.35, 24.37, 27.00, 29.68, 32.22,
2075                34.37, 36.65, 38.63, 40.31, 42.16, 43.78, 44.98, 47.08,
2076                49.09, 50.78, 52.15};
2077
2078
2079        private static final double[] tableZtoP={
2080                1.0, 9.20e-01,8.41e-01,7.64e-01,6.89e-01,6.17e-01,5.49e-01,4.84e-01,4.24e-01,3.68e-01,
2081                3.17e-01,2.71e-01,2.30e-01,1.94e-01,1.62e-01,1.34e-01,1.10e-01,8.91e-02,7.19e-02,5.74e-02,
2082                4.55e-02,3.57e-02,2.78e-02,2.14e-02,1.64e-02,1.24e-02,9.32e-03,6.93e-03,5.11e-03,3.73e-03,
2083                2.70e-03,1.94e-03,1.37e-03,9.67e-04,6.74e-04,4.65e-04,3.18e-04,2.16e-04,1.45e-04,9.62e-05,
2084                6.33e-05,4.13e-05,2.67e-05,1.71e-05,1.08e-05,6.80e-06,4.22e-06,2.60e-06,1.59e-06,9.58e-07,
2085                5.73e-07,3.40e-07,1.99e-07,1.16e-07,6.66e-08,3.80e-08,2.14e-08,1.20e-08,6.63e-09,3.64e-09,
2086                1.97e-09,1.06e-09,5.65e-10,2.98e-10,1.55e-10,8.03e-11,4.11e-11,2.08e-11,1.05e-11,5.20e-12,
2087                2.56e-12,1.25e-12,6.02e-13,2.88e-13,1.36e-13,6.38e-14,2.96e-14,1.36e-14,6.19e-15,2.79e-15,
2088                1.24e-15,5.50e-16,2.40e-16,1.04e-16,4.46e-17,1.90e-17,7.97e-18,3.32e-18,1.37e-18,5.58e-19,
2089                2.26e-19,9.03e-20,3.58e-20,1.40e-20,5.46e-21,2.10e-21,7.99e-22,3.02e-22,1.13e-22,4.16e-23,
2090                1.52e-23,5.52e-24,1.98e-24,7.05e-25,2.48e-25,8.64e-26,2.98e-26,1.02e-26,3.44e-27,1.15e-27,
2091                3.82e-28,1.25e-28,4.08e-29,1.31e-29,4.18e-30,1.32e-30,4.12e-31,1.27e-31,3.90e-32,1.18e-32,
2092                3.55e-33,1.06e-33,3.11e-34,9.06e-35,2.61e-35,7.47e-36,2.11e-36,5.91e-37,1.64e-37,4.50e-38,
2093                1.22e-38,3.29e-39,8.77e-40,2.31e-40,6.05e-41,1.56e-41,4.00e-42,1.02e-42,2.55e-43,6.33e-44,
2094                1.56e-44,3.80e-45,9.16e-46,2.19e-46,5.17e-47,1.21e-47,2.81e-48,6.45e-49,1.46e-49,3.30e-50};
2095        private static final double[] tablePtoZ ={
2096                0.00,0.73,1.24,1.64,1.99,2.30,2.58,2.83,3.07,3.29,
2097                3.50,3.70,3.89,4.07,4.25,4.42,4.58,4.74,4.89,5.04,
2098                5.19,5.33,5.46,5.60,5.73,5.86,5.99,6.11,6.23,6.35,
2099                6.47,6.58,6.70,6.81,6.92,7.02,7.13,7.24,7.34,7.44,
2100                7.54,7.64,7.74,7.84,7.93,8.03,8.12,8.21,8.30,8.40,
2101                8.49,8.57,8.66,8.75,8.84,8.92,9.01,9.09,9.17,9.25,
2102                9.34,9.42,9.50,9.58,9.66,9.73,9.81,9.89,9.97,10.04,
2103                10.12,10.19,10.27,10.34,10.41,10.49,10.56,10.63,10.70,10.77,
2104                10.84,10.91,10.98,11.05,11.12,11.19,11.26,11.32,11.39,11.46,
2105                11.52,11.59,11.66,11.72,11.79,11.85,11.91,11.98,12.04,12.10,
2106                12.17,12.23,12.29,12.35,12.42,12.48,12.54,12.60,12.66,12.72,
2107                12.78,12.84,12.90,12.96,13.02,13.07,13.13,13.19,13.25,13.31,
2108                13.36,13.42,13.48,13.53,13.59,13.65,13.70,13.76,13.81,13.87,
2109                13.92,13.98,14.03,14.09,14.14,14.19,14.25,14.30,14.35,14.41,
2110                14.46,14.51,14.57,14.62,14.67,14.72,14.77,14.83,14.88,14.93};
2111
2112        /** copy data from this class into AFPChain container object.
2113         *
2114         * @param afpChain
2115         * @param ca1
2116         * @param ca2
2117         */
2118         public void convertAfpChain(AFPChain afpChain, Atom[] ca1, Atom[] ca2) {
2119
2120                 afpChain.setBlockNum(1);
2121                 //afpChain.setAlignScore(z);
2122                 Matrix[] m ;
2123
2124                 if ( r != null ) {
2125                         m = new Matrix[1];
2126                         m[0] = r;
2127                 } else  {
2128                         m = new Matrix[0];
2129                 }
2130
2131                 Atom[] as ;
2132                 if ( t != null) {
2133                         as = new Atom[1];
2134                         as[0] = t;
2135                 } else {
2136                         as = new Atom[0];
2137                 }
2138
2139                 afpChain.setBlockRotationMatrix(m);
2140                 afpChain.setBlockShiftVector(as);
2141
2142                 int nse1 = ca1.length;
2143                 int nse2 = ca2.length;
2144                 //System.out.println("dist1 :" + dist1.length + " " + dist2.length);
2145
2146                 if ( nse1 > 0 && dist1.length > 0 )
2147                         afpChain.setDisTable1(new Matrix(dist1));
2148                 else
2149                         afpChain.setDisTable1 (Matrix.identity(3, 3));
2150                 if ( nse2 > 0 && dist2.length > 0 )
2151                         afpChain.setDisTable2(new Matrix(dist2));
2152                 else
2153                         afpChain.setDisTable2(Matrix.identity(3, 3));
2154
2155
2156                 char[] alnseq1 = new char[nse1+nse2+1];
2157                 char[] alnseq2 = new char[nse1+nse2+1] ;
2158                 char[] alnsymb = new char[nse1+nse2+1];
2159
2160                 int[][][] optAln = new int[1][2][nAtom];
2161                 afpChain.setOptAln(optAln);
2162
2163                 int pos = 0;
2164                 int nrIdent = 0;
2165                 int nrSim   = 0;
2166                 for(int ia=0; ia<lcmp; ia++) {
2167
2168                         // no gap
2169                         if(align_se1[ia]!=-1 && align_se2[ia]!=-1) {
2170                                 //System.out.println("ia " + ia + " pos " + pos + " "  + align_se1[ia] + " " + align_se2[ia]);
2171                                 optAln[0][0][pos] = align_se1[ia];
2172                                 optAln[0][1][pos] = align_se2[ia];
2173
2174                                 char l1 = getOneLetter(ca1[align_se1[ia]].getGroup());
2175                                 char l2 = getOneLetter(ca2[align_se2[ia]].getGroup());
2176
2177                                 alnseq1[ia] = Character.toUpperCase(l1);
2178                                 alnseq2[ia] = Character.toUpperCase(l2);
2179                                 alnsymb[ia] = ' ';
2180                                 if ( l1 == l2) {
2181                                         nrIdent++;
2182                                         nrSim++;
2183                                         alnsymb[ia] = '|';
2184                                 } else if ( AFPAlignmentDisplay.aaScore(l1, l2) > 0){
2185                                         nrSim++;
2186                                         alnsymb[ia] = ':';
2187                                 }
2188
2189                                 pos++;
2190
2191                         } else {
2192                                 // there is a gap at this position
2193                                 alnsymb[ia] = ' ';
2194                                 if (align_se1[ia] == -1 ) {
2195                                         alnseq1[ia] = '-';
2196                                 } else {
2197                                         char l1 = getOneLetter(ca1[align_se1[ia]].getGroup());
2198                                         alnseq1[ia] = Character.toUpperCase(l1);
2199                                 }
2200                                 if ( align_se2[ia] == -1 ) {
2201                                         alnseq2[ia] = '-';
2202                                 } else {
2203                                         char l2 = getOneLetter(ca2[align_se2[ia]].getGroup());
2204                                         alnseq2[ia] = Character.toUpperCase(l2);
2205                                 }
2206
2207                         }
2208                 }
2209
2210
2211                 afpChain.setAlnseq1(alnseq1);
2212                 afpChain.setAlnseq2(alnseq2);
2213                 afpChain.setAlnsymb(alnsymb);
2214
2215
2216                 // CE uses the aligned pairs as reference not the whole alignment including gaps...
2217                 if ( pos > 0) {
2218                         afpChain.setIdentity(nrIdent*1.0/pos);
2219                         afpChain.setSimilarity(nrSim*1.0/pos);
2220                 } else {
2221                         afpChain.setIdentity(0);
2222                         afpChain.setSimilarity(0);
2223                 }
2224
2225                 //AFPAlignmentDisplay.getAlign( afpChain,ca1,ca2);
2226
2227         }
2228
2229         public int getnAtom() {
2230                 return nAtom;
2231         }
2232
2233
2234
2235         public void setnAtom(int nAtom) {
2236                 this.nAtom = nAtom;
2237         }
2238
2239
2240
2241         public int getLcmp() {
2242                 return lcmp;
2243         }
2244
2245
2246
2247         public void setLcmp(int lcmp) {
2248                 this.lcmp = lcmp;
2249         }
2250
2251
2252
2253         public int[] getAlign_se1() {
2254                 return align_se1;
2255         }
2256
2257
2258
2259         public void setAlign_se1(int[] alignSe1) {
2260                 align_se1 = alignSe1;
2261         }
2262
2263
2264
2265         public int[] getAlign_se2() {
2266                 return align_se2;
2267         }
2268
2269
2270
2271         public void setAlign_se2(int[] alignSe2) {
2272                 align_se2 = alignSe2;
2273         }
2274
2275         /**
2276          * Caution: this matrix is overwriten with very different data at several
2277          * points in the alignment algorithm. After
2278          * {@link #initSumOfDistances(int, int, int, int, Atom[], Atom[]) initSumOfDistances}
2279          * is run, this will hold the distance matrix between AFPs.
2280          * @return mat
2281          */
2282         public double[][] getMatMatrix() {
2283                 return mat;
2284         }
2285
2286         public void setMatMatrix(double[][] matrix){
2287                 mat = matrix;
2288         }
2289
2290         /**
2291          * Gets the rotation matrix from the last call to
2292          * {@link #calc_rmsd(Atom[], Atom[], int, boolean) calc_rmsd}.
2293          * @return The rotatiokn matrix
2294          */
2295         public Matrix getRotationMatrix() {
2296                 return r;
2297         }
2298
2299         /**
2300          * Gets the shift from the last call to
2301          * {@link #calc_rmsd(Atom[], Atom[], int, boolean) calc_rmsd}.
2302          * @return The shift
2303          */
2304         public Atom getShift() {
2305                 return t;
2306         }
2307
2308        public double[][] getDist1() {
2309                return dist1;
2310        }
2311
2312        public void setDist1(double[][] dist1) {
2313                this.dist1 = dist1;
2314        }
2315
2316        public double[][] getDist2() {
2317                return dist2;
2318        }
2319
2320        public void setDist2(double[][] dist2) {
2321                this.dist2 = dist2;
2322        }
2323
2324
2325}