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

package de.ugoe.cs.autoquest.tasktrees.alignment.algorithms;

import java.util.ArrayList;
import java.util.logging.Level;

import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.SubstitutionMatrix;
import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.Constants;
import de.ugoe.cs.util.console.Console;

public class SmithWaterman implements AlignmentAlgorithm {

        /**
         * The first input
         */
        private int[] input1;

        /**
         * The second input String
         */
        private int[] input2;

        /**
         * The lengths of the input
         */
        private int length1, length2;

        /**
         * The score matrix. The true scores should be divided by the normalization
         * factor.
         */
        private MatrixEntry[][] matrix;

        private ArrayList<NumberSequence> alignment;

        /**
         * Substitution matrix to calculate scores
         */
        private SubstitutionMatrix submat;

        @Override
        public void align(NumberSequence input1, NumberSequence input2,
                        SubstitutionMatrix submat, float threshold) {
                this.input1 = input1.getSequence();
                this.input2 = input2.getSequence();
                length1 = input1.size();
                length2 = input2.size();
                this.submat = submat;

                // System.out.println("Starting SmithWaterman algorithm with a "
                // + submat.getClass() + " Substitution Matrix: " +
                // submat.getClass().getCanonicalName());

                matrix = new MatrixEntry[length1 + 1][length2 + 1];
                alignment = new ArrayList<NumberSequence>();

                for (int i = 0; i < (length1 + 1); i++) {
                        for (int j = 0; j < (length2 + 1); j++) {
                                matrix[i][j] = new MatrixEntry();
                        }
                }

                buildMatrix();
                traceback();

        }

        /**
         * Build the score matrix using dynamic programming.
         */
        private void buildMatrix() {
                if (submat.getGapPenalty() >= 0) {
                        throw new Error("Indel score must be negative");
                }

                // it's a gap
                matrix[0][0].setScore(0);
                matrix[0][0].setPrevious(null); // starting point

                // the first column
                for (int j = 1; j < length2; j++) {
                        matrix[0][j].setScore(0);
                        matrix[0][j].setPrevious(matrix[0][j - 1]);
                        matrix[0][j].setYvalue(input2[j]);
                        matrix[0][j].setXvalue(Constants.GAP_SYMBOL);
                }
                // the first row

                for (int j = 1; j < length1; j++) {
                        matrix[j][0].setScore(0);
                        matrix[j][0].setPrevious(matrix[j - 1][0]);
                        matrix[j][0].setXvalue(input1[j]);
                        matrix[j][0].setYvalue(Constants.GAP_SYMBOL);
                }

                for (int i = 1; i < length1; i++) {

                        for (int j = 1; j < length2; j++) {
                                final double diagScore = matrix[i - 1][j - 1].getScore()
                                                + similarity(i, j);
                                final double upScore = matrix[i][j - 1].getScore()
                                                + submat.getGapPenalty();
                                final double leftScore = matrix[i - 1][j].getScore()
                                                + submat.getGapPenalty();

                                matrix[i][j].setScore(Math.max(diagScore,
                                                Math.max(upScore, Math.max(leftScore, 0))));

                                // find the directions that give the maximum scores.
                                // TODO: Multiple directions are ignored, we choose the first
                                // maximum score
                                // True if we had a match
                                if (diagScore == matrix[i][j].getScore()) {
                                        matrix[i][j].setPrevious(matrix[i - 1][j - 1]);
                                        matrix[i][j].setXvalue(input1[i - 1]);
                                        matrix[i][j].setYvalue(input2[j - 1]);
                                }
                                // true if we took an event from sequence x and not from y
                                if (leftScore == matrix[i][j].getScore()) {
                                        matrix[i][j].setXvalue(input1[i - 1]);
                                        matrix[i][j].setYvalue(Constants.GAP_SYMBOL);
                                        matrix[i][j].setPrevious(matrix[i - 1][j]);
                                }
                                // true if we took an event from sequence y and not from x
                                if (upScore == matrix[i][j].getScore()) {
                                        matrix[i][j].setXvalue(Constants.GAP_SYMBOL);
                                        matrix[i][j].setYvalue(input2[j - 1]);
                                        matrix[i][j].setPrevious(matrix[i][j - 1]);
                                }
                                if (0 == matrix[i][j].getScore()) {
                                        matrix[i][j].setPrevious(matrix[0][0]);
                                        matrix[i][j].setXvalue(-2);
                                        matrix[i][j].setYvalue(-2);
                                }
                        }

                        // Set the complete score cell

                }
        }

