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
022package org.biojava.bio.dp.twohead;
023
024import java.util.Arrays;
025import java.util.List;
026
027import org.biojava.bio.dp.BackPointer;
028import org.biojava.bio.symbol.AlphabetManager;
029import org.biojava.bio.symbol.IllegalSymbolException;
030import org.biojava.bio.symbol.Symbol;
031import org.biojava.bio.symbol.SymbolList;
032
033/**
034 * A LIGHT implementation of PairDPCursor.
035 * <p>
036 * This object manages memory that is linear on the length of the shortest
037 * sequence. It does not maintain any data beyond that necessary for the next
038 * round of calcCell invocations.
039 *
040 * @author     Matthew Pocock 
041 * @author     David Huen (fixes for magical state)
042 */
043public class LightPairDPCursor implements PairDPCursor {
044  //    protected CrossProductAlphabet alpha;
045  private int[] pos;
046  private boolean flip;
047  private SymbolList[] seqs;
048  private double[][][] columns;
049  private double[][][] emissions;
050  /**
051   *  Description of the Field 
052   */
053  protected BackPointer[][][] bPointers;
054  /**
055   *  Description of the Field 
056   */
057  protected int numStates;
058  /**
059   *  Description of the Field 
060   */
061  protected double[] zeroCol;
062  /**
063   *  Description of the Field 
064   */
065  protected BackPointer[] emptyBP;
066  /**
067   *  Description of the Field 
068   */
069  protected int[] depth;
070  /**
071   *  Description of the Field 
072   */
073  protected EmissionCache eCache;
074
075
076  /**
077   *  Constructor for the LightPairDPCursor object 
078   *
079   * @param  seq1                        First sequence in this twohead DP.
080   * @param  seq2                        Second sequence in this twohead DP.
081   * @param  depth1                      The number of bases of context required in 
082   *                                     first sequence to compute DP matrix at a point (= max advance in first sequence + 1).
083   * @param  depth2                      The number of bases of context required in 
084   *                                     second sequence to compute DP matrix at a point (= max advance in second sequence + 1).
085   * @param  numStates                   Total number of states in model.
086   * @param  eCache                      Emission cache to be used with this run.
087   * @exception  IllegalSymbolException  Description of Exception 
088   */
089  public LightPairDPCursor(
090      SymbolList seq1, 
091      SymbolList seq2, 
092  //      CrossProductAlphabet alpha,
093      int depth1, 
094      int depth2, 
095      int numStates, 
096      EmissionCache eCache
097  ) throws IllegalSymbolException {
098    this.numStates = numStates;
099    //      this.alpha = alpha;
100    this.zeroCol = new double[this.numStates];
101    // don't touch this, please...
102    for (int i = 0; i < zeroCol.length; ++i) {
103      this.zeroCol[i] = Double.NaN;
104    }
105    this.emptyBP = new BackPointer[numStates];
106    this.pos = new int[2];
107    this.pos[0] = 0;
108    this.pos[1] = 0;
109    this.seqs = new SymbolList[2];
110    this.seqs[0] = seq1;
111    this.seqs[1] = seq2;
112    this.depth = new int[2];
113    this.depth[0] = depth1;
114    this.depth[1] = depth2;
115    this.eCache = eCache;
116
117    this.flip = this.seqs[1].length() > this.seqs[0].length();
118    //System.out.println("flip=" + this.flip);
119    if (flip) {
120      this.columns = 
121          new double[depth[0]][seqs[1].length() + 2][numStates];
122      this.emissions = 
123          new double[depth[0]][seqs[1].length() + 2][];
124      this.bPointers = 
125          new BackPointer[depth[0]][seqs[1].length() + 2][numStates];
126    }
127    else {
128      this.columns = 
129          new double[depth[1]][seqs[0].length() + 2][numStates];
130      this.emissions = 
131          new double[depth[1]][seqs[0].length() + 2][];
132      this.bPointers = 
133          new BackPointer[depth[1]][seqs[0].length() + 2][numStates];
134    }
135
136    for (int i = 0; i < columns.length; i++) {
137      double[][] ci = columns[i];
138      for (int j = 0; j < ci.length; j++) {
139        double[] cj = ci[j];
140        for (int k = 0; k < cj.length; k++) {
141          cj[k] = Double.NaN;
142        }
143      }
144    }
145
146    // compute the emissions vector at origin of matrix
147    calcEmissions(emissions[0]);
148  }
149
150
151  /**
152   *  Gets the Depth attribute of the LightPairDPCursor object 
153   *
154   * @return    The Depth value 
155   */
156  public int[] getDepth() {
157    return depth;
158  }
159
160
161  /**
162   *  Are there further Cells to be computed? 
163   *
164   * @return    Description of the Returned Value 
165   */
166  public boolean hasNext() {
167    int i = flip ? 0 : 1;
168    return 
169        pos[i] <= seqs[i].length() + 1;
170  }
171
172
173  /**
174   * Returns the minimal context of the DP matrix
175   * necessary to compute the value of a single point
176   * in that matrix.
177   * <p>
178   * The Cell [][] array has the origin as the
179   * current point to be evaluated and successive
180   * rows/columns represent rows/columns backwards
181   * within the DP matrix. 
182   *
183   * @return    An array representing the immediate context around the element to be computed. 
184   */
185  public Cell[][] press() {
186    Cell[][] cells = new Cell[depth[0]][depth[1]];
187    for (int i = 0; i < cells.length; i++) {
188      Cell[] ci = cells[i];
189      for (int j = 0; j < ci.length; j++) {
190        ci[j] = new Cell();
191      }
192    }
193    return cells;
194  }
195
196
197  /**
198   *  Description of the Method 
199   *
200   * @param  cells                       Description of Parameter 
201   * @exception  IllegalSymbolException  Description of Exception 
202   */
203  public void next(Cell[][] cells) throws IllegalSymbolException {
204    //System.out.println("Pos=" + pos[0] + ", " + pos[1] + " " + flip);
205    if (flip) {
206      for (int i = 0; i < depth[0]; i++) {
207        int ii = pos[0] - i;
208        boolean outI = (ii < 0) || (ii > seqs[0].length() + 1);
209        Cell[] cellI = cells[i];
210        // lucky in this case - can pre-fetch cells
211        double[][] columnsI = columns[i];
212        double[][] emissionsI = emissions[i];
213        BackPointer[][] bPointersI = bPointers[i];
214        for (int j = 0; j < depth[1]; j++) {
215          int jj = pos[1] - j;
216          boolean outJ = (jj < 0) || (jj > seqs[1].length() + 1);
217          //System.out.println("at " + i + "->" + ii + ", " + j + "->" + jj);
218          Cell c = cellI[j];
219          if (outI || outJ) {
220            c.scores = zeroCol;
221            c.emissions = zeroCol;
222            c.backPointers = emptyBP;
223          }
224          else {
225            c.scores = columnsI[jj];
226            c.emissions = emissionsI[jj];
227            c.backPointers = bPointersI[jj];
228          }
229        }
230      }
231      if (pos[1] <= seqs[1].length()) {
232        pos[1]++;
233      }
234      else {
235        pos[1] = 0;
236        pos[0]++;
237
238        // advance arrays
239        int depth = this.depth[0];
240        double[][] tempC = columns[depth - 1];
241        double[][] tempE = emissions[depth - 1];
242        BackPointer[][] tempBP = bPointers[depth - 1];
243        for (int i = 1; i < depth; i++) {
244          columns[i] = columns[i - 1];
245          emissions[i] = emissions[i - 1];
246          bPointers[i] = bPointers[i - 1];
247        }
248        columns[0] = tempC;
249        emissions[0] = tempE;
250        bPointers[0] = tempBP;
251
252        calcEmissions(tempE);
253      }
254    }
255    else {
256      // flip == false
257      for (int i = 0; i < depth[1]; i++) {
258        int ii = pos[1] - i;
259        boolean outI = (ii < 0) || (ii > seqs[1].length() + 1);
260        //Cell[] cellI = cells[i]; // can't pre-fetch as loop wrong way
261        double[][] columnsI = columns[i];
262        double[][] emissionsI = emissions[i];
263        BackPointer[][] bPointersI = bPointers[i];
264        for (int j = 0; j < depth[0]; j++) {
265          int jj = pos[0] - j;
266          boolean outJ = (jj < 0) || (jj > seqs[0].length() + 1);
267          //System.out.println("at " + i + "->" + ii + ", " + j + "->" + jj);
268          Cell c = cells[j][i];
269          if (outI || outJ) {
270            c.scores = zeroCol;
271            c.emissions = zeroCol;
272            c.backPointers = emptyBP;
273          }
274          else {
275            c.scores = columnsI[jj];
276            c.emissions = emissionsI[jj];
277            c.backPointers = bPointersI[jj];
278          }
279        }
280      }
281      if (pos[0] <= seqs[0].length()) {
282        pos[0]++;
283      }
284      else {
285        pos[0] = 0;
286        pos[1]++;
287
288        // advance arrays
289        int depth = this.depth[1];
290        double[][] tempC = columns[depth - 1];
291        double[][] tempE = emissions[depth - 1];
292        BackPointer[][] tempBP = bPointers[depth - 1];
293        for (int i = 1; i < depth; i++) {
294          columns[i] = columns[i - 1];
295          emissions[i] = emissions[i - 1];
296          bPointers[i] = bPointers[i - 1];
297        }
298        columns[0] = tempC;
299        emissions[0] = tempE;
300        bPointers[0] = tempBP;
301
302        calcEmissions(tempE);
303      }
304    }
305  }
306
307
308  /**
309   *  Description of the Method 
310   *
311   * @param  em                          Description of Parameter 
312   * @exception  IllegalSymbolException  Description of Exception 
313   */
314  private void calcEmissions(double[][] em)
315       throws IllegalSymbolException {
316    if (flip) {
317      //System.out.println("Calculating emissions at " + pos[0] + ", " + pos[1]);
318      Symbol[] symL = new Symbol[2];
319      List symList = Arrays.asList(symL);
320      int i = pos[0];
321      symL[0] = (i < 1 || i > seqs[0].length())
322           ? AlphabetManager.getGapSymbol()
323           : seqs[0].symbolAt(i);
324      for (int j = 0; j <= seqs[1].length() + 1; j++) {
325        symL[1] = (j < 1 || j > seqs[1].length())
326             ? AlphabetManager.getGapSymbol()
327             : seqs[1].symbolAt(j);
328        em[j] = eCache.getEmissions(symList, !(( i < 1 && j < 1 ) || ((i > seqs[0].length()) && (j > seqs[1].length())) ) );
329//        System.out.println("symbol " + symL[0].getName() + ", " + symL[1].getName() + "->" + em[j] + " " + (( i < 1 && j < 1 ) || ((i > seqs[0].length()) && (j > seqs[1].length())) ) );
330      }
331    }
332    else {
333      //System.out.println("Calculating emissions at " + pos[0] + ", " + pos[1]);
334      Symbol[] symL = new Symbol[2];
335      List symList = Arrays.asList(symL);
336      int j = pos[1];
337      symL[1] = (j < 1 || j > seqs[1].length())
338           ? AlphabetManager.getGapSymbol()
339           : seqs[1].symbolAt(j);
340      for (int i = 0; i <= seqs[0].length() + 1; i++) {
341        symL[0] = (i < 1 || i > seqs[0].length())
342             ? AlphabetManager.getGapSymbol()
343             : seqs[0].symbolAt(i);
344//        System.out.println("symbol " + symL[0].getName() + ", " + symL[1].getName() + " " + (( i < 1 && j < 1 ) || ((i > seqs[0].length()) && (j > seqs[1].length())) ) );
345        em[i] = eCache.getEmissions(symList, !(( i < 1 && j < 1 ) || ((i > seqs[0].length()) && (j > seqs[1].length())) ) );
346      }
347    }
348  }
349}
350