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

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

Refactoring and code cleanup
File size: 8.7 KB
RevLine 
[1558]1package de.ugoe.cs.autoquest.tasktrees.alignment.algorithms;
2
3import java.util.ArrayList;
[1589]4import java.util.Iterator;
5import java.util.LinkedList;
[1575]6import java.util.logging.Level;
[1558]7
[1572]8import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.SubstitutionMatrix;
[1578]9import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.Constants;
[1575]10import de.ugoe.cs.util.console.Console;
[1558]11
[1586]12public class SmithWatermanRepeated implements AlignmentAlgorithm {
[1558]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
[1572]35
[1585]36        private ArrayList<NumberSequence> alignment;
[1558]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];
[1572]57                alignment = new ArrayList<NumberSequence>();
[1558]58               
[1559]59                for (int i = 0; i <= length1+1; i++) {
[1558]60                        for(int j = 0; j< length2; j++) {
61                                matrix[i][j] = new MatrixEntry();
62                        }
63                }
64       
65
66                buildMatrix();
[1574]67                traceback();
[1558]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         */
[1568]80        private double similarity(int i, int j) {
[1578]81                return submat.getScore(input1[i - 1], input2[j - 1]);
[1558]82        }
83
84        /**
[1559]85         * Build the score matrix using dynamic programming.
[1558]86         */
87        private void buildMatrix() {
88                if (submat.getGapPenalty() >= 0) {
89                        throw new Error("Indel score must be negative");
90                }
91               
[1559]92                // it's a gap
[1558]93                matrix[0][0].setScore(0);
94                matrix[0][0].setPrevious(null); // starting point
95
96                // the first column
97                for (int j = 1; j < length2; j++) {
98                        matrix[0][j].setScore(0);
[1575]99                        //We don't need to go back to [0][0] if we reached matrix[0][x], so just end here
100                        //matrix[0][j].setPrevious(matrix[0][j-1]);
101                        matrix[0][j].setPrevious(null);
[1558]102                }
103               
104               
105               
[1587]106                for (int i = 1; i < length1 + 2; i++) {
[1558]107               
108                        // Formula for first row:
109                        // F(i,0) = max { F(i-1,0), F(i-1,j)-T  j=1,...,m
110                       
[1568]111                        double firstRowLeftScore = matrix[i-1][0].getScore();
[1559]112                        //for sequences of length 1
[1568]113                        double tempMax;
[1559]114                        int maxRowIndex;
115                        if(length2 == 1) {
116                                tempMax = matrix[i-1][0].getScore();
117                                maxRowIndex = 0;
118                        } else {
119                                tempMax = matrix[i-1][1].getScore();
120                                maxRowIndex = 1;
121                                //position of the maximal score of the previous row
122                               
123                                for(int j = 2; j < length2;j++) {
124                                        if(matrix[i-1][j].getScore() > tempMax) {
125                                                tempMax = matrix[i-1][j].getScore();
126                                                maxRowIndex = j;
127                                        }
[1558]128                                }
[1559]129
[1558]130                        }
[1559]131                                       
[1558]132                        tempMax -= scoreThreshold;
133                        matrix[i][0].setScore(Math.max(firstRowLeftScore, tempMax));
[1559]134                        if(tempMax ==matrix[i][0].getScore()){
135                                matrix[i][0].setPrevious(matrix[i-1][maxRowIndex]);
136                        }
[1558]137                       
[1559]138                        if(firstRowLeftScore == matrix[i][0].getScore()) {
139                                matrix[i][0].setPrevious(matrix[i-1][0]);
140                        }
[1558]141                       
[1559]142                        //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
143                        if(i<length1+1)
144                        {
145                                matrix[i][0].setXvalue(input1[i-1]);
[1578]146                                matrix[i][0].setYvalue(Constants.UNMATCHED_SYMBOL);
[1559]147                        }
148                        else {
149                        //End after we calculated final score
150                                return;
151                        }
152                       
153                       
[1558]154                        for (int j = 1; j < length2; j++) {
[1568]155                                double diagScore = matrix[i - 1][j - 1].getScore() + similarity(i, j);
156                                double upScore = matrix[i][j - 1].getScore() + submat.getGapPenalty();
157                                double leftScore = matrix[i - 1][j].getScore() + submat.getGapPenalty();
[1558]158
159                                matrix[i][j].setScore(Math.max(diagScore,Math.max(upScore, Math.max(leftScore,matrix[i][0].getScore()))));
160
161                                // find the directions that give the maximum scores.
[1578]162                                // TODO: Multiple directions are ignored, we choose the first maximum score
[1559]163                                //True if we had a match
[1558]164                                if (diagScore == matrix[i][j].getScore()) {
165                                        matrix[i][j].setPrevious(matrix[i-1][j-1]);
[1559]166                                        matrix[i][j].setXvalue(input1[i-1]);
167                                        matrix[i][j].setYvalue(input2[j-1]);
[1558]168                                }
[1559]169                                //true if we took an event from sequence x and not from y
[1558]170                                if (leftScore == matrix[i][j].getScore()) {
[1559]171                                        matrix[i][j].setXvalue(input1[i-1]);
[1578]172                                        matrix[i][j].setYvalue(Constants.GAP_SYMBOL);
[1558]173                                        matrix[i][j].setPrevious(matrix[i-1][j]);
174                                }
[1559]175                                //true if we took an event from sequence y and not from x
[1558]176                                if (upScore == matrix[i][j].getScore()) {
[1578]177                                        matrix[i][j].setXvalue(Constants.GAP_SYMBOL);
[1559]178                                        matrix[i][j].setYvalue(input2[j-1]);
[1558]179                                        matrix[i][j].setPrevious(matrix[i][j-1]);
180                                }
[1559]181                                //true if we ended a matching region
[1558]182                                if (matrix[i][0].getScore() == matrix[i][j].getScore()) {
[1559]183                                        matrix[i][j].setPrevious(matrix[i][0]);
184                                        matrix[i][j].setXvalue(input1[i-1]);
[1578]185                                        matrix[i][j].setYvalue(Constants.UNMATCHED_SYMBOL);
[1558]186                                }
187                        }
[1559]188                       
189                        //Set the complete score cell
190                       
[1558]191                }
192        }
193
194        /**
195         * Get the maximum value in the score matrix.
196         */
197        public double getMaxScore() {
198                double maxScore = 0;
199
200                // skip the first row and column
201                for (int i = 1; i <= length1; i++) {
202                        for (int j = 1; j < length2; j++) {
203                                if (matrix[i][j].getScore() > maxScore) {
204                                        maxScore = matrix[i][j].getScore();
205                                }
206                        }
207                }
208
209                return maxScore;
210        }
211
[1586]212        /* (non-Javadoc)
213         * @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignmentScore()
[1558]214         */
[1586]215        @Override
[1568]216        public double getAlignmentScore() {
[1559]217                return matrix[length1+1][0].getScore();
[1558]218        }
219
[1572]220        public void traceback() {
[1559]221                MatrixEntry tmp = matrix[length1+1][0];
[1589]222                LinkedList<Integer> aligned1 = new LinkedList<Integer>();
223                LinkedList<Integer> aligned2 = new LinkedList<Integer>();
224                do {
[1572]225                       
[1589]226                        aligned1.add(new Integer(tmp.getXvalue()));
227                        aligned2.add(new Integer(tmp.getYvalue()));
228
[1572]229                        tmp = tmp.getPrevious();
[1589]230
231                } while (tmp != null);
232               
233                // reverse order of the alignment
234                int reversed1[] = new int[aligned1.size()];
235                int reversed2[] = new int[aligned2.size()];
236
237                int count = 0;
238                for (Iterator<Integer> it = aligned1.descendingIterator(); it.hasNext();) {
[1572]239                        count++;
[1589]240                        reversed1[reversed1.length - count] = it.next();
241                       
[1572]242                }
[1589]243                count = 0;
244                for (Iterator<Integer> it = aligned2.descendingIterator(); it.hasNext();) {
245                        count++;
246                        reversed2[reversed2.length - count] = it.next();
247                }
248
[1572]249                NumberSequence ns1 = new NumberSequence(reversed1.length);
250                NumberSequence ns2 = new NumberSequence(reversed2.length);
251                ns1.setSequence(reversed1);
252                ns2.setSequence(reversed2);
[1589]253
[1572]254                alignment.add(ns1);
255                alignment.add(ns2);
256        }
257       
[1589]258
259       
[1572]260        public void printAlignment() {
261                MatrixEntry tmp = matrix[length1+1][0];
[1559]262                String aligned1 = "";
263                String aligned2 = "";
264                int count = 0;
265                do
266                {
267                        String append1="";
268                        String append2="";
269                                       
[1578]270                        if(tmp.getXvalue() == Constants.GAP_SYMBOL) {
[1559]271                                append1 = "  ___";
272                        }
[1578]273                        else if(tmp.getXvalue() == Constants.UNMATCHED_SYMBOL) {
[1559]274                                append1 = "  ...";
275                        }
276                        else {
277                                append1 = String.format("%5d", tmp.getXvalue());
278                        }
[1558]279
[1578]280                        if(tmp.getYvalue() == Constants.GAP_SYMBOL) {
[1559]281                                append2 = "  ___";
282                        }
[1578]283                        else if(tmp.getYvalue() == Constants.UNMATCHED_SYMBOL) {
[1559]284                                append2 = "  ...";
285                        }
286                        else {
287                                append2 = String.format("%5d", tmp.getYvalue());
288                        }
289                        if(count != 0)
290                        {
291                                aligned1 = append1 + aligned1;
292                                aligned2 = append2 + aligned2;
293                        }
294                       
295                        tmp = tmp.getPrevious();
296                        count++;
297                       
298                } while(tmp != null);
299                System.out.println(aligned1);
300                System.out.println(aligned2);
301        }
[1558]302       
[1559]303
304       
[1558]305        /**
306         * print the dynmaic programming matrix
307         */
308        public void printDPMatrix() {
309                System.out.print("          ");
310                for (int i = 1; i <= length1; i++)
311                        System.out.format("%5d", input1[i - 1]);
312                System.out.println();
[1587]313                for (int j = 0; j <= length2; j++) {
[1558]314                        if (j > 0)
315                                System.out.format("%5d ",input2[j - 1]);
316                        else{
317                                System.out.print("      ");
318                        }
[1559]319                        for (int i = 0; i <= length1 + 1; i++) {
320                                if((i<length1+1) || (i==length1+1 && j==0))     {
321                                        System.out.format("%4.1f ",matrix[i][j].getScore());
322                                }
323                               
324                        }                       
[1558]325                        System.out.println();
326                }
327        }
328
329
[1586]330        /* (non-Javadoc)
331         * @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignment()
332         */
333        @Override
[1585]334        public ArrayList<NumberSequence> getAlignment() {
[1572]335                return alignment;
336        }
337
[1585]338        public void setAlignment(ArrayList<NumberSequence> alignment) {
[1572]339                this.alignment = alignment;
340        }
341
[1558]342}
Note: See TracBrowser for help on using the repository browser.