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

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

Removed parameters from alignmentalgorihm factory constructor and changed interface by adding a new align() method, which now gets all the data via parameter
File size: 7.3 KB
Line 
1package de.ugoe.cs.autoquest.tasktrees.alignment.algorithms;
2
3import java.util.ArrayList;
4import java.util.logging.Level;
5
6import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.SubstitutionMatrix;
7import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.Constants;
8import de.ugoe.cs.util.console.Console;
9
10public class SmithWaterman implements AlignmentAlgorithm {
11
12        /**
13         * The first input
14         */
15        private int[] input1;
16
17        /**
18         * The second input String
19         */
20        private int[] input2;
21
22        /**
23         * The lengths of the input
24         */
25        private int length1, length2;
26
27        /**
28         * The score matrix. The true scores should be divided by the normalization
29         * factor.
30         */
31        private MatrixEntry[][] matrix;
32
33        private ArrayList<NumberSequence> alignment;
34
35        /**
36         * Substitution matrix to calculate scores
37         */
38        private SubstitutionMatrix submat;
39
40
41        /**
42         * Compute the similarity score of substitution The position of the first
43         * character is 1. A position of 0 represents a gap.
44         *
45         * @param i
46         *            Position of the character in str1
47         * @param j
48         *            Position of the character in str2
49         * @return Cost of substitution of the character in str1 by the one in str2
50         */
51        private double similarity(int i, int j) {
52                return submat.getScore(input1[i - 1], input2[j - 1]);
53        }
54
55        /**
56         * Build the score matrix using dynamic programming.
57         */
58        private void buildMatrix() {
59                if (submat.getGapPenalty() >= 0) {
60                        throw new Error("Indel score must be negative");
61                }
62
63                // it's a gap
64                matrix[0][0].setScore(0);
65                matrix[0][0].setPrevious(null); // starting point
66
67                // the first column
68                for (int j = 1; j < length2; j++) {
69                        matrix[0][j].setScore(0);
70                        matrix[0][j].setPrevious(matrix[0][j - 1]);
71                        matrix[0][j].setYvalue(input2[j]);
72                        matrix[0][j].setXvalue(Constants.GAP_SYMBOL);
73                }
74                // the first row
75
76                for (int j = 1; j < length1; j++) {
77                        matrix[j][0].setScore(0);
78                        matrix[j][0].setPrevious(matrix[j - 1][0]);
79                        matrix[j][0].setXvalue(input1[j]);
80                        matrix[j][0].setYvalue(Constants.GAP_SYMBOL);
81                }
82
83                for (int i = 1; i < length1; i++) {
84
85                        for (int j = 1; j < length2; j++) {
86                                double diagScore = matrix[i - 1][j - 1].getScore()
87                                                + similarity(i, j);
88                                double upScore = matrix[i][j - 1].getScore()
89                                                + submat.getGapPenalty();
90                                double leftScore = matrix[i - 1][j].getScore()
91                                                + submat.getGapPenalty();
92
93                                matrix[i][j].setScore(Math.max(diagScore,
94                                                Math.max(upScore, Math.max(leftScore, 0))));
95
96                                // find the directions that give the maximum scores.
97                                // TODO: Multiple directions are ignored, we choose the first
98                                // maximum score
99                                // True if we had a match
100                                if (diagScore == matrix[i][j].getScore()) {
101                                        matrix[i][j].setPrevious(matrix[i - 1][j - 1]);
102                                        matrix[i][j].setXvalue(input1[i - 1]);
103                                        matrix[i][j].setYvalue(input2[j - 1]);
104                                }
105                                // true if we took an event from sequence x and not from y
106                                if (leftScore == matrix[i][j].getScore()) {
107                                        matrix[i][j].setXvalue(input1[i - 1]);
108                                        matrix[i][j].setYvalue(Constants.GAP_SYMBOL);
109                                        matrix[i][j].setPrevious(matrix[i - 1][j]);
110                                }
111                                // true if we took an event from sequence y and not from x
112                                if (upScore == matrix[i][j].getScore()) {
113                                        matrix[i][j].setXvalue(Constants.GAP_SYMBOL);
114                                        matrix[i][j].setYvalue(input2[j - 1]);
115                                        matrix[i][j].setPrevious(matrix[i][j - 1]);
116                                }
117                                if (0 == matrix[i][j].getScore()) {
118                                        matrix[i][j].setPrevious(matrix[0][0]);
119                                        matrix[i][j].setXvalue(-2);
120                                        matrix[i][j].setYvalue(-2);
121                                }
122                        }
123
124                        // Set the complete score cell
125
126                }
127        }
128
129        /**
130         * Get the maximum value in the score matrix.
131         */
132        public double getMaxScore() {
133                double maxScore = 0;
134
135                // skip the first row and column
136                for (int i = 1; i <= length1; i++) {
137                        for (int j = 1; j < length2; j++) {
138                                if (matrix[i][j].getScore() > maxScore) {
139                                        maxScore = matrix[i][j].getScore();
140                                }
141                        }
142                }
143
144                return maxScore;
145        }
146
147        /*
148         * (non-Javadoc)
149         *
150         * @see
151         * de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm
152         * #getAlignmentScore()
153         */
154        @Override
155        public double getAlignmentScore() {
156                return getMaxScore();
157        }
158
159        public void traceback() {
160                MatrixEntry tmp = matrix[length1][length2];
161                int aligned1[] = new int[length1 + length2 + 2];
162                int aligned2[] = new int[length1 + length2 + 2];
163                int count = 0;
164                do {
165                        if (length1 + length2 + 2 == count) {
166                                Console.traceln(Level.WARNING,
167                                                "Traceback longer than both sequences summed up!");
168                                break;
169                        }
170                        aligned1[count] = tmp.getXvalue();
171                        aligned2[count] = tmp.getYvalue();
172
173                        tmp = tmp.getPrevious();
174                        count++;
175
176                } while (tmp != null);
177                count--;
178                // reverse order of the alignment
179                int reversed1[] = new int[count];
180                int reversed2[] = new int[count];
181
182                for (int i = count - 1; i > 0; i--) {
183                        reversed1[reversed1.length - i] = aligned1[i];
184                        reversed2[reversed2.length - i] = aligned2[i];
185                }
186
187                NumberSequence ns1 = new NumberSequence(reversed1.length);
188                NumberSequence ns2 = new NumberSequence(reversed2.length);
189                ns1.setSequence(reversed1);
190                ns2.setSequence(reversed2);
191
192                alignment.add(ns1);
193                alignment.add(ns2);
194        }
195
196        public void printAlignment() {
197                MatrixEntry tmp = matrix[length1 + 1][0];
198                String aligned1 = "";
199                String aligned2 = "";
200                int count = 0;
201                do {
202                        String append1 = "";
203                        String append2 = "";
204
205                        if (tmp.getXvalue() == Constants.GAP_SYMBOL) {
206                                append1 = "  ___";
207                        } else {
208                                append1 = String.format("%5d", tmp.getXvalue());
209                        }
210
211                        if (tmp.getYvalue() == Constants.GAP_SYMBOL) {
212                                append2 = "  ___";
213                        } else {
214                                append2 = String.format("%5d", tmp.getYvalue());
215                        }
216                        if (count != 0) {
217                                aligned1 = append1 + aligned1;
218                                aligned2 = append2 + aligned2;
219                        }
220
221                        tmp = tmp.getPrevious();
222                        count++;
223
224                } while (tmp != null);
225                System.out.println(aligned1);
226                System.out.println(aligned2);
227        }
228
229        /**
230         * print the dynmaic programming matrix
231         */
232        public void printDPMatrix() {
233                System.out.print("          ");
234                for (int i = 1; i <= length1; i++)
235                        System.out.format("%5d", input1[i - 1]);
236                System.out.println();
237                for (int j = 0; j <= length2; j++) {
238                        if (j > 0)
239                                System.out.format("%5d ", input2[j - 1]);
240                        else {
241                                System.out.print("      ");
242                        }
243                        for (int i = 0; i <= length1; i++) {
244                                        System.out.format("%4.1f ", matrix[i][j].getScore());
245                        }
246                        System.out.println();
247                }
248        }
249
250        /*
251         * (non-Javadoc)
252         *
253         * @see
254         * de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm
255         * #getAlignment()
256         */
257        @Override
258        public ArrayList<NumberSequence> getAlignment() {
259                return alignment;
260        }
261
262        public void setAlignment(ArrayList<NumberSequence> alignment) {
263                this.alignment = alignment;
264        }
265
266        @Override
267        public ArrayList<ArrayList<NumberSequence>> getMatches() {
268                // TODO Auto-generated method stub
269                return null;
270        }
271
272        @Override
273        public void align(int[] input1, int[] input2, SubstitutionMatrix submat,
274                        float threshold) {
275                this.input1 = input1;
276                this.input2 = input2;
277                length1 = input1.length;
278                length2 = input2.length;
279                this.submat = submat;
280
281                // System.out.println("Starting SmithWaterman algorithm with a "
282                // + submat.getClass() + " Substitution Matrix: " +
283                // submat.getClass().getCanonicalName());
284
285                matrix = new MatrixEntry[length1 + 1][length2 + 1];
286                alignment = new ArrayList<NumberSequence>();
287
288                for (int i = 0; i < length1+1; i++) {
289                        for (int j = 0; j < length2+1; j++) {
290                                matrix[i][j] = new MatrixEntry();
291                        }
292                }
293
294                buildMatrix();
295                traceback();
296               
297        }
298
299}
Note: See TracBrowser for help on using the repository browser.