package bgu.bio.algorithms.graphs;
import gnu.trove.list.array.TIntArrayList;
import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.Random;
import org.junit.Assert;
import org.junit.Test;
import bgu.bio.adt.graphs.Tree;
import bgu.bio.adt.graphs.WeightedGraph;
import bgu.bio.adt.rna.NodeLabel;
import bgu.bio.adt.rna.RNASpecificTree;
import bgu.bio.adt.rna.costs.RnaSpecificCostFunction;
import bgu.bio.adt.rna.costs.RnaSpecificSmoothCost;
import bgu.bio.adt.rna.costs.RnaSpecificTypeRelatedTreePruning;
import bgu.bio.adt.tuples.IntPair;
import bgu.bio.adt.tuples.IntPairComparator;
import bgu.bio.algorithms.alignment.AffineGapGlobalSequenceAlignmentNoMatrix;
import bgu.bio.algorithms.alignment.SequenceAlignment;
import bgu.bio.algorithms.graphs.hsa.HSA;
import bgu.bio.algorithms.graphs.hsa.matchers.BipartiteCavityMatcher;
import bgu.bio.algorithms.graphs.hsa.matchers.CostFlowNetwork;
import bgu.bio.algorithms.graphs.hsa.matchers.MatcherFactory;
import bgu.bio.algorithms.graphs.hsa.matchers.OrderedBipartiteCavityMatcher;
import bgu.bio.algorithms.graphs.hsa.matchers.OrderedMatcherFactory;
import bgu.bio.algorithms.graphs.hsa.matchers.UnorderedBipartiteCavityMatcherMinCostFlow;
import bgu.bio.algorithms.graphs.hsa.matchers.UnorderedMatcherFactory;
import bgu.bio.util.AffineGapScoringMatrix;
import bgu.bio.util.MathOperations;
import bgu.bio.util.ScoringMatrix;
import bgu.bio.util.alphabet.RNABasePairAlphabet;
import bgu.bio.util.alphabet.RnaAlphabet;
public class TestHSA {
@Test
public void testRNAShuffleExample() {
RnaSpecificTypeRelatedTreePruning prune = new RnaSpecificTypeRelatedTreePruning(
false, true);
RnaSpecificSmoothCost smooth = new RnaSpecificSmoothCost();
SequenceAlignment aligner = new AffineGapGlobalSequenceAlignmentNoMatrix(
10, 10, RnaAlphabet.getInstance(), new AffineGapScoringMatrix(
"matrix" + File.separator
+ "interval-RIBOSUM85-60.matrix",
RnaAlphabet.getInstance()));
RnaSpecificCostFunction cost = new RnaSpecificCostFunction(
new ScoringMatrix("matrix" + File.separator
+ "bp-RIBOSUM85-60.matrix",
RNABasePairAlphabet.getInstance()), aligner);
final RNASpecificTree t = new RNASpecificTree("T");
t.buildFromViennaFormat(
"UUUAGCAGUUGGGAGCGCCUGAAGGGAGGUCCUNCGAUCGAAACCA".toCharArray(),
"((..((........))(((.....)))....(.......)))....".toCharArray());
final RNASpecificTree s = new RNASpecificTree("S");
s.buildFromViennaFormat(
"GAUAGAUGGUCUGGGUGUCGCCAGAUCGGGGUNCAAUUCCCCGUCGCCA"
.toCharArray(),
"((..(......)(((.....)))...(((((.......)))))))...."
.toCharArray());
Random rand = new Random(123l);
final RNASpecificTree sR = new RNASpecificTree("S-Random");
sR.buildFromViennaFormat(
"GAUAGAUGGUCUGGGUGUCGCCAGAUCGGGGUNCAAUUCCCCGUCGCCA"
.toCharArray(),
"((..(......)(((.....)))...(((((.......)))))))...."
.toCharArray());
sR.shuffle(rand);
prune.calculateCost(t);
prune.calculateCost(s);
prune.calculateCost(sR);
smooth.calculateCost(t);
smooth.calculateCost(s);
smooth.calculateCost(sR);
double ans1 = 0, ans2 = 0;
// reuse
HSA hsaEngine2 = new HSA(1, 1, cost, new UnorderedMatcherFactory(),
true);
hsaEngine2.computeHSA(t, s);
ans1 = hsaEngine2.computeHSA(t, sR);
// no reuse
hsaEngine2 = new HSA(1, 1, cost, new UnorderedMatcherFactory(), true);
hsaEngine2.computeHSA(t, s);
hsaEngine2.setFactory(new UnorderedMatcherFactory());
ans2 = hsaEngine2.computeHSA(t, sR);
Assert.assertEquals("Got a problem with the matchers", true,
MathOperations.equals(ans1, ans2));
}
@Test
public void testUnrootedExample() {
RnaSpecificTypeRelatedTreePruning pruneRooted = new RnaSpecificTypeRelatedTreePruning(
false, true);
RnaSpecificTypeRelatedTreePruning pruneUnrooted = new RnaSpecificTypeRelatedTreePruning(
false, false);
RnaSpecificSmoothCost smooth = new RnaSpecificSmoothCost();
SequenceAlignment aligner = new AffineGapGlobalSequenceAlignmentNoMatrix(
10, 10, RnaAlphabet.getInstance(), new AffineGapScoringMatrix(
"matrix" + File.separator
+ "interval-RIBOSUM85-60.matrix",
RnaAlphabet.getInstance()));
RnaSpecificCostFunction cost = new RnaSpecificCostFunction(
new ScoringMatrix("matrix" + File.separator
+ "bp-RIBOSUM85-60.matrix",
RNABasePairAlphabet.getInstance()), aligner);
final RNASpecificTree t = new RNASpecificTree("T");
t.buildFromViennaFormat(
"GGGGGGGGuuuAAAAAAAAAAcccUUUUUUUUUUuuuGGGGGGcccUUUUUUuuuCCCCCCCC"
.toCharArray(),
"((((((((...((((((((((...))))))))))...((((((...))))))...))))))))"
.toCharArray());
final RNASpecificTree s = new RNASpecificTree("S");
s.buildFromViennaFormat(
"UUUUUUUUUUuuuGGGGGGcccUUUUUUuuuCCCCCCCCcccGGGGGGGGuuuAAAAAAAAAA"
.toCharArray(),
"((((((((((...((((((...))))))...((((((((...))))))))...))))))))))"
.toCharArray());
smooth.calculateCost(t);
smooth.calculateCost(s);
double ans1 = 0, ans2 = 0;
// reuse
HSA hsaEngine = new HSA(1, 1, cost, new OrderedMatcherFactory(), true);
pruneRooted.calculateCost(t);
pruneRooted.calculateCost(s);
ans1 = hsaEngine.computeHSA(t, s);
pruneUnrooted.calculateCost(t);
pruneUnrooted.calculateCost(s);
TIntArrayList[] tb = new TIntArrayList[3];
ans2 = hsaEngine.computeHSA(t, s, tb);
for (int i = 0; i < tb[0].size(); i++) {
System.out.println(tb[0].getQuick(i) + " -> " + tb[1].getQuick(i));
}
Assert.assertEquals("Should be unrooted alignment", false,
MathOperations.equals(ans1, ans2));
}
@Test
public void testOrderedBipartiteMatching() {
// 0,1,2,3,4,5,6,7
// X: A,u,v,B,C,w,D,p
// Y: x,C,y,D,A,B
// Cavity alignment 7-5 score: -2 + 2*6 + 3 = 13
// A,u,v,B,_,C,w,D,_
// _,_,_,_,x,C,y,D,A
// Opt alignment score: 4*(-1) + 4*2 + 3 = 7
// A,u,v,B,_,C,w,D,p
// A,_,_,B,x,C,y,D,_
int sizeX = 8;
int sizeY = 6;
BipartiteCavityMatcher matcher = new OrderedBipartiteCavityMatcher(
sizeX, sizeY);
for (int y = 0; y < sizeY; ++y) {
matcher.setDelCostY(y, 2);
}
for (int x = 0; x < sizeX; ++x) {
matcher.setDelCostX(x, 2);
for (int y = 0; y < sizeY; ++y) {
matcher.setMatchCost(x, y, 5);
}
}
matcher.setMatchCost(0, 4, -1);
matcher.setMatchCost(3, 5, -1);
matcher.setMatchCost(4, 1, -1);
matcher.setMatchCost(6, 3, -1);
matcher.setMatchCost(5, 2, 3);
Assert.assertEquals("Cavity matching answer is wrong", 13,
matcher.minCostCavityMatching(7, 5), 0.001);
TIntArrayList[] matching = matcher.getCavityMatching(7, 5);
IntPairComparator comp = new IntPairComparator();
ArrayList<IntPair> ans = new ArrayList<IntPair>();
ArrayList<IntPair> exp = new ArrayList<IntPair>();
ans.add(new IntPair(6, 3));
ans.add(new IntPair(5, 2));
ans.add(new IntPair(4, 1));
for (int i = 0; i < matching[0].size(); ++i) {
exp.add(new IntPair(matching[0].get(i), matching[1].get(i)));
}
Collections.sort(ans, comp);
Collections.sort(exp, comp);
isSame(ans, exp);
Assert.assertEquals("Minimum matching answer is wrong", 7,
matcher.minCostMatching(), 0.001);
ans = new ArrayList<IntPair>();
exp = new ArrayList<IntPair>();
matching = matcher.getCurrMatching();
ans.add(new IntPair(3, 5));
ans.add(new IntPair(0, 4));
ans.add(new IntPair(6, 3));
ans.add(new IntPair(5, 2));
ans.add(new IntPair(4, 1));
for (int i = 0; i < matching[0].size(); ++i) {
exp.add(new IntPair(matching[0].get(i), matching[1].get(i)));
}
Collections.sort(ans, comp);
Collections.sort(exp, comp);
isSame(ans, exp);
}
@Test
public void testSelfSim() {
RnaSpecificTypeRelatedTreePruning prune = new RnaSpecificTypeRelatedTreePruning(
true, false);
RnaSpecificSmoothCost smooth = new RnaSpecificSmoothCost();
SequenceAlignment aligner = new AffineGapGlobalSequenceAlignmentNoMatrix(
10, 10, RnaAlphabet.getInstance(), new AffineGapScoringMatrix(
"matrix" + File.separator
+ "interval-RIBOSUM85-60.matrix",
RnaAlphabet.getInstance()));
RnaSpecificCostFunction cost = new RnaSpecificCostFunction(
new ScoringMatrix("matrix" + File.separator
+ "bp-RIBOSUM85-60.matrix",
RNABasePairAlphabet.getInstance()), aligner);
final RNASpecificTree t = new RNASpecificTree("T");
/*t.buildFromViennaFormat(
"UUUAGCAGUUGGGAGCGCCUGAAGGGAGGUCCUNCGAUCGAAACCA".toCharArray(),
"((..((........))(((.....)))....(.......)))....".toCharArray());
*/
t.buildFromViennaFormat(
"GCCUGAAGGGACGCCUGAAGGGA".toCharArray(),
"(((.....))).(((.....)))".toCharArray());
prune.calculateCost(t);
smooth.calculateCost(t);
try {
t.toDotFile("/home/milon/tmp/1.dot", true);
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
double ans1 = 0, ans2 = 0;
// reuse
HSA hsaEngine2 = new HSA(1, 1, cost, new OrderedMatcherFactory(),
true);
TIntArrayList[] trace = new TIntArrayList[3];
ans1 = hsaEngine2.computeHSA(t, t, trace);
Assert.assertEquals("trace should be the same", trace[0], trace[1]);
ans2 = hsaEngine2.computeHSA(t, trace);
System.out.println(ans1);
System.out.println(ans2);
Assert.assertEquals(
"Can't be the same score as regular comparison. should choose between non-trivial changes",
false, MathOperations.equals(ans1, ans2));
}
private void isSame(ArrayList<IntPair> ans, ArrayList<IntPair> exp) {
Assert.assertEquals("Wrong amount of matches", exp.size(), ans.size());
for (int i = 0; i < ans.size(); i++) {
Assert.assertEquals("Incorrect match at index 0", exp.get(i)
.getFirst(), ans.get(i).getFirst());
Assert.assertEquals("Incorrect match at index 1", exp.get(i)
.getSecond(), ans.get(i).getSecond());
}
}
private static void testHighDegHsa() {
int[][] tEdges = { { 1, 2, 3, 4, 5, 6, 7, 8, 9 }, { 0 }, { 0 }, { 0 },
{ 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, };
double[][] tWeights = { { 10, 10, 10, 10, 10, 10, 10, 10, 10 }, { 10 },
{ 20 }, { 30 }, { 40 }, { 50 }, { 60 }, { 70 }, { 80 }, { 90 }, };
double[] tSmooth = { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 };
NodeLabel[] tLabels = { new NodeLabel("0"), new NodeLabel("1"),
new NodeLabel("2"), new NodeLabel("3"), new NodeLabel("4"),
new NodeLabel("5"), new NodeLabel("6"), new NodeLabel("7"),
new NodeLabel("8"), new NodeLabel("9"), };
Tree t = new Tree(tEdges, tWeights, tSmooth, tLabels);
int[][] sEdges = { { 1, 2, 3, 4, 5 }, { 0 }, { 0 }, { 0 }, { 0 },
{ 0 }, };
double[][] sWeights = { { 10, 20, 30, 40, 50 }, { 10 }, { 20 }, { 30 },
{ 40 }, { 50 }, };
double[] sSmooth = { 2, 2, 2, 2, 2, 2 };
NodeLabel[] sLabels = { new NodeLabel("0"), new NodeLabel("1"),
new NodeLabel("2"), new NodeLabel("3"), new NodeLabel("4"),
new NodeLabel("5"), };
Tree s = new Tree(sEdges, sWeights, sSmooth, sLabels);
CostFunction w = new CostFunction() {
@Override
public double cost(NodeLabel l1, NodeLabel l2) {
double diff = Math.abs(l1.getLabelValue().get(0)
- l2.getLabelValue().get(0));
if (diff <= 3)
return diff - 1;
return Double.POSITIVE_INFINITY;
}
};
TIntArrayList[] alignment = new TIntArrayList[3];
UnorderedMatcherFactory matcherFactory = new UnorderedMatcherFactory();
HSA hsa = new HSA(w, matcherFactory);
System.out.println("Cost t-s: " + hsa.computeHSA(t, s, alignment));
for (int i = 0; i < alignment[0].size(); ++i) {
System.out.print("(" + alignment[0].get(i) + ","
+ alignment[1].get(i) + "), ");
}
System.out.println("\nCost s-t: " + hsa.computeHSA(s, t, alignment));
for (int i = 0; i < alignment[0].size(); ++i) {
System.out.print("(" + alignment[0].get(i) + ","
+ alignment[1].get(i) + "), ");
}
System.out.println("\nCost t-t: " + hsa.computeHSA(t, t, alignment));
for (int i = 0; i < alignment[0].size(); ++i) {
System.out.print("(" + alignment[0].get(i) + ","
+ alignment[1].get(i) + "), ");
}
System.out.println("\nCost s-s: " + hsa.computeHSA(s, s, alignment));
for (int i = 0; i < alignment[0].size(); ++i) {
System.out.print("(" + alignment[0].get(i) + ","
+ alignment[1].get(i) + "), ");
}
long start = System.currentTimeMillis();
int iters = 10000;
for (int i = 0; i < iters; ++i) {
hsa.computeHSA(t, s);
}
System.out.println("\nTime of " + iters + " iterations for t-s: "
+ (System.currentTimeMillis() - start));
start = System.currentTimeMillis();
for (int i = 0; i < iters; ++i) {
hsa.computeHSA(s, t);
}
System.out.println("\nTime of " + iters + " iterations for s-t: "
+ (System.currentTimeMillis() - start));
start = System.currentTimeMillis();
for (int i = 0; i < iters; ++i) {
hsa.computeHSA(t, t);
}
System.out.println("\nTime of " + iters + " iterations for t-t: "
+ (System.currentTimeMillis() - start));
start = System.currentTimeMillis();
for (int i = 0; i < iters; ++i) {
hsa.computeHSA(s, s);
}
System.out.println("\nTime of " + iters + " iterations for s-s: "
+ (System.currentTimeMillis() - start));
String flowType;
// if (matcherFactory.make(3, 3) instanceof
// UnorderedBipartiteCavityMatcherMinCostMaxFlow) {
// flowType = "min-cost-max-flow";
// } else
flowType = "min-cost-flow";
System.out.println("Flow type: " + flowType);
}
private static void testHsa() {
// t:
// 5-a
// |
// 4-b
// / | \
// / | \
// 2-c 0-h 3-g
// / \
// / \
// 8-i 1-d
// / \
// 6-e 7-f
// s:
// 0-a
// |
// 3-b
// / \
// / \
// 2-j 4-d
// / / | \
// 1-c / | \
// | 7-f 5-i 6-e
// 8-k
// Matchings: (5,0)-a, (4,3)-b, (2,1)-c, (1,4)-d, (6,6)-e, (7,7)-f.
// Cost: -3*6 = -18.
// Smoothings: t - 0, s - 2. Cost: 1.1 + 1.2 = 2.3.
// Prunings: t - (4,3), (0,8). s - (1,8), (4,5). Cost: 4.
// Total alignment cost: -11.7
int[][] tEdges = { { 4, 1, 8 }, { 0, 7, 6 }, { 4 }, { 4 },
{ 5, 3, 0, 2 }, { 4 }, { 1 }, { 1 }, { 0 } };
double[][] tWeights = { { 40, 10, 1 }, { 10, 70, 60 }, { 40 }, { 40 },
{ 50, 1, 10, 20 }, { 40 }, { 10 }, { 10 }, { 10 } };
double[] tSmooth = { 1.1, 2, 2, 2, 2, 2, 2, 2, 2 };
NodeLabel[] tLabels = { new NodeLabel("h"), new NodeLabel("d"),
new NodeLabel("c"), new NodeLabel("g"), new NodeLabel("b"),
new NodeLabel("a"), new NodeLabel("e"), new NodeLabel("f"),
new NodeLabel("i"), };
Tree t = new Tree(tEdges, tWeights, tSmooth, tLabels);
int[][] sEdges = { { 3 }, { 2, 8 }, { 3, 1 }, { 0, 4, 2 },
{ 3, 6, 5, 7 }, { 4 }, { 4 }, { 4 }, { 1 } };
double[][] sWeights = { { 30 }, { 20, 1 }, { 30, 10 }, { 10, 40, 20 },
{ 30, 60, 1, 70 }, { 40 }, { 40 }, { 40 }, { 10 } };
double[] sSmooth = { 2, 2, 1.2, 2, 2, 2, 2, 2, 2 };
NodeLabel[] sLabels = { new NodeLabel("a"), new NodeLabel("c"),
new NodeLabel("j"), new NodeLabel("b"), new NodeLabel("d"),
new NodeLabel("i"), new NodeLabel("e"), new NodeLabel("f"),
new NodeLabel("k"), };
CostFunction w = new CostFunction() {
@Override
public double cost(NodeLabel l1, NodeLabel l2) {
if (Arrays.equals(l1.getLabelValue().toArray(), l2
.getLabelValue().toArray()))
return -3;
else
return 5;
}
};
Tree s = new Tree(sEdges, sWeights, sSmooth, sLabels);
checkTrees(t, w, s);
// checking for rooting problems:
final double INF = MathOperations.INFINITY;
w = new CostFunction() {
@Override
public double cost(NodeLabel l1, NodeLabel l2) {
if (Arrays.equals(l1.getLabelValue().toArray(), l2
.getLabelValue().toArray()))
return -3;
else if (l1.getLabelValue().get(0) == 'b'
|| l2.getLabelValue().get(0) == 'b') {
return INF;
} else
return 5;
}
};
tWeights = new double[][] { { INF, 10, 1 }, { INF, 70, 60 }, { INF },
{ INF }, { 50, 1, 10, 20 }, { INF }, { INF }, { INF }, { INF } };
t = new Tree(tEdges, tWeights, tSmooth, tLabels);
sWeights = new double[][] { { INF }, { INF, 1 }, { INF, 10 },
{ 10, 40, 20 }, { INF, 60, 1, 70 }, { INF }, { INF }, { INF },
{ INF } };
s = new Tree(sEdges, sWeights, sSmooth, sLabels);
checkTrees(t, w, s);
}
private static void checkTrees(Tree t, CostFunction w, Tree s) {
TIntArrayList[] alignment = new TIntArrayList[3];
MatcherFactory matcherFactory = new UnorderedMatcherFactory();
HSA hsa = new HSA(w, matcherFactory);
double cost = hsa.computeHSA(t, s, alignment);
System.out.println("Cost: " + cost);
for (int i = 0; i < alignment[0].size(); ++i) {
System.out.print("(" + alignment[0].get(i) + ","
+ alignment[1].get(i) + "), ");
}
System.out.println();
// checking for pooling problems:
double currCost;
int i = 0;
int iters = 0;
for (; i < iters; ++i) {
hsa = new HSA(w, matcherFactory);
currCost = hsa.computeHSA(t, s);
if (!MathOperations.equals(cost, currCost)) {
System.err.println("Different costs at repeated runs: " + cost
+ ", " + currCost);
break;
}
}
if (i == iters) {
System.out.println("Repeated runs OK.");
}
// checking for symmetry problems:
alignment = new TIntArrayList[3];
currCost = hsa.computeHSA(s, t, alignment);
System.out.println("Cost: " + cost);
for (int j = 0; j < alignment[0].size(); ++j) {
System.out.print("(" + alignment[0].get(j) + ","
+ alignment[1].get(j) + "), ");
}
System.out.println();
// checking for pooling problems:
i = 0;
for (; i < iters; ++i) {
hsa = new HSA(w, matcherFactory);
currCost = hsa.computeHSA(s, t);
if (!MathOperations.equals(cost, currCost)) {
System.err.println("Different costs at repeated runs: " + cost
+ ", " + currCost);
break;
}
}
if (i == iters) {
System.out.println("Repeated runs OK.");
}
}
private static void testBipartiteMatching() {
// Full matching: (0,2), (1,0); costs: -20 + -30 + 2 + 1 + 1 = -46.
// 0, 0: 0.0 *** (1, 2),
// 0, 1: -25.0 *** (1, 0), (2, 2),
// 0, 2: -26.0 *** Cost: -16.0, Matching: (1, 0),
// 0, 3: -25.0 *** Cost: -19.0, Matching: (1, 0), (2, 2),
// 1, 0: -16.0 *** Cost: -6.0, Matching: (0, 2),
// 1, 1: -23.0 *** Cost: -17.0, Matching: (0, 2), (2, 0),
// 1, 2: -6.0 *** Cost: 4.0, Matching: (0, 1), (2, 0),
// 1, 3: -23.0 *** Cost: -17.0, Matching: (0, 2), (2, 0),
// 2, 0: -23.0 *** Cost: -16.0, Matching: (0, 2), (1, 3),
// 2, 1: -49.0 *** Cost: -46.0, Matching: (0, 2), (1, 0),
// 2, 2: -32.0 *** Cost: -25.0, Matching: (0, 1), (1, 0),
// 2, 3: -49.0 *** Cost: -46.0, Matching: (0, 2), (1, 0),
BipartiteCavityMatcher matcher = new UnorderedBipartiteCavityMatcherMinCostFlow(
3, 4);
matcher.setAllMatchCost(0, -4, -3, -20, -2);
matcher.setAllMatchCost(1, -30, -4, -4, -4);
matcher.setAllMatchCost(2, -4, 4, 4, 4);
matcher.setAllDelCostX(5, 5, 2);
matcher.setAllDelCostY(5, 1, 5, 1);
checkCavityMatchings(matcher);
// check forced matchings:
double inf = Double.POSITIVE_INFINITY;
matcher.setDelCostX(0, inf);
matcher.setDelCostY(2, inf);
System.out.println("\n");
System.out.println("Checking forced match of x0 and y2: ");
checkCavityMatchings(matcher);
matcher.setDelCostY(2, 5);
matcher.setDelCostY(0, inf);
System.out.println("\n");
System.out.println("Checking forced match of x0 and y0: ");
checkCavityMatchings(matcher);
matcher.setDelCostX(0, 5);
matcher.setDelCostY(0, 5);
System.out.println("\n");
System.out.println("Checking forced match of x2: ");
matcher.setDelCostX(2, inf);
checkCavityMatchings(matcher);
System.out.println("\n");
System.out.println("Checking group replacement: ");
matcher = new UnorderedBipartiteCavityMatcherMinCostFlow(4, 3);
matcher.setAllMatchCost(0, -4, -30, -4);
matcher.setAllMatchCost(1, -3, -4, 4);
matcher.setAllMatchCost(2, -20, -4, 4);
matcher.setAllMatchCost(3, -2, -4, 4);
matcher.setAllDelCostX(5, 1, 5, 1);
matcher.setAllDelCostY(5, 5, 2);
checkCavityMatchings(matcher);
}
private static void checkCavityMatchings(BipartiteCavityMatcher matcher) {
matcher.clearCavityMatchings();
matcher.setDelCostX(0, matcher.getDelCostX(0)); // forcing an update
matcher.minCostMatching();
System.out.println(matcher);
for (int x = 0; x < matcher.getXSize(); ++x) {
for (int y = 0; y < matcher.getYSize(); ++y) {
System.out.println(x + ", " + y + ": "
+ matcher.minCostCavityMatching(x, y) + " *** "
+ matcher);
}
}
matcher.clearCavityMatchings();
matcher.processAllCavityMatchingY(0);
System.out.println("\n");
for (int y = 0; y < matcher.getYSize(); ++y) {
System.out.println(0 + ", " + y + ": "
+ matcher.minCostCavityMatching(0, y));
}
System.out.println("\n");
matcher.processAllCavityMatchingY(1);
for (int y = 0; y < matcher.getYSize(); ++y) {
System.out.println(1 + ", " + y + ": "
+ matcher.minCostCavityMatching(1, y));
}
System.out.println("\n");
matcher.clearCavityMatchings();
matcher.processAllPairsCavityMatching();
for (int x = 0; x < matcher.getXSize(); ++x) {
for (int y = 0; y < matcher.getYSize(); ++y) {
System.out.println(x + ", " + y + ": "
+ matcher.minCostCavityMatching(x, y));
}
}
System.out.println("\n");
System.out.print(matcher.minCostMatching());
}
private static void testMinCostFlow() {
int[][] neighbors = { { 1, 2 }, { 3, 4, 0 }, { 3, 4, 0 }, { 5, 1, 2 },
{ 5, 1, 2 }, { 3, 4 } };
double inf = Double.POSITIVE_INFINITY;
double[][] weights = { { 0, 0 }, { -6, -1, inf }, { -10, -6, inf },
{ 0, inf, inf }, { 0, inf, inf }, { inf, inf }, };
int[][] capacities = { { 11, 11 }, { 10, 10, 0 }, { 10, 10, 0 },
{ 11, 0, 0 }, { 11, 0, 0 }, { 0, 0 }, };
CostFlowNetwork network = new CostFlowNetwork(neighbors, weights,
capacities, 0, 5);
network.computeMinCostFlow();
}
private static void testShortestPaths() {
int[][] neighbors = { { 1, 5 }, { 2 }, { 3, 5 }, {}, { 3 }, { 4, 1 } };
double[][] weights = { { 4, 2 }, { 3 }, { 0, 0 }, {}, { 3 }, { 2, 1 } };
WeightedGraph graph = new WeightedGraph(neighbors, weights);
ShortestPaths sssp = new ShortestPaths();
double[] delta = sssp.computeShortestPaths(graph, 0);
}
public static void testTree() {
int[][] neighbors = { { 4 }, { 6 }, { 4 }, { 5 }, { 2, 6, 7, 0 },
{ 6, 3, 9, 8 }, { 1, 4, 5 }, { 4 }, { 5 }, { 5 }, };
Tree tree = new Tree(neighbors, null, null, null);
}
public static void testRNAExample() {
RnaSpecificTypeRelatedTreePruning prune = new RnaSpecificTypeRelatedTreePruning(
false, true);
RnaSpecificSmoothCost smooth = new RnaSpecificSmoothCost();
SequenceAlignment aligner = new AffineGapGlobalSequenceAlignmentNoMatrix(
10, 10, RnaAlphabet.getInstance(), new AffineGapScoringMatrix(
"matrix" + File.separator
+ "interval-RIBOSUM85-60.matrix",
RnaAlphabet.getInstance()));
RnaSpecificCostFunction cost = new RnaSpecificCostFunction(
new ScoringMatrix("matrix" + File.separator
+ "bp-RIBOSUM85-60.matrix",
RNABasePairAlphabet.getInstance()), aligner);
RNASpecificTree t = new RNASpecificTree();
// t.buildFromViennaFormat("GAGCCGUAUGCGAUGAAAGUCGCACGUACGGUUC".toCharArray(),
// "((((((((.((((((..))))))...))))))))".toCharArray());
// t.buildFromViennaFormat("GAGCCGUAUGCGAUGAAAGUCGCACGUACGGUUC".toCharArray(),
// "((((((((.((((((..))))))...))))))))".toCharArray());
// t.buildFromViennaFormat("GAGCCGUAUGCGAUGAAAGUCGCACGUACGGUUC".toCharArray(),
// "((((((((.((((((..))))))...))))))))".toCharArray());
t.buildFromViennaFormat(
"NNNUAUAGUUUGAGUUCGAUUGCGCUUCGUAUGUUGCGUCUACGUAAAAACGCUCAGUUUAAAUUAUAACUGCAAAAAAUAAUAACAAUUCUUACGCUUUAGCUGCCUAAUAAGCGCUUAACGUAGAUCCUCCCAGGAUCGUCCAUGUUCUGGAUCUGGGUCCUAAAUUUAGUGGACUUACGCUCAAAGCUUCCACCUGGAGUUGCGAGAAGAGACUAAUCAGGUUAGUCAUUGCUGGGUGCCCUGUCAUACGGCGUUUGCAAUGAUGAAAUUUAAAUAGUAUGAAUAUGAGCGUAGAUAUCCGAGGGGCAAUAUGCUUAGACGCNNN"
.toCharArray(),
"...((((.......((((((..(((((........((((..........))))...........................................(.((((.(((.......))).)))).)..(((((..(((((((.......)))))))....)))))....................(((..((((.........))))..)))..................(((((...(((.((((.........))))).))..)))))......................))))).....)).))))......))))............"
.toCharArray());
prune.calculateCost(t);
smooth.calculateCost(t);
RNASpecificTree s = new RNASpecificTree();
// s.buildFromViennaFormat("GGACCGAUGGUAGUGUCUUCGGAUGCGAGAGUAGGUC".toCharArray(),
// ".((((((((((((((((....))))))))))))))))".toCharArray());
// s.buildFromViennaFormat("GAGGGGAUGAAAAUCCCCUCGAGGGGAUGAAAAUCCCCUCGAGGGGAUGAAAAUCCCCUCGAGGGGAUGAAAUCCCCUC".toCharArray(),
// "(((((((......)))))))(((((((......)))))))(((((((......)))))))(((((((.....)))))))".toCharArray());
// s.buildFromViennaFormat("CCGGgGGAUCACCACGGCGGgGGAUCACCACGG".toCharArray(),
// "((((.((....))))))(((.((....))))).".toCharArray());
s.buildFromViennaFormat(
"GGCGACACGGAUUCCAGUGCAUAUCUUAGUGAUACUCCAGUUAACUCCAUACUUUCCCUGCAAUACGCUAUUCGCCUCAGAUGUAUUUGGGUGGUUGCUCCACUAAAGCCCAGGAAUAUCCAGCCAGUUACAUUUGAGGCCAUUUGGGCUUAAGCGUAUUCCAUGGAAAGUUUUCUCCCCACAUUUCGGAAAUUAAAUUCCGAGCCAGCAAGAAAAUCUUCUCUGUUACAAUUUGACAUGGCUAAAAACUGUACUAAUCAAAAUGAAAAAUGUUUCUCUUGGGCGUAAUCUCAUACAAUGAUUACCCUUAAAGAUCGAACAUUUAAACAAUAAUAUUUGAUAUGAUAUUUUCAAUUUCUAUGCUAUGCCAAAGUGUCUGACAUAAUCAAACAUUUGCACAUUCUUUGACCAAGAAUAGUCAGCAAAUUGUAUUUUCAAUCAAUGCAGACCAUAUGUUCCAGUUUCGGAGAUUUUUUGCUGCCAAACGGAAUACUUAUAAAAACCCACAUUCUAUUUACAUCACUAAGAAGAGCAUUGCAAUCUGUUUAGCC"
.toCharArray(),
"((((..(((((((.((((((...((((((((((.................(((((((.((.((((((((....((((((((((((.((((...((((.(((...........))).)).))..)))).))))))))))))............)))))))).)).)))))))........................(((((((...((((((((((.((((.((((..((((((((.((((((((.(((.((.....((((((.....((((((((.(((((((.(((((...........))))))))..))))..))))))))...........))))))..........)).))).)).)))))).))))(.(((((.(((.........((((((((((((.......))))).))..)))))...............)))))))))..))))..))))....)))))))))))))).....)))))))..........................))))))))))...)))))).)))))))..))))"
.toCharArray());
prune.calculateCost(s);
smooth.calculateCost(s);
TIntArrayList[] alignment = new TIntArrayList[3];
MatcherFactory matcherFactory = new UnorderedMatcherFactory();
// MatcherFactory matcherFactory = new OrderedMatcherFactory();
HSA hsa = new HSA(cost, matcherFactory);
System.out.println("Cost: " + hsa.computeHSA(s, t, alignment));
for (int i = 0; i < alignment[0].size(); ++i) {
System.out.print("(" + alignment[0].get(i) + ","
+ alignment[1].get(i) + "), ");
}
// System.out.println("\nCost: " + hsa.computeHSA(s, t,alignment));
System.out.println("\nCost: " + hsa.computeHSA(t, s, alignment));
for (int i = 0; i < alignment[0].size(); ++i) {
System.out.print("(" + alignment[0].get(i) + ","
+ alignment[1].get(i) + "), ");
}
long start = System.currentTimeMillis();
int iters = 30;
for (int i = 0; i < iters; ++i) {
hsa.computeHSA(t, s);
}
System.out.println("\nTime of " + iters + " iterations: "
+ (System.currentTimeMillis() - start));
String flowType;
// if (matcherFactory.make(3, 3) instanceof
// UnorderedBipartiteCavityMatcherMinCostMaxFlow) {
// flowType = "min-cost-max-flow";
// } else
flowType = "min-cost-flow";
System.out.println("Flow type: " + flowType);
// System.out.println("\nCost: " + hsa.computeHSA(t, s, alignment));
// for (int i=0; i<alignment[0].size(); ++i){
// System.out.print("(" + alignment[0].get(i) +"," + alignment[1].get(i)
// + "), ");
// }
}
}