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

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

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