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