Skip to content

Commit 79a0ba6

Browse files
benedictpatenbenedictpaten
committed
Changing insert rules
1 parent 8197e0e commit 79a0ba6

File tree

4 files changed

+39
-23
lines changed

4 files changed

+39
-23
lines changed

‎impl/referenceProblem2.c‎

Lines changed: 31 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -511,7 +511,7 @@ static stList *getInsertPointsNext(int64_t n, stList *edges, reference *ref) {
511511
return insertPoints;
512512
}
513513

514-
static void getInsertionPoints(int64_t n, stList *previousEdges, stList *nextEdges, reference *ref, refAdjList *aL,
514+
static void getInsertionPoints(int64_t n, stList *previousEdges, stList *nextEdges, reference *ref, refAdjList *aL, refAdjList *dAL,
515515
stList *insertPoints) {
516516
stList *previousInsertPoints = getInsertPointsPrevious(n, previousEdges, ref);
517517
stList_appendAll(insertPoints, previousInsertPoints);
@@ -532,8 +532,24 @@ static void getInsertionPoints(int64_t n, stList *previousEdges, stList *nextEdg
532532
}
533533
if (reference_getFirst(ref, insertPoint_adjNode(iPP)) == reference_getFirst(ref, insertPoint_adjNode(iPN))) {
534534
bool equivalentInsertPoints = llabs(reference_getNext(ref, insertPoint_adjNode(iPP))) == llabs(insertPoint_adjNode(iPN)) ? 0 : 1; //I don't think the absolute values are necessary
535-
if (refAdjList_getWeight(aL, n, insertPoint_adjNode(iPP)) > refAdjList_getWeight(aL, -n,
536-
insertPoint_adjNode(iPN))) {
535+
double nWL = refAdjList_getWeight(aL, insertPoint_adjNode(iPP), n); //new weight left
536+
double nWR = refAdjList_getWeight(aL, -n, insertPoint_adjNode(iPN)); //new weight right
537+
assert(nWL > 0);
538+
assert(nWR > 0);
539+
bool left = nWL > nWR;
540+
if(equivalentInsertPoints) {
541+
assert(reference_getNext(ref, insertPoint_adjNode(iPP)) != INT64_MAX);
542+
assert(reference_getPrevious(ref, insertPoint_adjNode(iPN)) != INT64_MAX);
543+
double eWL = refAdjList_getWeight(dAL, insertPoint_adjNode(iPP), reference_getNext(ref, insertPoint_adjNode(iPP))); //existing weight left
544+
double eWR = refAdjList_getWeight(dAL, -reference_getPrevious(ref, insertPoint_adjNode(iPN)), insertPoint_adjNode(iPN)); //existing weight right
545+
if(eWL == 0) {
546+
left = eWR > 0 || nWL > nWR;
547+
}
548+
else {
549+
left = !(nWR == 0 || nWR > nWL);
550+
}
551+
}
552+
if (left) {
537553
stList_append(
538554
insertPoints,
539555
insertPoint_construct(n, insertPoint_adjNode(iPP), 1,
@@ -559,16 +575,16 @@ static void getInsertionPoints(int64_t n, stList *previousEdges, stList *nextEdg
559575
stList_destruct(nextInsertPoints);
560576
}
561577

562-
static insertPoint *getABestInsertNode(int64_t n, refAdjList *aL, reference *ref) {
578+
static insertPoint *getABestInsertNode(int64_t n, refAdjList *aL, refAdjList *dAL, reference *ref) {
563579
assert(!reference_inGraph(ref, n));
564580
insertPoint *bestIP = NULL;
565581
//Get list of edges to nodes already in the reference
566582
stList *previousEdges = getRelevantEdges(aL, ref, n);
567583
stList *nextEdges = getRelevantEdges(aL, ref, -n);
568584
//Now do the hardwork of determining the best insertion point
569585
stList *insertPoints = stList_construct3(0, free);
570-
getInsertionPoints(n, previousEdges, nextEdges, ref, aL, insertPoints);
571-
getInsertionPoints(-n, nextEdges, previousEdges, ref, aL, insertPoints);
586+
getInsertionPoints(n, previousEdges, nextEdges, ref, aL, dAL, insertPoints);
587+
getInsertionPoints(-n, nextEdges, previousEdges, ref, aL, dAL, insertPoints);
572588
//Get best insertion point
573589
for (int64_t i = 0; i < stList_length(insertPoints); i++) {
574590
insertPoint *iP = stList_get(insertPoints, i);
@@ -588,9 +604,9 @@ static insertPoint *getABestInsertNode(int64_t n, refAdjList *aL, reference *ref
588604
return bestIP;
589605
}
590606

591-
static void insertNode(int64_t n, refAdjList *aL, reference *ref) {
607+
static void insertNode(int64_t n, refAdjList *aL, refAdjList *dAL, reference *ref) {
592608
if (!reference_inGraph(ref, n)) { //Have a node to add
593-
insertPoint *bestIP = getABestInsertNode(n, aL, ref);
609+
insertPoint *bestIP = getABestInsertNode(n, aL, dAL, ref);
594610
if (bestIP == NULL) { //Make up a location
595611
st_logDebug("Got a node with no edges linking it into the graph\n");
596612
reference_insertNode(ref, reference_getFirstOfInterval(ref, 0), n);
@@ -699,13 +715,13 @@ static bool connectedNodes_empty(connectedNodes *cN) {
699715
return stSortedSet_size(cN->byNode) == 0;
700716
}
701717

702-
static insertPoint *connectedNodes_popBestInsert(connectedNodes *cN, refAdjList *aL, reference *ref, double wiggle) {
718+
static insertPoint *connectedNodes_popBestInsert(connectedNodes *cN, refAdjList *aL, refAdjList *dAL, reference *ref, double wiggle) {
703719
assert(wiggle <= 1.0);
704720
int64_t i = 0;
705721
while (1) {
706722
ConnectedNodeEdge *cNE = stSortedSet_getLast(cN->byWeight);
707723
stSortedSet_remove(cN->byWeight, cNE);
708-
insertPoint *iP = getABestInsertNode(refEdge_to((refEdge *) cNE), aL, ref);
724+
insertPoint *iP = getABestInsertNode(refEdge_to((refEdge *) cNE), aL, dAL, ref);
709725
assert(iP != NULL);
710726
cNE->inconsistentAdjacencyWeight = cNE->weightOfEdgesInGraph - insertPoint_score(iP);
711727
assert(insertPoint_score(iP) / cNE->weightOfEdgesInGraph <= 1.0001);
@@ -725,14 +741,14 @@ static insertPoint *connectedNodes_popBestInsert(connectedNodes *cN, refAdjList
725741
* Actual algorithms to make the reference.
726742
*/
727743

728-
void makeReferenceGreedily2(refAdjList *aL, reference *ref, double wiggle) {
744+
void makeReferenceGreedily2(refAdjList *aL, refAdjList *dAL, reference *ref, double wiggle) {
729745
assert(reference_getIntervalNumber(ref) > 0 || refAdjList_getNodeNumber(aL) == 0);
730746
connectedNodes *cN = connectedNodes_construct(aL, ref);
731747
//Iterate over the nodes to check any nodes that are not in the reference
732748
int64_t i = 0;
733749
for (int64_t n = 1; n <= refAdjList_getNodeNumber(aL); n++) {
734750
while (!connectedNodes_empty(cN)) {
735-
insertPoint *iP = connectedNodes_popBestInsert(cN, aL, ref, wiggle);
751+
insertPoint *iP = connectedNodes_popBestInsert(cN, aL, dAL, ref, wiggle);
736752
assert(iP != NULL);
737753
reference_insertNode2(ref, iP);
738754
connectedNodes_addNode(cN, insertPoint_node(iP));
@@ -741,7 +757,7 @@ void makeReferenceGreedily2(refAdjList *aL, reference *ref, double wiggle) {
741757
}
742758
if(!reference_inGraph(ref, n)) {
743759
i++;
744-
insertNode(n, aL, ref);
760+
insertNode(n, aL, dAL, ref);
745761
connectedNodes_addNode(cN, n);
746762
connectedNodes_addNode(cN, -n);
747763
}
@@ -751,14 +767,14 @@ void makeReferenceGreedily2(refAdjList *aL, reference *ref, double wiggle) {
751767
connectedNodes_destruct(cN);
752768
}
753769

754-
void updateReferenceGreedily(refAdjList *aL, reference *ref, int64_t permutations) {
770+
void updateReferenceGreedily(refAdjList *aL, refAdjList *dAL, reference *ref, int64_t permutations) {
755771
for (int64_t i = 0; i < permutations; i++) {
756772
for (int64_t j = 1; j <= refAdjList_getNodeNumber(aL); j++) {
757773
int64_t n = st_randomInt(1, refAdjList_getNodeNumber(aL) + 1);
758774
assert(reference_inGraph(ref, n));
759775
reference_removeNode(ref, n);
760776
if (!reference_inGraph(ref, n)) {
761-
insertNode(n, aL, ref);
777+
insertNode(n, aL, dAL, ref);
762778
assert(reference_inGraph(ref, n));
763779
}
764780
}

‎inc/stReferenceProblem2.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -108,9 +108,9 @@ void reference_log(reference *ref);
108108
* Reference algorithms
109109
*/
110110

111-
void makeReferenceGreedily2(refAdjList *aL, reference *ref, double wiggle);
111+
void makeReferenceGreedily2(refAdjList *aL, refAdjList *dAL, reference *ref, double wiggle);
112112

113-
void updateReferenceGreedily(refAdjList *aL, reference *ref, int64_t permutations);
113+
void updateReferenceGreedily(refAdjList *aL, refAdjList *dAL, reference *ref, int64_t permutations);
114114

115115
/*
116116
* Create a topological sort of each reference interval, trying to place nodes that are connected by direct adjacencies next to one another.

‎tests/referenceProblemTest2.c‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -182,12 +182,12 @@ static void testMakeReferenceGreedily(CuTest *testCase) {
182182
for (int64_t i = 0; i < 100; i++) {
183183
setup();
184184
time_t startTime = time(NULL);
185-
makeReferenceGreedily2(aL, ref, 0.99);
185+
makeReferenceGreedily2(aL, dAL, ref, 0.99);
186186
int64_t badAdjacencyCount = getBadAdjacencyCount(dAL, ref);
187187
st_logInfo("Greedy it took %" PRIi64 " seconds, score: %Lf of possible: %Lf, bad adjacency count: %" PRIi64 "\n", time(NULL) - startTime,
188188
getReferenceScore(aL, ref), refAdjList_getMaxPossibleScore(aL), badAdjacencyCount);
189189
checkIsValidReference(testCase);
190-
updateReferenceGreedily(aL, ref, 10);
190+
updateReferenceGreedily(aL, dAL, ref, 10);
191191
long double greedyPermutationScore = getReferenceScore(aL, ref);
192192
int64_t badAdjacencyCountGreedyPermutations = getBadAdjacencyCount(dAL, ref);
193193
st_logInfo("Greedy with update permutations, it took %" PRIi64 " seconds, score: %Lf of possible: %Lf, bad adjacency count: %" PRIi64 "\n",
@@ -268,8 +268,8 @@ static void testADBDCExample(CuTest *testCase) {
268268
fn(theta, _3D, _5B, adjacencyLength, DL, BL, n);
269269
fn(theta, _3D, _3B, adjacencyLength, DL, BL, n);
270270

271-
makeReferenceGreedily2(aL, ref, 0.99);
272-
updateReferenceGreedily(aL, ref, 100);
271+
makeReferenceGreedily2(aL, aL, ref, 0.99);
272+
updateReferenceGreedily(aL, aL, ref, 100);
273273
st_logInfo("Running reference example problem, score: %f of possible: %f\n", getReferenceScore(aL, ref),
274274
refAdjList_getMaxPossibleScore(aL));
275275
reference_log(ref);

‎tests/testBin/referenceMedianProblemTest2.c‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -86,8 +86,8 @@ int main(int argc, char *argv[]) {
8686
* Compute the ordering
8787
*/
8888

89-
makeReferenceGreedily2(aL, ref, 0.99);
90-
updateReferenceGreedily(aL, ref, permutationNumber);
89+
makeReferenceGreedily2(aL, aL, ref, 0.99);
90+
updateReferenceGreedily(aL, aL, ref, permutationNumber);
9191

9292
/*
9393
* Print out the median genome,

0 commit comments

Comments
 (0)