--- /dev/null
+package fr.orsay.lri.varna.models.treealign;
+import java.util.ArrayList;
+
+
+
+/**
+ *
+ *
+ * @author Raphael Champeimont
+ */
+public class ExampleDistance3 implements TreeAlignLabelDistanceSymmetric<RNANodeValue2> {
+ public double f(RNANodeValue2 v1, RNANodeValue2 v2) {
+ if (v1 == null) {
+ if (v2 == null) {
+ return 0;
+ } else if (!v2.isSingleNode()) { // v2 is a list of bases
+ // We insert all bases, with a cost of 1 for each base.
+ return v2.getNodes().size();
+ } else { // v2 is a single node
+ return 2;
+ }
+ } else if (!v1.isSingleNode()) { // v1 is a list of bases
+ if (v2 == null) {
+ return v1.getNodes().size();
+ } else if (!v2.isSingleNode()) { // v2 is a list of bases
+ // We compute the sequence distance
+ return alignSequenceNodes(v1, v2).getDistance();
+ } else { // v2 is a single node
+ return 2 + v1.getNodes().size();
+ }
+ } else { // v1 is a single node
+ // all the same as when v1 == null
+ if (v2 == null) {
+ return 2;
+ } else if (!v2.isSingleNode()) { // v2 is a list of bases
+ return 2 + v2.getNodes().size();
+ } else { // v2 is a single node
+ String l1 = v1.getNode().getLeftNucleotide();
+ String r1 = v1.getNode().getRightNucleotide();
+ String l2 = v2.getNode().getLeftNucleotide();
+ String r2 = v2.getNode().getRightNucleotide();
+ // We have cost(subst((x,y) to (x',y'))) = 1
+ // when x != x' and y != y'.
+ // It means it is less than substituting 2 non-paired bases
+ return (!l1.equals(l2) ? 0.5 : 0)
+ + (!r1.equals(r2) ? 0.5 : 0);
+ }
+ }
+ }
+
+
+ public class SequenceAlignResult {
+ private double distance;
+ private int[][] alignment;
+
+ public double getDistance() {
+ return distance;
+ }
+ public void setDistance(double distance) {
+ this.distance = distance;
+ }
+
+ /** The result array is a matrix of height 2
+ * and width at most length(sequence A) + length(sequence B).
+ * with result[0] is the alignment for A
+ * and result[1] the alignment for B.
+ * The alignment consists int the indexes of the original
+ * bases positions, with -1 when there is no match.
+ */
+ public int[][] getAlignment() {
+ return alignment;
+ }
+ public void setAlignment(int[][] alignment) {
+ this.alignment = alignment;
+ }
+
+ }
+
+ /**
+ * Align two sequences contained in nodes.
+ * Both nodes have to be non-single nodes, otherwise an
+ * RNANodeValue2WrongTypeException exception will be thrown.
+ */
+ public SequenceAlignResult alignSequenceNodes(RNANodeValue2 v1, RNANodeValue2 v2) {
+ char[] A = v1.computeSequence();
+ char[] B = v2.computeSequence();
+ return alignSequences(A, B);
+ }
+
+ /**
+ * Align sequences using the Needleman-Wunsch algorithm.
+ * Time: O(A.length * B.length)
+ * Space: O(A.length * B.length)
+ * Space used by the returned object: O(A.length + B.length)
+ */
+ public SequenceAlignResult alignSequences(char[] A, char[] B) {
+ SequenceAlignResult result = new SequenceAlignResult();
+
+ final int la = A.length;
+ final int lb = B.length;
+ double[][] F = new double[la+1][lb+1];
+ int[][] decisions = new int[la+1][lb+1];
+ final double d = 1; // insertion/deletion cost
+ final double substCost = 1; // substitution cost
+ for (int i=0; i<=la; i++)
+ F[i][0] = d*i;
+ for (int j=0; j<=lb; j++)
+ F[0][j] = d*j;
+ for (int i=1; i<=la; i++)
+ for (int j=1; j<=lb; j++)
+ {
+ double min;
+ int decision;
+ double match = F[i-1][j-1] + (A[i-1] == B[j-1] ? 0 : substCost);
+ double delete = F[i-1][j] + d;
+ if (match < delete) {
+ decision = 1;
+ min = match;
+ } else {
+ decision = 2;
+ min = delete;
+ }
+ double insert = F[i][j-1] + d;
+ if (insert < min) {
+ decision = 3;
+ min = insert;
+ }
+ F[i][j] = min;
+ decisions[i][j] = decision;
+ }
+
+ result.setDistance(F[la][lb]);
+
+ int[][] alignment = computeAlignment(F, decisions, A, B);
+ result.setAlignment(alignment);
+
+ return result;
+ }
+
+ private int[][] computeAlignment(double[][] F, int[][] decisions, char[] A, char[] B) {
+ // At worst the alignment will be of length (A.length + B.length)
+ ArrayList<Integer> AlignmentA = new ArrayList<Integer>(A.length + B.length);
+ ArrayList<Integer> AlignmentB = new ArrayList<Integer>(A.length + B.length);
+ int i = A.length;
+ int j = B.length;
+ while (i > 0 && j > 0)
+ {
+ int decision = decisions[i][j];
+ switch (decision) {
+ case 1:
+ AlignmentA.add(i-1);
+ AlignmentB.add(j-1);
+ i = i - 1;
+ j = j - 1;
+ break;
+ case 2:
+ AlignmentA.add(i-1);
+ AlignmentB.add(-1);
+ i = i - 1;
+ break;
+ case 3:
+ AlignmentA.add(-1);
+ AlignmentB.add(j-1);
+ j = j - 1;
+ break;
+ default:
+ throw (new Error("Bug in ExampleDistance3: decision = " + decision));
+ }
+ }
+ while (i > 0)
+ {
+ AlignmentA.add(i-1);
+ AlignmentB.add(-1);
+ i = i - 1;
+ }
+ while (j > 0)
+ {
+ AlignmentA.add(-1);
+ AlignmentB.add(j-1);
+ j = j - 1;
+ }
+
+ // Convert the ArrayLists to the right format:
+ // We need to reverse the list and change the numbering of
+ int l = AlignmentA.size();
+ int[][] result = new int[2][l];
+ for (i=0; i<l; i++) {
+ result[0][i] = AlignmentA.get(l-1-i);
+ result[1][i] = AlignmentB.get(l-1-i);
+ }
+ return result;
+
+ }
+}
git://source.jalview.org / jalview.git / blobdiff