source: branches/ralph/src/main/java/de/ugoe/cs/autoquest/tasktrees/alignment/algorithms/SmithWatermanRepeated.java @ 1596

Last change on this file since 1596 was 1596, checked in by rkrimmel, 10 years ago

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