Context Navigation

source: branches/autoquest-core-tasktrees-alignment/src/main/java/de/ugoe/cs/autoquest/tasktrees/alignment/algorithms/SmithWatermanRepeated.java @ 1654

Last change on this file since 1654 was 1654, checked in by rkrimmel, 10 years ago
Adding Debug output to find this freaking damn error.
File size: 10.0 KB

Line
1	package de.ugoe.cs.autoquest.tasktrees.alignment.algorithms;
2
3	import java.util.ArrayList;
4	import java.util.Iterator;
5	import java.util.LinkedList;
6
7
8	import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.SubstitutionMatrix;
9	import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.Constants;
10
11
12	public class SmithWatermanRepeated implements AlignmentAlgorithm {
13
14	/**
15	* The first input
16	*/
17	private int[] input1;
18
19	/**
20	* The second input String
21	*/
22	private int[] input2;
23
24	/**
25	* The lengths of the input
26	*/
27	private int length1, length2;
28
29	/**
30	* The score matrix. The true scores should be divided by the normalization
31	* factor.
32	*/
33	private MatrixEntry[][] matrix;
34
35
36	private ArrayList<NumberSequence> alignment;
37
38	private float scoreThreshold;
39
40	/**
41	* Substitution matrix to calculate scores
42	*/
43	private SubstitutionMatrix submat;
44
45	public SmithWatermanRepeated() {
46
47	}
48
49	/**
50	* Compute the similarity score of substitution The position of the first
51	* character is 1. A position of 0 represents a gap.
52	*
53	* @param i
54	* Position of the character in str1
55	* @param j
56	* Position of the character in str2
57	* @return Cost of substitution of the character in str1 by the one in str2
58	*/
59	private double similarity(int i, int j) {
60	return submat.getScore(input1[i - 1], input2[j - 1]);
61	}
62
63	/**
64	* Build the score matrix using dynamic programming.
65	*/
66	private void buildMatrix() {
67	if (submat.getGapPenalty() >= 0) {
68	throw new Error("Indel score must be negative");
69	}
70
71	// it's a gap
72	matrix[0][0].setScore(0);
73	matrix[0][0].setPrevious(null); // starting point
74	matrix[0][0].setXvalue(Constants.UNMATCHED_SYMBOL);
75	matrix[0][0].setYvalue(Constants.UNMATCHED_SYMBOL);
76
77	// the first column
78	for (int j = 1; j < length2; j++) {
79	matrix[0][j].setScore(0);
80	//We don't need to go back to [0][0] if we reached matrix[0][x], so just end here
81	//matrix[0][j].setPrevious(matrix[0][j-1]);
82	matrix[0][j].setPrevious(null);
83	}
84
85
86
87	for (int i = 1; i < length1 + 2; i++) {
88
89	// Formula for first row:
90	// F(i,0) = max { F(i-1,0), F(i-1,j)-T j=1,...,m
91
92	double firstRowLeftScore = matrix[i-1][0].getScore();
93	//for sequences of length 1
94	double tempMax;
95	int maxRowIndex;
96	if(length2 == 1) {
97	tempMax = matrix[i-1][0].getScore();
98	maxRowIndex = 0;
99	} else {
100	tempMax = matrix[i-1][1].getScore();
101	maxRowIndex = 1;
102	//position of the maximal score of the previous row
103
104	for(int j = 2; j <= length2;j++) {
105	if(matrix[i-1][j].getScore() > tempMax) {
106	tempMax = matrix[i-1][j].getScore();
107	maxRowIndex = j;
108	}
109	}
110
111	}
112
113	tempMax -= scoreThreshold;
114	matrix[i][0].setScore(Math.max(firstRowLeftScore, tempMax));
115	if(tempMax == matrix[i][0].getScore()){
116	matrix[i][0].setPrevious(matrix[i-1][maxRowIndex]);
117	}
118
119	if(firstRowLeftScore == matrix[i][0].getScore()) {
120	matrix[i][0].setPrevious(matrix[i-1][0]);
121	}
122
123	//The last additional score is not related to a character in the input sequence, it's the total score. Therefore we don't need to save something for it
124	//and can end here
125	if(i<length1+1) {
126	matrix[i][0].setXvalue(input1[i-1]);
127	matrix[i][0].setYvalue(Constants.UNMATCHED_SYMBOL);
128	}
129	else {
130	return;
131	}
132
133
134
135	for (int j = 1; j <= length2; j++) {
136	double diagScore = matrix[i - 1][j - 1].getScore() + similarity(i, j);
137	double upScore = matrix[i][j - 1].getScore() + submat.getGapPenalty();
138	double leftScore = matrix[i - 1][j].getScore() + submat.getGapPenalty();
139
140	matrix[i][j].setScore(Math.max(diagScore,Math.max(upScore, Math.max(leftScore,matrix[i][0].getScore()))));
141
142	// find the directions that give the maximum scores.
143	// TODO: Multiple directions are ignored, we choose the first maximum score
144	//True if we had a match
145	if (diagScore == matrix[i][j].getScore()) {
146	matrix[i][j].setPrevious(matrix[i-1][j-1]);
147	matrix[i][j].setXvalue(input1[i-1]);
148	matrix[i][j].setYvalue(input2[j-1]);
149	}
150	//true if we took an event from sequence x and not from y
151	if (leftScore == matrix[i][j].getScore()) {
152	matrix[i][j].setXvalue(input1[i-1]);
153	matrix[i][j].setYvalue(Constants.GAP_SYMBOL);
154	matrix[i][j].setPrevious(matrix[i-1][j]);
155	}
156	//true if we took an event from sequence y and not from x
157	if (upScore == matrix[i][j].getScore()) {
158	matrix[i][j].setXvalue(Constants.GAP_SYMBOL);
159	matrix[i][j].setYvalue(input2[j-1]);
160	matrix[i][j].setPrevious(matrix[i][j-1]);
161	}
162	//true if we ended a matching region
163	if (matrix[i][0].getScore() == matrix[i][j].getScore()) {
164	matrix[i][j].setPrevious(matrix[i][0]);
165	matrix[i][j].setXvalue(input1[i-1]);
166	matrix[i][j].setYvalue(Constants.UNMATCHED_SYMBOL);
167	}
168	}
169
170	//Set the complete score cell
171
172	}
173	}
174
175	/**
176	* Get the maximum value in the score matrix.
177	*/
178	public double getMaxScore() {
179	double maxScore = 0;
180
181	// skip the first row and column
182	for (int i = 1; i <= length1; i++) {
183	for (int j = 1; j <= length2; j++) {
184	if (matrix[i][j].getScore() > maxScore) {
185	maxScore = matrix[i][j].getScore();
186	}
187	}
188	}
189
190	return maxScore;
191	}
192
193	/* (non-Javadoc)
194	* @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignmentScore()
195	*/
196	@Override
197	public double getAlignmentScore() {
198	return matrix[length1+1][0].getScore();
199	}
200
201	public void traceback() {
202	MatrixEntry tmp = matrix[length1+1][0].getPrevious();
203	LinkedList<Integer> aligned1 = new LinkedList<Integer>();
204	LinkedList<Integer> aligned2 = new LinkedList<Integer>();
205	while (tmp.getPrevious() != null) {
206
207	aligned1.add(new Integer(tmp.getXvalue()));
208	aligned2.add(new Integer(tmp.getYvalue()));
209
210	tmp = tmp.getPrevious();
211	}
212
213	// reverse order of the alignment
214	int reversed1[] = new int[aligned1.size()];
215	int reversed2[] = new int[aligned2.size()];
216
217	int count = 0;
218	for (Iterator<Integer> it = aligned1.iterator(); it.hasNext();) {
219	count++;
220	reversed1[reversed1.length - count] = it.next();
221
222	}
223	count = 0;
224	for (Iterator<Integer> it = aligned2.iterator(); it.hasNext();) {
225	count++;
226	reversed2[reversed2.length - count] = it.next();
227	}
228
229
230	NumberSequence ns1 = new NumberSequence(reversed1.length);
231	NumberSequence ns2 = new NumberSequence(reversed2.length);
232	ns1.setSequence(reversed1);
233	ns2.setSequence(reversed2);
234	ns1.setId(alignment.get(0).getId());
235	ns2.setId(alignment.get(1).getId());
236
237	alignment.set(0, ns1);
238	alignment.set(1, ns2);
239	}
240
241
242
243	public void printAlignment() {
244	int[] tmp1 = alignment.get(0).getSequence();
245	int[] tmp2 = alignment.get(1).getSequence();
246	for (int i=0; i< tmp1.length;i++) {
247	if(tmp1[i] == Constants.GAP_SYMBOL) {
248	System.out.print(" ___");
249	}
250	else if(tmp1[i] == Constants.UNMATCHED_SYMBOL) {
251	System.out.print(" ...");
252	}
253	else {
254	System.out.format("%5d", tmp1[i]);
255	}
256
257	}
258	System.out.println();
259	for (int i=0; i< tmp2.length;i++) {
260	if(tmp2[i] == Constants.GAP_SYMBOL) {
261	System.out.print(" ___");
262	}
263	else if(tmp2[i] == Constants.UNMATCHED_SYMBOL) {
264	System.out.print(" ...");
265	}
266	else {
267	System.out.format("%5d", tmp2[i]);
268	}
269
270	}
271	System.out.println();
272
273
274
275	}
276
277	public ArrayList<Match> getMatches() {
278	ArrayList<Match> result = new ArrayList<Match>();
279
280	//both alignment sequences should be equally long
281	int i = 0;
282	int[] seq1 = alignment.get(0).getSequence();
283	int[] seq2 = alignment.get(1).getSequence();
284	int start = 0;
285	while (i < seq1.length){
286	if(seq2[i] != Constants.UNMATCHED_SYMBOL) {
287	start = i;
288	int count = 0;
289	while(i < seq2.length && seq2[i] != Constants.UNMATCHED_SYMBOL) {
290	i++;
291	count++;
292	}
293	//I am really missing memcpy here? How does one do this better in java?
294	int[] tmp1 = new int[count];
295	int[] tmp2 = new int[count];
296	for (int j = 0; j<count;j++) {
297	tmp1[j] = seq1[start+j];
298	tmp2[j] = seq2[start+j];
299	}
300	NumberSequence tmpns1 = new NumberSequence(count);
301	NumberSequence tmpns2 = new NumberSequence(count);
302	tmpns1.setSequence(tmp1);
303	tmpns2.setSequence(tmp2);
304	Match tmpal = new Match();
305	tmpal.setFirstSequence(tmpns1);
306	tmpal.setSecondSequence(tmpns2);
307	//tmpal.addOccurence(new MatchOccurence(start,alignment.get(0).getId()));
308	//tmpal.addOccurence(new MatchOccurence(start,alignment.get(1).getId()));
309	result.add(tmpal);
310	}
311	i++;
312	}
313	return result;
314	}
315
316	/**
317	* print the dynmaic programming matrix
318	*/
319	public void printDPMatrix() {
320	System.out.print(" ");
321	for (int i = 1; i <= length1; i++)
322	System.out.format("%5d", input1[i - 1]);
323	System.out.println();
324	for (int j = 0; j <= length2; j++) {
325	if (j > 0)
326	System.out.format("%5d ",input2[j - 1]);
327	else{
328	System.out.print(" ");
329	}
330	for (int i = 0; i <= length1 + 1; i++) {
331	if((i<length1+1) \|\| (i==length1+1 && j==0)) {
332	System.out.format("%4.1f ",matrix[i][j].getScore());
333	}
334
335	}
336	System.out.println();
337	}
338	}
339
340
341	/* (non-Javadoc)
342	* @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignment()
343	*/
344	@Override
345	public ArrayList<NumberSequence> getAlignment() {
346	return alignment;
347	}
348
349	public void setAlignment(ArrayList<NumberSequence> alignment) {
350	this.alignment = alignment;
351	}
352
353	@Override
354	public void align(NumberSequence input1, NumberSequence input2, SubstitutionMatrix submat,
355	float threshold) {
356
357	alignment = new ArrayList<NumberSequence>();
358	alignment.add(input1);
359	alignment.add(input2);
360
361	this.input1=input1.getSequence();
362	this.input2=input2.getSequence();
363
364	length1 = input1.size();
365	length2 = input2.size();
366	this.submat = submat;
367
368	//System.out.println("Starting SmithWaterman algorithm with a "
369	// + submat.getClass() + " Substitution Matrix: " + submat.getClass().getCanonicalName());
370	this.scoreThreshold = threshold;
371
372	matrix = new MatrixEntry[length1+2][length2+1];
373
374
375	for (int i = 0; i <= length1+1; i++) {
376	for(int j = 0; j<= length2; j++) {
377	matrix[i][j] = new MatrixEntry();
378	}
379	}
380
381	//Console.traceln(Level.INFO,"Generating DP Matrix");
382	buildMatrix();
383	//Console.traceln(Level.INFO,"Doing traceback");
384	traceback();
385	}
386
387	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: