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 */
021
022
023package org.biojava.bio.dp.onehead;
024
025import java.io.Serializable;
026import java.util.ArrayList;
027import java.util.HashMap;
028import java.util.List;
029import java.util.Map;
030
031import org.biojava.bio.BioError;
032import org.biojava.bio.BioException;
033import org.biojava.bio.alignment.Alignment;
034import org.biojava.bio.alignment.SimpleAlignment;
035import org.biojava.bio.dist.Distribution;
036import org.biojava.bio.dp.BackPointer;
037import org.biojava.bio.dp.DP;
038import org.biojava.bio.dp.DPMatrix;
039import org.biojava.bio.dp.DotState;
040import org.biojava.bio.dp.EmissionState;
041import org.biojava.bio.dp.IllegalTransitionException;
042import org.biojava.bio.dp.MagicalState;
043import org.biojava.bio.dp.MarkovModel;
044import org.biojava.bio.dp.ScoreType;
045import org.biojava.bio.dp.SimpleStatePath;
046import org.biojava.bio.dp.State;
047import org.biojava.bio.dp.StatePath;
048import org.biojava.bio.symbol.AlphabetManager;
049import org.biojava.bio.symbol.DoubleAlphabet;
050import org.biojava.bio.symbol.GappedSymbolList;
051import org.biojava.bio.symbol.IllegalAlphabetException;
052import org.biojava.bio.symbol.IllegalSymbolException;
053import org.biojava.bio.symbol.SimpleGappedSymbolList;
054import org.biojava.bio.symbol.SimpleSymbolList;
055import org.biojava.bio.symbol.Symbol;
056import org.biojava.bio.symbol.SymbolList;
057
058/**
059 * An implementation of DP that aligns a single sequence against a single model.
060 *
061 * @author Matthew Pocock
062 * @author Thomas Down
063 * @author Samiul Hasan
064 * @author Lukas Kall
065 */
066public class SingleDP extends DP implements Serializable {
067  protected final HashMap emissionsProb;
068  protected final HashMap emissionsOdds;
069  protected final HashMap emissionsNull;
070
071  public SingleDP(MarkovModel model)
072  throws IllegalSymbolException, IllegalTransitionException, BioException {
073    super(model);
074    emissionsProb = new HashMap();
075    emissionsOdds = new HashMap();
076    emissionsNull = new HashMap();
077  }
078
079  public void update() {
080    // System.out.println("Updating emissions as underlying model has changed!");
081    super.update();
082    // workaround for bug in vm
083    if(emissionsProb != null) {
084      emissionsProb.clear();
085    }
086    if(emissionsOdds != null) {
087      emissionsOdds.clear();
088    }
089    if(emissionsNull != null) {
090      emissionsNull.clear();
091    }
092  }
093  
094    /**
095     * This method is public for the benefit of training algorithms,
096     * and in the future we should look at a better way of exposing
097     * the emissions cache.
098     */
099
100  public double [] getEmission(Symbol sym, ScoreType scoreType)
101  throws IllegalSymbolException {
102    Map emissions;
103    if(scoreType == ScoreType.PROBABILITY) {
104      emissions = emissionsProb;
105    } else if(scoreType == ScoreType.ODDS) {
106      emissions = emissionsOdds;
107    } else if(scoreType == ScoreType.NULL_MODEL) {
108      emissions = emissionsNull;
109    } else {
110      throw new BioError("Unknown ScoreType object: " + scoreType);
111    }
112    double [] em = (double []) emissions.get(sym);
113    if(em == null) {
114      int dsi = getDotStatesIndex();
115      em = new double[dsi];
116      State [] states = getStates();
117      if(sym == AlphabetManager.getGapSymbol()) {
118        em[0] = 0;
119      } else {
120        em[0] = Double.NEGATIVE_INFINITY;
121      }
122      for(int i = 1; i < dsi; i++) {
123        EmissionState es = (EmissionState) states[i];
124        Distribution dis = es.getDistribution();
125        em[i] = Math.log(scoreType.calculateScore(dis, sym));
126      }
127      emissions.put(sym, em);
128      /*System.out.println("Emissions for " + sym);
129      for(int i = 0; i < em.length; i++) {
130        System.out.println("\t" + states[i] + "\t-> " + em[i]);
131      }*/
132    }
133    return em;
134  }
135  
136  public double forward(SymbolList [] seq, ScoreType scoreType)
137  throws IllegalSymbolException, IllegalAlphabetException, IllegalSymbolException {
138    if(seq.length != 1) {
139      throw new IllegalArgumentException("seq must be 1 long, not " + seq.length);
140    }
141    lockModel();
142    DPCursor dpCursor = new SmallCursor(
143      getStates(),
144      seq[0],
145      seq[0].iterator()
146    );
147    double score = forward(dpCursor, scoreType);
148    unlockModel();
149
150    return score;
151  }
152  
153  public double backward(SymbolList [] seq, ScoreType scoreType)
154  throws IllegalSymbolException, IllegalAlphabetException, IllegalSymbolException {
155    if(seq.length != 1) {
156      throw new IllegalArgumentException("seq must be 1 long, not " + seq.length);
157    }
158    lockModel();
159    DPCursor dpCursor = new SmallCursor(
160      getStates(),
161      seq[0],
162      new ReverseIterator(seq[0])
163    );
164    double score = backward(dpCursor, scoreType);
165    unlockModel();
166    
167    return score;
168  }
169
170  public DPMatrix forwardMatrix(SymbolList [] seq, ScoreType scoreType)
171  throws IllegalSymbolException, IllegalAlphabetException, IllegalSymbolException {
172    if(seq.length != 1) {
173      throw new IllegalArgumentException("seq must be 1 long, not " + seq.length);
174    }
175    
176    lockModel();
177    SingleDPMatrix matrix = new SingleDPMatrix(this, seq[0]);
178    DPCursor dpCursor = new MatrixCursor(matrix, seq[0].iterator(), +1);
179    matrix.setScore(forward(dpCursor, scoreType));
180    unlockModel();
181    
182    return matrix;
183  }
184  
185  public DPMatrix backwardMatrix(SymbolList [] seq, ScoreType scoreType)
186  throws IllegalSymbolException, IllegalAlphabetException, IllegalSymbolException {
187    if(seq.length != 1) {
188      throw new IllegalArgumentException("seq must be 1 long, not " + seq.length);
189    }
190    
191    lockModel();
192    SingleDPMatrix matrix = new SingleDPMatrix(this, seq[0]);
193    DPCursor dpCursor = new MatrixCursor(matrix, new ReverseIterator(seq[0]), -1);
194    matrix.setScore(backward(dpCursor, scoreType));
195    unlockModel();
196    
197    return matrix;
198  }
199  
200  public DPMatrix forwardMatrix(SymbolList [] seq, DPMatrix matrix, ScoreType scoreType)
201  throws IllegalArgumentException, IllegalSymbolException,
202  IllegalAlphabetException, IllegalSymbolException {
203    if(seq.length != 1) {
204      throw new IllegalArgumentException("seq must be 1 long, not " + seq.length);
205    }
206    
207    lockModel();
208    SingleDPMatrix sm = (SingleDPMatrix) matrix;
209    DPCursor dpCursor = new MatrixCursor(sm, seq[0].iterator(), +1);
210    sm.setScore(forward(dpCursor, scoreType));
211    unlockModel();
212    
213    return sm;
214  }
215
216  public DPMatrix backwardMatrix(SymbolList [] seq, DPMatrix matrix, ScoreType scoreType)
217  throws IllegalArgumentException, IllegalSymbolException,
218  IllegalAlphabetException, IllegalSymbolException {
219    if(seq.length != 1) {
220      throw new IllegalArgumentException("seq must be 1 long, not " + seq.length);
221    }
222    
223    lockModel();
224    SingleDPMatrix sm = (SingleDPMatrix) matrix;
225    DPCursor dpCursor = new MatrixCursor(sm, new ReverseIterator(seq[0]), -1);
226    sm.setScore(backward(dpCursor, scoreType));
227    unlockModel();
228    
229    return sm;
230  }
231
232  protected double forward(DPCursor dpCursor, ScoreType scoreType)
233  throws IllegalSymbolException {
234    forward_initialize(dpCursor, scoreType);
235    forward_recurse(dpCursor, scoreType);
236    return forward_termination(dpCursor, scoreType);
237  }
238
239  protected double backward(DPCursor dpCursor, ScoreType scoreType)
240  throws IllegalSymbolException {
241    backward_initialize(dpCursor, scoreType);
242    backward_recurse(dpCursor, scoreType);
243    return backward_termination(dpCursor, scoreType);
244  }
245
246  protected void forward_initialize(DPCursor dpCursor, ScoreType scoreType)
247    throws IllegalSymbolException {
248    double [] v = dpCursor.currentCol();
249    State [] states = getStates();
250    
251    for (int l = 0; l < getDotStatesIndex(); l++) {
252      if(states[l] == getModel().magicalState()) {
253        //prob 1
254        v[l] = 0.0;
255      } else {
256        //prob 0
257        v[l] = Double.NEGATIVE_INFINITY;
258      }
259    }
260    
261    int [][] transitions = getForwardTransitions();
262    double [][] transitionScore = getForwardTransitionScores(scoreType);
263    double [] currentCol = dpCursor.currentCol();
264    //l over dots
265    for (int l = getDotStatesIndex(); l < states.length; l++) {
266      double score = 0.0;
267      int [] tr = transitions[l];
268      double [] trs = transitionScore[l];
269        
270      int ci = 0;
271      while(
272        ci < tr.length  && (
273        currentCol[tr[ci]] == Double.NEGATIVE_INFINITY ||
274        currentCol[tr[ci]] == Double.NaN ||
275        currentCol[tr[ci]] == Double.POSITIVE_INFINITY
276        )
277      ) {
278        ci++;
279      }
280      double constant = (ci < tr.length) ? currentCol[tr[ci]] : 0.0;
281        
282      for(int kc = 0; kc < tr.length; kc++) {
283        int k = tr[kc];
284
285        if(
286          currentCol[k] != Double.NEGATIVE_INFINITY &&
287          currentCol[k] != Double.NaN &&
288          currentCol[k] != Double.POSITIVE_INFINITY
289        ) {
290          double t = trs[kc];
291          score += Math.exp(t + (currentCol[k] - constant));
292        } else {
293        }
294      }
295      currentCol[l] = Math.log(score) + constant;
296    }
297  }
298
299  protected void backward_initialize(DPCursor dpCursor, ScoreType scoreType)
300    throws IllegalSymbolException {
301    double [] v = dpCursor.currentCol();
302    State [] states = getStates();
303
304    for (int l = 0; l < states.length; l++) {
305      if(states[l] == getModel().magicalState()) {
306        v[l] = 0.0;
307      } else {
308        v[l] = Double.NEGATIVE_INFINITY;
309      }
310    }
311  }
312
313  private void forward_recurse(DPCursor dpCursor, ScoreType scoreType)
314    throws IllegalSymbolException {
315    State [] states = getStates();
316    int [][] transitions = getForwardTransitions();
317    double [][] transitionScore = getForwardTransitionScores(scoreType);
318
319    // int _index = 0;
320    while (dpCursor.canAdvance()) {
321      // _index++;
322      // System.out.println("\n*** Index=" + _index + " ***");
323      dpCursor.advance();
324      Symbol sym = dpCursor.currentRes();
325      double [] emissions = getEmission(sym, scoreType);
326//      System.out.println("Consuming " + sym.getName());
327      double [] currentCol = dpCursor.currentCol();
328      double [] lastCol = dpCursor.lastCol();
329      for (int l = 0; l < getDotStatesIndex(); l++) { //any -> emission
330        double weight = emissions[l];
331        if (weight == Double.NEGATIVE_INFINITY) {
332          // System.out.println("*");
333          currentCol[l] = Double.NEGATIVE_INFINITY;
334        } else {
335          double score = 0.0;
336          int [] tr = transitions[l];
337          double [] trs = transitionScore[l];
338          // System.out.println("l=" + states[l].getName());
339          int ci = 0;
340          while (
341            ci < tr.length &&
342            (lastCol[tr[ci]] == Double.NEGATIVE_INFINITY
343            || lastCol[tr[ci]] == Double.NaN
344            || lastCol[tr[ci]] == Double.POSITIVE_INFINITY)
345          ) {
346            ci++;
347          }
348          double constant = (ci < tr.length) ? lastCol[tr[ci]] : 0.0;
349
350          for (int kc = 0; kc < tr.length; kc++) {
351            int k = tr[kc];
352            // System.out.println("k=" + states[k].getName());
353            if (lastCol[k] != Double.NEGATIVE_INFINITY) {
354              double t = trs[kc];
355              
356              if(states[l]== getModel().magicalState())  {
357                  // System.out.print("magic " + "lastCol[k]=" + lastCol[k] + " , ");
358                  // System.out.println("t=" + t);
359              }
360              
361              score += Math.exp(t + (lastCol[k] - constant));
362            } else {
363              // System.out.println("-");
364            }
365          }
366          // new_l = emission_l(sym) * sum_k(transition(k, l) * old_k)
367          currentCol[l] = (weight + Math.log(score)) + constant;
368          
369          // System.out.println("currentCol[" + states[l].getName() + "]=" + currentCol[l]);
370
371          if(states[l] == getModel().magicalState())  {
372              // System.out.print("magic\t");
373              //System.out.print("Weight " + weight + "\t");
374              // System.out.print(", score " + score + " = " + Math.log(score) + "\t");
375              // System.out.println(", constant " + constant);             
376          }
377        }
378      }
379      for(int l = getDotStatesIndex(); l < states.length; l++) { // all dot states from emissions
380        double score = 0.0;
381        int [] tr = transitions[l];
382        double [] trs = transitionScore[l];
383        
384        int ci = 0;
385        while(
386          ci < tr.length  && (
387          currentCol[tr[ci]] == Double.NEGATIVE_INFINITY
388          || currentCol[tr[ci]] == Double.NaN
389          || currentCol[tr[ci]] == Double.POSITIVE_INFINITY)
390        ) {
391          ci++;
392        }
393        double constant = (ci < tr.length) ? currentCol[tr[ci]] : 0.0;
394        //System.out.println("constant: " + constant);
395        //System.out.println("read from state: " + ((ci < tr.length) ? states[tr[ci]].getName() : "none"));
396        for(int kc = 0; kc < tr.length; kc++) {
397          int k = tr[kc];
398
399          if(currentCol[k] != Double.NEGATIVE_INFINITY 
400          && currentCol[k] !=Double.NaN
401          && currentCol[k] != Double.POSITIVE_INFINITY) {
402            double t = trs[kc];
403            score += Math.exp(t + (currentCol[k] - constant));
404          } else {
405          }
406        }
407        currentCol[l] = Math.log(score) + constant;
408        //System.out.println("currentCol[" + states[l].getName() + "]=" + currentCol[l]);
409      }
410    }
411  }
412
413  protected void backward_recurse(DPCursor dpCursor, ScoreType scoreType)
414    throws IllegalSymbolException {
415    State [] states = getStates();
416    int stateCount = states.length;
417    int [][] transitions = getBackwardTransitions();
418    double [][] transitionScore = getBackwardTransitionScores(scoreType);
419    double [] prevScores = new double[getDotStatesIndex()];
420
421    while (dpCursor.canAdvance()) {
422      dpCursor.advance();
423      Symbol sym = dpCursor.lastRes();
424      double [] emissions = getEmission(sym, scoreType);
425      double [] currentCol = dpCursor.currentCol();
426      double [] lastCol = dpCursor.lastCol();
427      for(int k = getDotStatesIndex() - 1; k >= 0; k--) {
428        prevScores[k] = emissions[k];
429      }
430      
431//System.out.println(sym.getName());
432      for (int k = stateCount-1; k >= 0; k--) {
433//System.out.println("State " + k + " of " + stateCount + ", " + transitions.length);
434//System.out.println(states[k].getName());
435        int [] tr = transitions[k];
436        double [] trs = transitionScore[k];
437        double score = 0.0;
438        int ci = 0;
439        while (
440          ci < tr.length &&
441          lastCol[tr[ci]] == Double.NEGATIVE_INFINITY
442        ) {
443          ci++;
444        }
445        double constant = (ci < tr.length) ? lastCol[tr[ci]] : 0.0;
446//System.out.println("Chosen constant: " + constant);
447        for (int lc = tr.length-1; lc >= 0; lc--) { // any->emission
448          int l = tr[lc];
449          if(l >= getDotStatesIndex()) {
450            continue;
451          }
452//System.out.println(states[k].getName() + " -> " + states[l].getName());
453          double weight = prevScores[l];
454//System.out.println("weight = " + weight);
455          if (
456            lastCol[l] != Double.NEGATIVE_INFINITY &&
457            weight != Double.NEGATIVE_INFINITY
458          ) {
459            double t = trs[lc];
460            score += Math.exp(t + weight + (lastCol[l] - constant));
461          }
462        }
463//System.out.println("Score = " + score);
464        for(int lc = tr.length-1; lc >= 0; lc--) { // any->dot
465          int l = tr[lc];
466          if(l < getDotStatesIndex() || l <= k) {
467            break;
468          }
469          /*System.out.println(
470            "Processing dot-state transition " +
471            states[k].getName() + " -> " + states[l].getName()
472          );*/
473          if(currentCol[l] != Double.NEGATIVE_INFINITY) {
474            score += Math.exp(trs[lc] + (currentCol[l] - constant));
475          }
476        }
477//System.out.println("Score = " + score);
478        currentCol[k] = Math.log(score) + constant;
479//System.out.println("currentCol = " + currentCol[k]);
480      }
481    }
482  }
483
484  private double forward_termination(DPCursor dpCursor, ScoreType scoreType)
485    throws IllegalSymbolException {
486    double [] scores = dpCursor.currentCol();
487    State [] states = getStates();
488
489    int l = 0;
490    while (states[l] != getModel().magicalState())
491      l++;
492
493    return scores[l];
494  }
495
496  protected double backward_termination(DPCursor dpCursor, ScoreType scoreType)
497    throws IllegalSymbolException {
498    double [] scores = dpCursor.currentCol();
499    State [] states = getStates();
500
501    int l = 0;
502    while (states[l] != getModel().magicalState())
503      l++;
504
505    return scores[l];
506  }
507  
508  public StatePath viterbi(SymbolList [] symList, ScoreType scoreType)
509  throws IllegalSymbolException {
510    SymbolList r = symList[0];
511    DPCursor dpCursor = new SmallCursor(getStates(), r, r.iterator());
512    return viterbi(dpCursor, scoreType);
513  }
514
515  private StatePath viterbi(DPCursor dpCursor, ScoreType scoreType)
516  throws IllegalSymbolException {
517    lockModel();
518    
519    State [] states = getStates();
520
521    int [][] transitions = getForwardTransitions();
522    double [][] transitionScore = getForwardTransitionScores(scoreType);
523    int stateCount = states.length;
524
525    BackPointer [] oldPointers = new BackPointer[stateCount];
526    BackPointer [] newPointers = new BackPointer[stateCount];
527
528    // initialize
529    {
530      double [] vc = dpCursor.currentCol();
531      double [] vl = dpCursor.lastCol();
532      for (int l = 0; l < getDotStatesIndex(); l++) {
533        if(states[l] == getModel().magicalState()) {
534          //System.out.println("Initializing start state to 0.0");
535          vc[l] = vl[l] = 0.0;
536          oldPointers[l] = newPointers[l] = new BackPointer(states[l]);
537        } else {
538          vc[l] = vl[l] = Double.NEGATIVE_INFINITY;
539        }
540      }
541      for (int l = getDotStatesIndex(); l < stateCount; l++) {
542        int [] tr = transitions[l];
543        double [] trs = transitionScore[l];
544        double transProb = Double.NEGATIVE_INFINITY;
545        double trans = Double.NEGATIVE_INFINITY;
546        int prev = -1;
547        for (int kc = 0; kc < tr.length; kc++) {
548          int k = tr[kc];
549          double t = trs[kc];
550          double s = vc[k];
551          double p = t + s;
552          if (p > transProb) {
553            transProb = p;
554            prev = k;
555            trans = t;
556          }
557        }
558        if(prev != -1) {
559          vc[l] = vl[l] = transProb;
560          oldPointers[l] = newPointers[l] = new BackPointer(
561            states[l],
562            newPointers[prev],
563            trans
564          );
565        } else {
566          vc [l] = vl[l] = Double.NEGATIVE_INFINITY;
567          oldPointers[l] = newPointers[l] = null;
568        }
569      }          
570    }
571
572    // viterbi
573    while (dpCursor.canAdvance()) { // symbol i
574      dpCursor.advance();
575      Symbol sym = dpCursor.currentRes();
576      double [] emissions = getEmission(sym, scoreType);
577      //System.out.println(sym.getName());
578      double [] currentCol = dpCursor.currentCol();
579      double [] lastCol = dpCursor.lastCol();
580      for (int l = 0; l < states.length; l++) { // don't move from magical state
581        double emission;
582        if(l < getDotStatesIndex()) {
583          emission = emissions[l];
584        } else {
585          emission = 0.0;
586        }
587        int [] tr = transitions[l];
588        //System.out.println("Considering " + tr.length + " alternatives");
589        double [] trs = transitionScore[l];
590        if (emission == Double.NEGATIVE_INFINITY) {
591          //System.out.println(states[l].getName() + ": impossible emission");
592          currentCol[l] = Double.NEGATIVE_INFINITY;
593          newPointers[l] = null;
594        } else {
595          double transProb = Double.NEGATIVE_INFINITY;
596          double trans = Double.NEGATIVE_INFINITY;
597          int prev = -1;
598          for (int kc = 0; kc < tr.length; kc++) {
599            int k = tr[kc];
600            double t = trs[kc];
601            double s = (l < getDotStatesIndex()) ? lastCol[k] : currentCol[k];
602            double p = t + s;
603            /*System.out.println("Looking at scores from " + states[k].getName());
604            System.out.println("Old = " + lastCol[k]);
605            System.out.println("New = " + currentCol[k]);
606            System.out.println(
607              "Considering " + states[k].getName() + " -> " +
608              states[l].getName() + ", " + t + " + " + s + " = " + p
609            );*/
610            if (p > transProb) {
611              transProb = p;
612              prev = k;
613              trans = t;
614            }
615          }
616          if(prev != -1) {
617            currentCol[l] = transProb + emission;
618            /*System.out.println(
619              states[prev].getName() + "->" + states[l].getName() + ", " +
620              (trans + emission)
621            );*/
622            newPointers[l] = new BackPointer(
623              states[l],
624              (l < getDotStatesIndex()) ? oldPointers[prev] : newPointers[prev],
625              trans + emission
626            );
627            /*System.out.println("Succesfully completed " + states[l].getName());
628            System.out.println("Old = " + lastCol[l]);
629            System.out.println("New = " + currentCol[l]);*/
630          } else {
631            //System.out.println(states[l].getName() + ": Nowhere to come from");
632            currentCol[l] = Double.NEGATIVE_INFINITY;
633            newPointers[l] = null;
634          }
635        }
636      }
637      
638      BackPointer [] bp = newPointers;
639      newPointers = oldPointers;
640      oldPointers = bp;
641    }
642
643    // find max in last row
644    BackPointer best = null;
645    double bestScore = Double.NaN;
646    for (int l = 0; l < stateCount; l++) {
647      if (states[l] == getModel().magicalState()) {
648        best = oldPointers[l].back;
649        bestScore = dpCursor.currentCol()[l];
650        break;
651      }
652    }
653
654    int len = 0;
655    BackPointer b2 = best;
656    int dotC = 0;
657    int emC = 0;
658    // trace back ruit to check out size of path
659    while(b2.back != b2) {
660      len++;
661      if(b2.state instanceof EmissionState) {
662        emC++;
663      } else {
664        dotC++;
665      }
666      b2 = b2.back;
667    };
668
669    Map aMap = new HashMap();
670    aMap.put(dpCursor.symList(), dpCursor.symList());
671    Alignment ali = new SimpleAlignment(aMap);
672    GappedSymbolList symView = new SimpleGappedSymbolList(ali);
673    double [] scores = new double[len];
674    List stateList = new ArrayList(len);
675    for (int j = 0; j < len; j++) {
676      stateList.add(null);
677    }
678
679    b2 = best;
680    int ri = dpCursor.symList().length()+1;
681    int lc = len;
682    int gaps = 0;
683    while(b2.back != b2) {
684      lc--;
685      //System.out.println("At " + lc + " state=" + b2.state.getName() + ", score=" + b2.score + ", back=" + b2.back);
686      if(b2.state instanceof MagicalState) {
687        b2 = b2.back;
688        continue;
689      }
690      stateList.set(lc, b2.state);
691      if(b2.state instanceof DotState) {
692        symView.addGapInSource(ri);
693        gaps++;
694      } else {
695        ri--;
696      }
697      scores[lc] = b2.score;
698      b2 = b2.back;
699    }
700
701    /*System.out.println("Counted " + emC + " emissions and " + dotC + " dots");
702    System.out.println("Counted backpointers. Alignment of length " + len);
703    System.out.println("Counted states " + stateList.size());
704    System.out.println("Input list had length " + dpCursor.symList().length());
705    System.out.println("Added gaps: " + gaps);
706    System.out.println("Gapped view has length " + symView.length());*/
707
708    unlockModel();
709    return new SimpleStatePath(
710      bestScore,
711      symView,
712      new SimpleSymbolList(getModel().stateAlphabet(), stateList),
713      DoubleAlphabet.fromArray(scores)
714    );
715  }
716}