        /*
         * (non-Javadoc)
         * 
         * @see
         * de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm
         * #getAlignment()
         */
        @Override
        public ArrayList<NumberSequence> getAlignment() {
                return alignment;
        }

        /*
         * (non-Javadoc)
         * 
         * @see
         * de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm
         * #getAlignmentScore()
         */
        @Override
        public double getAlignmentScore() {
                return getMaxScore();
        }

        @Override
        public ArrayList<Match> getMatches() {
                // TODO Auto-generated method stub
                return null;
        }

        /**
         * Get the maximum value in the score matrix.
         */
        @Override
        public double getMaxScore() {
                double maxScore = 0;

                // skip the first row and column
                for (int i = 1; i <= length1; i++) {
                        for (int j = 1; j < length2; j++) {
                                if (matrix[i][j].getScore() > maxScore) {
                                        maxScore = matrix[i][j].getScore();
                                }
                        }
                }

                return maxScore;
        }

        @Override
        public void printAlignment() {
                MatrixEntry tmp = matrix[length1 + 1][0];
                String aligned1 = "";
                String aligned2 = "";
                int count = 0;
                do {
                        String append1 = "";
                        String append2 = "";

                        if (tmp.getXvalue() == Constants.GAP_SYMBOL) {
                                append1 = "  ___";
                        } else {
                                append1 = String.format("%5d", tmp.getXvalue());
                        }

                        if (tmp.getYvalue() == Constants.GAP_SYMBOL) {
                                append2 = "  ___";
                        } else {
                                append2 = String.format("%5d", tmp.getYvalue());
                        }
                        if (count != 0) {
                                aligned1 = append1 + aligned1;
                                aligned2 = append2 + aligned2;
                        }

                        tmp = tmp.getPrevious();
                        count++;

                } while (tmp != null);
                System.out.println(aligned1);
                System.out.println(aligned2);
        }

        /**
         * print the dynmaic programming matrix
         */
        @Override
        public void printDPMatrix() {
                System.out.print("          ");
                for (int i = 1; i <= length1; i++) {
                        System.out.format("%5d", input1[i - 1]);
                }
                System.out.println();
                for (int j = 0; j <= length2; j++) {
                        if (j > 0) {
                                System.out.format("%5d ", input2[j - 1]);
                        } else {
                                System.out.print("      ");
                        }
                        for (int i = 0; i <= length1; i++) {
                                System.out.format("%4.1f ", matrix[i][j].getScore());
                        }
                        System.out.println();
                }
        }

        public void setAlignment(ArrayList<NumberSequence> alignment) {
                this.alignment = alignment;
        }

        /**
         * Compute the similarity score of substitution The position of the first
         * character is 1. A position of 0 represents a gap.
         * 
         * @param i
         *            Position of the character in str1
         * @param j
         *            Position of the character in str2
         * @return Cost of substitution of the character in str1 by the one in str2
         */
        private double similarity(int i, int j) {
                return submat.getScore(input1[i - 1], input2[j - 1]);
        }

        public void traceback() {
                MatrixEntry tmp = matrix[length1][length2];
                final int aligned1[] = new int[length1 + length2 + 2];
                final int aligned2[] = new int[length1 + length2 + 2];
                int count = 0;
                do {
                        if ((length1 + length2 + 2) == count) {
                                Console.traceln(Level.WARNING,
                                                "Traceback longer than both sequences summed up!");
                                break;
                        }
                        aligned1[count] = tmp.getXvalue();
                        aligned2[count] = tmp.getYvalue();

                        tmp = tmp.getPrevious();
                        count++;

                } while (tmp != null);
                count--;
                // reverse order of the alignment
                final int reversed1[] = new int[count];
                final int reversed2[] = new int[count];

                for (int i = count - 1; i > 0; i--) {
                        reversed1[reversed1.length - i] = aligned1[i];
                        reversed2[reversed2.length - i] = aligned2[i];
                }

                final NumberSequence ns1 = new NumberSequence(reversed1.length);
                final NumberSequence ns2 = new NumberSequence(reversed2.length);
                ns1.setSequence(reversed1);
                ns2.setSequence(reversed2);

                alignment.add(ns1);
                alignment.add(ns2);
        }

}