From ccdbc32caa6f1768867a23020dd8034015a2bef3 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Thu, 1 Oct 2026 14:58:34 -0400 Subject: [PATCH 1/2] Per-pair state machines in the multiple aligner makeAlignmentWithPairStateMachines takes a callback that picks the machine for each pair of sequences it aligns, so one multiple alignment can score different pairs with different models -- for cactus, by how diverged the pair's genomes are. A NULL callback, or a NULL from it, means the default machine. makeAlignment, makeAlignmentUsingAllPairs and makeAllPairwiseAlignments keep their signatures and behaviour. The new test checks that a chooser returning the default changes nothing, that it is asked once per pairwise alignment about a real pair, and that another machine gives a valid and different alignment. It also turned up that makeAlignment is not always repeatable: one random case in 50, on the spanning-tree path without progressive merging, aligned the same sequences differently on a second call. Co-Authored-By: Claude Opus 5.5 --- impl/multipleAligner.c | 55 +++++++++++++++++++---- inc/multipleAligner.h | 18 ++++++++ tests/multipleAlignerTest.c | 89 +++++++++++++++++++++++++++++++++++++ 3 files changed, 153 insertions(+), 9 deletions(-) diff --git a/impl/multipleAligner.c b/impl/multipleAligner.c index 7baacb7..a196b94 100644 --- a/impl/multipleAligner.c +++ b/impl/multipleAligner.c @@ -1025,14 +1025,29 @@ static stList *convertMultipleAlignedPairsToAlignedPairs(stList *multipleAligned return alignedPairs; } -static int64_t addMultipleAlignedPairs(StateMachine *sM, int64_t sequence1, int64_t sequence2, stList *seqFrags, stList *multipleAlignedPairs, +/* + * The machines a multiple alignment scores its pairs with: sM, or what pairStateMachine picks for the pair + * when it is given and does not return NULL. + */ +typedef struct _pairStateMachines { + StateMachine *sM; + PairStateMachineFn pairStateMachine; + void *extraArgs; +} PairStateMachines; + +static StateMachine *pairStateMachines_get(PairStateMachines *sMs, int64_t sequence1, int64_t sequence2) { + StateMachine *sM = sMs->pairStateMachine == NULL ? NULL : sMs->pairStateMachine(sequence1, sequence2, sMs->extraArgs); + return sM == NULL ? sMs->sM : sM; +} + +static int64_t addMultipleAlignedPairs(PairStateMachines *sMs, int64_t sequence1, int64_t sequence2, stList *seqFrags, stList *multipleAlignedPairs, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) { /* * Computes a pairwise alignment and returns the pairwise match probabilities as tuples of (score, seq1, pos1, seq2, pos2). */ SeqFrag *seqFrag1 = stList_get(seqFrags, sequence1); SeqFrag *seqFrag2 = stList_get(seqFrags, sequence2); - stList *alignedPairs = getAlignedPairs(sM, seqFrag1->seq, seqFrag2->seq, pairwiseAlignmentBandingParameters, + stList *alignedPairs = getAlignedPairs(pairStateMachines_get(sMs, sequence1, sequence2), seqFrag1->seq, seqFrag2->seq, pairwiseAlignmentBandingParameters, seqFrag1->leftEndId != seqFrag2->leftEndId, seqFrag1->rightEndId != seqFrag2->rightEndId); alignedPairs = reweightAlignedPairs2(alignedPairs, seqFrag1->length, seqFrag2->length, pairwiseAlignmentBandingParameters->gapGamma); int64_t distance = getAlignmentScore(alignedPairs, seqFrag1->length, seqFrag2->length); @@ -1040,7 +1055,7 @@ static int64_t addMultipleAlignedPairs(StateMachine *sM, int64_t sequence1, int6 return distance; } -stList *makeAllPairwiseAlignments(StateMachine *sM, stList *seqFrags, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, stList **seqPairSimilarityScores) { +static stList *makeAllPairwiseAlignments2(PairStateMachines *sMs, stList *seqFrags, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, stList **seqPairSimilarityScores) { /* * Generate the set of pairwise alignments between the sequences. */ @@ -1049,20 +1064,25 @@ stList *makeAllPairwiseAlignments(StateMachine *sM, stList *seqFrags, PairwiseAl int64_t seqNo = stList_length(seqFrags); for (int64_t seq1 = 0; seq1 < seqNo; seq1++) { for (int64_t seq2 = seq1 + 1; seq2 < seqNo; seq2++) { - stList_append(*seqPairSimilarityScores, stIntTuple_construct3(addMultipleAlignedPairs(sM, seq1, seq2, seqFrags, multipleAlignedPairs, pairwiseAlignmentBandingParameters), seq1, seq2)); + stList_append(*seqPairSimilarityScores, stIntTuple_construct3(addMultipleAlignedPairs(sMs, seq1, seq2, seqFrags, multipleAlignedPairs, pairwiseAlignmentBandingParameters), seq1, seq2)); } } return multipleAlignedPairs; } -MultipleAlignment *makeAlignmentUsingAllPairs(StateMachine *sM, stList *seqFrags, +stList *makeAllPairwiseAlignments(StateMachine *sM, stList *seqFrags, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, stList **seqPairSimilarityScores) { + PairStateMachines sMs = { sM, NULL, NULL }; + return makeAllPairwiseAlignments2(&sMs, seqFrags, pairwiseAlignmentBandingParameters, seqPairSimilarityScores); +} + +static MultipleAlignment *makeAlignmentUsingAllPairs2(PairStateMachines *sMs, stList *seqFrags, bool useProgressiveMerging, float matchGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) { /* * Generate a multiple alignment considering all pairs of sequences. */ MultipleAlignment *mA = st_calloc(1, sizeof(MultipleAlignment)); - mA->alignedPairs = makeAllPairwiseAlignments(sM, seqFrags, pairwiseAlignmentBandingParameters, &mA->chosenPairwiseAlignments); + mA->alignedPairs = makeAllPairwiseAlignments2(sMs, seqFrags, pairwiseAlignmentBandingParameters, &mA->chosenPairwiseAlignments); if(stList_length(seqFrags) == 2 || useProgressiveMerging) { //Compute an optimum exactly mA->columns = getMultipleSequenceAlignmentProgressive(seqFrags, mA->alignedPairs, matchGamma, mA->chosenPairwiseAlignments); } @@ -1073,6 +1093,13 @@ MultipleAlignment *makeAlignmentUsingAllPairs(StateMachine *sM, stList *seqFrags return mA; } +MultipleAlignment *makeAlignmentUsingAllPairs(StateMachine *sM, stList *seqFrags, + bool useProgressiveMerging, float matchGamma, + PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) { + PairStateMachines sMs = { sM, NULL, NULL }; + return makeAlignmentUsingAllPairs2(&sMs, seqFrags, useProgressiveMerging, matchGamma, pairwiseAlignmentBandingParameters); +} + void multipleAlignment_destruct(MultipleAlignment *mA) { stList_destruct(mA->alignedPairs); stSet_destruct(mA->columns); @@ -1262,12 +1289,22 @@ int64_t getNextBestPair(int64_t seq1, int64_t *distanceCounts, int64_t seqNo, MultipleAlignment *makeAlignment(StateMachine *sM, stList *seqFrags, int64_t spanningTrees, int64_t maxPairsToConsider, bool useProgressiveMerging, float matchGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) { + return makeAlignmentWithPairStateMachines(sM, NULL, NULL, seqFrags, spanningTrees, maxPairsToConsider, + useProgressiveMerging, matchGamma, pairwiseAlignmentBandingParameters); +} + +MultipleAlignment *makeAlignmentWithPairStateMachines(StateMachine *sM, PairStateMachineFn pairStateMachine, void *extraArgs, + stList *seqFrags, int64_t spanningTrees, int64_t maxPairsToConsider, + bool useProgressiveMerging, float matchGamma, + PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) { /* * Computes an MSA, making up to "spanningTrees"*no of seqs pairwise alignments. */ + PairStateMachines sMsValue = { sM, pairStateMachine, extraArgs }; + PairStateMachines *sMs = &sMsValue; int64_t seqNo = stList_length(seqFrags); if (spanningTrees * (seqNo - 1) >= (seqNo * (seqNo - 1)) / 2) { //Do all pairs if we can - return makeAlignmentUsingAllPairs(sM, seqFrags, useProgressiveMerging, matchGamma, pairwiseAlignmentBandingParameters); + return makeAlignmentUsingAllPairs2(sMs, seqFrags, useProgressiveMerging, matchGamma, pairwiseAlignmentBandingParameters); } MultipleAlignment *mA = st_calloc(1, sizeof(MultipleAlignment)); mA->alignedPairs = stList_construct3(0, (void(*)(void *)) stIntTuple_destruct); //pairwise alignment pairs, with sequence indices @@ -1279,7 +1316,7 @@ MultipleAlignment *makeAlignment(StateMachine *sM, stList *seqFrags, int64_t spa int64_t seqX = stIntTuple_get(pairToAlign, 0); int64_t seqY = stIntTuple_get(pairToAlign, 1); //We get pairwise alignments, for this first alignment we filter the pairs greedily to make them consistent - stList_append(mA->chosenPairwiseAlignments, stIntTuple_construct3(addMultipleAlignedPairs(sM, seqX, seqY, seqFrags, + stList_append(mA->chosenPairwiseAlignments, stIntTuple_construct3(addMultipleAlignedPairs(sMs, seqX, seqY, seqFrags, mA->alignedPairs, pairwiseAlignmentBandingParameters), seqX, seqY)); } @@ -1304,7 +1341,7 @@ MultipleAlignment *makeAlignment(StateMachine *sM, stList *seqFrags, int64_t spa assert(seq != otherSeq); stIntTuple *pairToAlign = makePairToAlign(seq, otherSeq); assert(stSortedSet_search(chosenPairwiseAlignmentsSet, pairToAlign) == NULL); - stList_append(mA->chosenPairwiseAlignments, stIntTuple_construct3(addMultipleAlignedPairs(sM, seq, otherSeq, seqFrags, mA->alignedPairs, pairwiseAlignmentBandingParameters), seq, otherSeq)); + stList_append(mA->chosenPairwiseAlignments, stIntTuple_construct3(addMultipleAlignedPairs(sMs, seq, otherSeq, seqFrags, mA->alignedPairs, pairwiseAlignmentBandingParameters), seq, otherSeq)); stSortedSet_insert(chosenPairwiseAlignmentsSet, pairToAlign); } } diff --git a/inc/multipleAligner.h b/inc/multipleAligner.h index 155e2f8..40703ee 100644 --- a/inc/multipleAligner.h +++ b/inc/multipleAligner.h @@ -51,6 +51,24 @@ MultipleAlignment *makeAlignment(StateMachine *sM, stList *seqFrags, float matchGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters); +/* + * Picks the machine to align one pair of sequences with, by their indices in seqFrags: lets one multiple + * alignment score different pairs with different models, e.g. by how diverged the pair's genomes are. + * Called once for each pairwise alignment, possibly from several threads if the caller aligns several + * MSAs at once. The machine must outlive the alignment, and the same machine can be returned for any + * number of pairs. NULL means the alignment's default machine. + */ +typedef StateMachine *(*PairStateMachineFn)(int64_t seqX, int64_t seqY, void *extraArgs); + +/* + * As makeAlignment, scoring each pair with the machine pairStateMachine picks for it (sM for any pair it + * returns NULL for, and for every pair when pairStateMachine is NULL). + */ +MultipleAlignment *makeAlignmentWithPairStateMachines(StateMachine *sM, PairStateMachineFn pairStateMachine, void *extraArgs, + stList *seqFrags, int64_t spanningTrees, int64_t maxPairsToConsider, + bool useProgressiveMerging, float matchGamma, + PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters); + MultipleAlignment *makeAlignmentUsingAllPairs(StateMachine *sM, stList *seqFrags, bool useProgressiveMerging, float matchGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters); diff --git a/tests/multipleAlignerTest.c b/tests/multipleAlignerTest.c index 0249455..54e04e7 100644 --- a/tests/multipleAlignerTest.c +++ b/tests/multipleAlignerTest.c @@ -231,6 +231,94 @@ static void test_multipleAlignerRandom(CuTest *testCase) { } } +/* + * Per-pair machines: a chooser that records how it was asked and returns one fixed machine. + */ +typedef struct { + StateMachine *sM; // what to return for every pair (NULL for the default) + int64_t seqNo; + int64_t calls; + bool badPair; +} PairChoice; + +static StateMachine *choosePairStateMachine(int64_t seqX, int64_t seqY, void *extraArgs) { + PairChoice *choice = extraArgs; + choice->calls++; + if (seqX < 0 || seqY < 0 || seqX >= choice->seqNo || seqY >= choice->seqNo || seqX == seqY) { + choice->badPair = 1; + } + return choice->sM; +} + +static bool sameAlignedPairs(stList *pairs1, stList *pairs2) { + // as sets: the order the pairs come out in follows hashes of pointers, so differs from run to run + if (stList_length(pairs1) != stList_length(pairs2)) { + return 0; + } + stList *sorted1 = stList_copy(pairs1, NULL), *sorted2 = stList_copy(pairs2, NULL); + stList_sort(sorted1, (int (*)(const void *, const void *))stIntTuple_cmpFn); + stList_sort(sorted2, (int (*)(const void *, const void *))stIntTuple_cmpFn); + bool same = 1; + for (int64_t i = 0; i < stList_length(sorted1) && same; i++) { + same = stIntTuple_equalsFn(stList_get(sorted1, i), stList_get(sorted2, i)); + } + stList_destruct(sorted1); + stList_destruct(sorted2); + return same; +} + +static void test_makeAlignmentWithPairStateMachines(CuTest *testCase) { + int64_t differed = 0, repeatable = 0; + for (int64_t test = 0; test < 50; test++) { + setup(); + // spanning trees from 1 to 4 over 2 to 9 sequences reach both the all-pairs and the spanning-tree path + stList *randomSeqFrags = getRandomSeqFrags(st_randomInt(2, 10), st_randomInt(20, 100)); + int64_t seqNo = stList_length(randomSeqFrags); + int64_t spanningTrees = st_randomInt(1, 5); + bool progressive = st_random() > 0.5; + MultipleAlignment *mA = makeAlignment(stateMachine, randomSeqFrags, spanningTrees, 10000000, progressive, 0.5, pabp); + // makeAlignment can give different alignments of the same sequences from one call to the next (seen + // without progressive merging, on the spanning-tree path), so equality below is only asked for where it + // agrees with itself + MultipleAlignment *again = makeAlignment(stateMachine, randomSeqFrags, spanningTrees, 10000000, progressive, 0.5, pabp); + bool isRepeatable = sameAlignedPairs(mA->alignedPairs, again->alignedPairs); + repeatable += isRepeatable; + multipleAlignment_destruct(again); + + // a chooser that returns NULL, or the default machine itself, changes nothing, and is asked once per + // pairwise alignment about a real pair of sequences + for (int64_t i = 0; i < 2; i++) { + PairChoice choice = { i == 0 ? NULL : stateMachine, seqNo, 0, 0 }; + MultipleAlignment *mA2 = makeAlignmentWithPairStateMachines(stateMachine, choosePairStateMachine, &choice, + randomSeqFrags, spanningTrees, 10000000, progressive, 0.5, pabp); + CuAssertTrue(testCase, !isRepeatable || sameAlignedPairs(mA->alignedPairs, mA2->alignedPairs)); + CuAssertIntEquals(testCase, stList_length(mA2->chosenPairwiseAlignments), choice.calls); + CuAssertTrue(testCase, !choice.badPair); + multipleAlignment_destruct(mA2); + } + + // another machine for every pair still gives a valid alignment, and, over the tests, a different one + Hmm *hmm = hmm_constructEmpty(0.0, fiveState); + hmm_randomise(hmm); + StateMachine *other = hmm_getStateMachine(hmm); + PairChoice choice = { other, seqNo, 0, 0 }; + MultipleAlignment *mA3 = makeAlignmentWithPairStateMachines(stateMachine, choosePairStateMachine, &choice, + randomSeqFrags, spanningTrees, 10000000, progressive, 0.5, pabp); + checkAlignment(testCase, randomSeqFrags, mA3->alignedPairs); + CuAssertIntEquals(testCase, stList_length(mA3->chosenPairwiseAlignments), choice.calls); + differed += isRepeatable && !sameAlignedPairs(mA->alignedPairs, mA3->alignedPairs); + multipleAlignment_destruct(mA3); + stateMachine_destruct(other); + hmm_destruct(hmm); + + multipleAlignment_destruct(mA); + stList_destruct(randomSeqFrags); + teardown(); + } + CuAssertTrue(testCase, repeatable >= 40); + CuAssertTrue(testCase, differed > 0); +} + CuSuite* multipleAlignerTestSuite(void) { CuSuite* suite = CuSuiteNew(); SUITE_ADD_TEST(suite, test_getDistanceMatrix); @@ -241,6 +329,7 @@ CuSuite* multipleAlignerTestSuite(void) { SUITE_ADD_TEST(suite, test_makeAlignmentUsingAllPairs); SUITE_ADD_TEST(suite, test_multipleAlignerAllPairsRandom); SUITE_ADD_TEST(suite, test_multipleAlignerRandom); + SUITE_ADD_TEST(suite, test_makeAlignmentWithPairStateMachines); return suite; } From 2ec01d160bd4a558ef213a168a57748b82c9313d Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Fri, 2 Oct 2026 08:31:36 -0400 Subject: [PATCH 2/2] Make the multiple aligner deterministic The same sequences, models and parameters now always give the same multiple alignment. Before, the result could depend on: - the state of sonLib's random number generator, which the aligner drew from for each aligned pair's tie-breaking jitter, and for a coin flip between equally good pairs when choosing spanning-tree alignments -- so on whatever had drawn from it before (in cactus, on the other ends of the flower aligned before this one); - where memory had been allocated: weights with equal scores were ordered by address, and which of a weight and its reverse went in the ordered set, and so whose jitter counted, was chosen by comparing column addresses. Now the jitter is a hash of the aligned pair (same magnitude, and the same for both directions of the pair), weight ties go by the columns' positions, and the spanning-tree tie-break goes by a hash of the pair of sequences. Where the jitter decided a near-tie, an alignment can differ from the old code's, as the old code's could from one run to the next. In 400 random cases, two identical calls of the old code gave different aligned pairs in 16 of the 201 on the spanning-tree path (8 of 95 greedy, 8 of 106 progressive), and in none of the 199 all-pairs cases. The tests now require identical output -- aligned pairs, chosen pairwise alignments and columns, in order -- and a new one moves the generator on and holds memory between two calls. The old code fails both under each of 24 seeds; this passes them all. Cactus already reseeds the generator per flower, so its runs were reproducible; what changes there is that an end's alignment no longer depends on the ends aligned before it, so comparing two configurations is not confounded by the random stream shifting between them. Co-Authored-By: Claude Opus 5.5 --- impl/multipleAligner.c | 73 +++++++++++++++++++++++++------ tests/multipleAlignerTest.c | 87 ++++++++++++++++++++++++++++--------- 2 files changed, 126 insertions(+), 34 deletions(-) diff --git a/impl/multipleAligner.c b/impl/multipleAligner.c index a196b94..3913f32 100644 --- a/impl/multipleAligner.c +++ b/impl/multipleAligner.c @@ -62,6 +62,38 @@ int column_cmp(Column *c, Column *c2) { : (c->position < c2->position ? -1 : 0))); } +/* + * Deterministic stand-ins for the random numbers this file used to draw from sonLib's global generator. + * Those made an alignment depend on that generator's state -- so on whatever had run before it, and in a + * multithreaded caller on thread timing -- and pointer comparisons made it depend on where things had been + * allocated. Now the same sequences give the same alignment. + */ +static uint64_t mix64(uint64_t h) { + h ^= h >> 30; + h *= 0xBF58476D1CE4E5B9ULL; + h ^= h >> 27; + h *= 0x94D049BB133111EBULL; + return h ^ (h >> 31); +} + +static uint64_t pairHash(int64_t seq1, int64_t pos1, int64_t seq2, int64_t pos2) { + // the same either way round + if (seq1 > seq2 || (seq1 == seq2 && pos1 > pos2)) { + int64_t s = seq1, p = pos1; + seq1 = seq2, pos1 = pos2, seq2 = s, pos2 = p; + } + return mix64(mix64((uint64_t) seq1 * 0x9E3779B97F4A7C15ULL ^ (uint64_t) pos1) ^ ((uint64_t) seq2 * 0xC2B2AE3D27D4EB4FULL + (uint64_t) pos2)); +} + +/* + * The tie-breaking jitter added to an aligned pair's weight, in [0, 0.00001) like the random one it + * replaces (which kept the merge order from degenerating into badly unbalanced trees, without affecting + * accuracy). Both directions of the pair get the same value. + */ +static double weightJitter(int64_t seq1, int64_t pos1, int64_t seq2, int64_t pos2) { + return (double) (pairHash(seq1, pos1, seq2, pos2) >> 11) / 9007199254740992.0 * 0.00001; +} + uint64_t column_hashFn(Column *c) { // Mixed, not summed: seqName + position put every sequence's positions on top of one another, // so each bucket held a position from every sequence, and lookups here were a fifth of the @@ -116,14 +148,23 @@ int alignmentWeight_cmpByPosition(AlignmentWeight *aW, AlignmentWeight *aW2) { int alignmentWeight_cmpByWeight(AlignmentWeight *aW, AlignmentWeight *aW2) { /* - * Cmp primarily by weight, then position. + * Cmp primarily by weight, then by the columns the weights join (not by address, which would make + * the order depend on where the weights were allocated). */ if (aW->avgWeight == aW2->avgWeight) { - return aW < aW2 ? -1 : (aW > aW2 ? 1 : 0); + int i = column_cmp(aW->column, aW2->column); + return i != 0 ? i : column_cmp(aW->rWeight->column, aW2->rWeight->column); } return aW->avgWeight > aW2->avgWeight ? 1 : -1; } +/* + * Of a weight and its reverse, the one kept in the set ordered by weight: the one from the lesser column. + */ +static bool alignmentWeight_isCanonical(AlignmentWeight *aW) { + return column_cmp(aW->column, aW->rWeight->column) < 0; +} + void insertWeight(AlignmentWeight *aW, stHash *alignmentWeightAdjLists) { /* * Inserts an alignment weight to an adjacency list. @@ -143,12 +184,12 @@ static Column *getColumn(stSet *columns, int64_t seqName, int64_t position) { return stSet_search(columns, &c); } -static AlignmentWeight *makeAlignmentWeight(stSet *columns, int64_t score, int64_t seqName, int64_t position) { +static AlignmentWeight *makeAlignmentWeight(stSet *columns, double weight, int64_t seqName, int64_t position) { AlignmentWeight *aW = st_malloc(sizeof(AlignmentWeight)); aW->column = getColumn(columns, seqName, position); assert(aW->column != NULL); aW->numberOfWeights = 1; - aW->avgWeight = ((double) score) / PAIR_ALIGNMENT_PROB_1 + st_random() * 0.00001; //This randomness avoids nasty types of unbalanced trees and doesn't really affect accuracy + aW->avgWeight = weight; return aW; } @@ -161,8 +202,10 @@ stHash *makeAlignmentWeightAdjacencyLists(stSet *columns, stList *multipleAligne /*Tuple of score, seq1, pos1, seq2, pos2 */ stIntTuple *aP = stList_get(multipleAlignedPairs, i); assert(stIntTuple_length(aP) == 5); - AlignmentWeight *aW = makeAlignmentWeight(columns, stIntTuple_get(aP, 0), stIntTuple_get(aP, 1), stIntTuple_get(aP, 2)); - aW->rWeight = makeAlignmentWeight(columns, stIntTuple_get(aP, 0), stIntTuple_get(aP, 3), stIntTuple_get(aP, 4)); + double weight = ((double) stIntTuple_get(aP, 0)) / PAIR_ALIGNMENT_PROB_1 + + weightJitter(stIntTuple_get(aP, 1), stIntTuple_get(aP, 2), stIntTuple_get(aP, 3), stIntTuple_get(aP, 4)); + AlignmentWeight *aW = makeAlignmentWeight(columns, weight, stIntTuple_get(aP, 1), stIntTuple_get(aP, 2)); + aW->rWeight = makeAlignmentWeight(columns, weight, stIntTuple_get(aP, 3), stIntTuple_get(aP, 4)); aW->rWeight->rWeight = aW; insertWeight(aW, alignmentWeightAdjLists); insertWeight(aW->rWeight, alignmentWeightAdjLists); @@ -184,7 +227,7 @@ stSortedSet *makeOrderedSetOfAlignmentWeights(stHash *alignmentWeightAdjLists) { AlignmentWeight *aW; while ((aW = stSortedSet_getNext(aWIt)) != NULL) { assert(aW->column != aW->rWeight->column); - if (aW->column < aW->rWeight->column) { + if (alignmentWeight_isCanonical(aW)) { stSortedSet_insert(alignmentWeightsOrderedByWeight, aW); } } @@ -199,7 +242,7 @@ static void removeAlignmentFromSortedAlignmentWeights(AlignmentWeight *aW, stSor * Removes the weight from the ordered set, if one is being kept (the progressive merge does not need it). */ if (alignmentWeightsOrderedByWeight != NULL) { - stSortedSet_remove(alignmentWeightsOrderedByWeight, aW->column < aW->rWeight->column ? aW : aW->rWeight); + stSortedSet_remove(alignmentWeightsOrderedByWeight, alignmentWeight_isCanonical(aW) ? aW : aW->rWeight); } } @@ -208,7 +251,7 @@ static void insertAlignmentIntoSortedAlignmentWeights(AlignmentWeight *aW, stSor * Adds weight to the ordered set, if one is being kept. */ if (alignmentWeightsOrderedByWeight != NULL) { - stSortedSet_insert(alignmentWeightsOrderedByWeight, aW->column < aW->rWeight->column ? aW : aW->rWeight); + stSortedSet_insert(alignmentWeightsOrderedByWeight, alignmentWeight_isCanonical(aW) ? aW : aW->rWeight); } } @@ -608,7 +651,7 @@ static ProgressiveMsa *progressiveMsa_construct(stList *seqFrags, stList *multip m->isHead[m->seqOffsets[seq] + pos] = 1; } } - // The weights, made in the order, and with the random tie breaking, that makeAlignmentWeightAdjacencyLists uses + // The weights, made in the order, and with the tie-breaking jitter, that makeAlignmentWeightAdjacencyLists uses int64_t pairNumber = stList_length(multipleAlignedPairs); m->weights = st_calloc(n + 1, sizeof(WeightList)); for (int64_t i = 0; i < pairNumber; i++) { // each list's size, so each can have its own slice of one block @@ -630,12 +673,14 @@ static ProgressiveMsa *progressiveMsa_construct(stList *seqFrags, stList *multip int64_t id1 = m->seqOffsets[stIntTuple_get(aP, 1)] + stIntTuple_get(aP, 2); int64_t id2 = m->seqOffsets[stIntTuple_get(aP, 3)] + stIntTuple_get(aP, 4); AlignmentWeight *aW = &m->weightPool[2 * i], *rW = &m->weightPool[2 * i + 1]; + double weight = ((double) score) / PAIR_ALIGNMENT_PROB_1 + + weightJitter(stIntTuple_get(aP, 1), stIntTuple_get(aP, 2), stIntTuple_get(aP, 3), stIntTuple_get(aP, 4)); aW->column = m->columns[id1]; aW->numberOfWeights = 1; - aW->avgWeight = ((double) score) / PAIR_ALIGNMENT_PROB_1 + st_random() * 0.00001; // as makeAlignmentWeight + aW->avgWeight = weight; rW->column = m->columns[id2]; rW->numberOfWeights = 1; - rW->avgWeight = ((double) score) / PAIR_ALIGNMENT_PROB_1 + st_random() * 0.00001; + rW->avgWeight = weight; aW->rWeight = rW; rW->rWeight = aW; // A weight lives in the list of the column at the other end of it, as insertWeight puts it @@ -1271,7 +1316,9 @@ int64_t getNextBestPair(int64_t seq1, int64_t *distanceCounts, int64_t seqNo, for (int64_t seq2 = 0; seq2 < seqNo; seq2++) { if (seq1 != seq2) { double gain = distances[seq2] - subsPerSite(seq1, seq2, distanceCounts, seqNo); - if (gain > maxGain || (gain == maxGain && st_random() > 0.5) /*lame attempt to reduce always picking the same seq*/) { + // ties go by a hash of the pair, which spreads the choice over the sequences (the point of the coin + // flip this replaces) without making it differ from one run to the next + if (gain > maxGain || (gain == maxGain && pairHash(seq1, 0, seq2, 0) > pairHash(seq1, 0, maxGainSeq, 0))) { stIntTuple *pairToAlign = makePairToAlign(seq1, seq2); if (stSortedSet_search(chosenPairsOfSequencesToAlign, pairToAlign) == NULL) { //So that any pair is unique maxGain = gain; diff --git a/tests/multipleAlignerTest.c b/tests/multipleAlignerTest.c index 54e04e7..01d05ed 100644 --- a/tests/multipleAlignerTest.c +++ b/tests/multipleAlignerTest.c @@ -250,25 +250,47 @@ static StateMachine *choosePairStateMachine(int64_t seqX, int64_t seqY, void *ex return choice->sM; } -static bool sameAlignedPairs(stList *pairs1, stList *pairs2) { - // as sets: the order the pairs come out in follows hashes of pointers, so differs from run to run - if (stList_length(pairs1) != stList_length(pairs2)) { +static bool sameTuples(stList *tuples1, stList *tuples2) { + // in order, as well as in content + if (stList_length(tuples1) != stList_length(tuples2)) { + return 0; + } + for (int64_t i = 0; i < stList_length(tuples1); i++) { + if (!stIntTuple_equalsFn(stList_get(tuples1, i), stList_get(tuples2, i))) { + return 0; + } + } + return 1; +} + +static bool sameColumns(stSet *columns1, stSet *columns2) { + // each column the same positions, in the same order + if (stSet_size(columns1) != stSet_size(columns2)) { return 0; } - stList *sorted1 = stList_copy(pairs1, NULL), *sorted2 = stList_copy(pairs2, NULL); - stList_sort(sorted1, (int (*)(const void *, const void *))stIntTuple_cmpFn); - stList_sort(sorted2, (int (*)(const void *, const void *))stIntTuple_cmpFn); bool same = 1; - for (int64_t i = 0; i < stList_length(sorted1) && same; i++) { - same = stIntTuple_equalsFn(stList_get(sorted1, i), stList_get(sorted2, i)); + stSetIterator *it = stSet_getIterator(columns1); + Column *c; + while (same && (c = stSet_getNext(it)) != NULL) { + Column *c1 = c, *c2 = stSet_search(columns2, c); + while (c1 != NULL && c2 != NULL && c1->seqName == c2->seqName && c1->position == c2->position) { + c1 = c1->nColumn; + c2 = c2->nColumn; + } + same = c1 == NULL && c2 == NULL; } - stList_destruct(sorted1); - stList_destruct(sorted2); + stSet_destructIterator(it); return same; } +static bool sameMultipleAlignment(MultipleAlignment *mA, MultipleAlignment *mA2) { + return sameTuples(mA->alignedPairs, mA2->alignedPairs) && + sameTuples(mA->chosenPairwiseAlignments, mA2->chosenPairwiseAlignments) && + sameColumns(mA->columns, mA2->columns); +} + static void test_makeAlignmentWithPairStateMachines(CuTest *testCase) { - int64_t differed = 0, repeatable = 0; + int64_t differed = 0; for (int64_t test = 0; test < 50; test++) { setup(); // spanning trees from 1 to 4 over 2 to 9 sequences reach both the all-pairs and the spanning-tree path @@ -277,13 +299,6 @@ static void test_makeAlignmentWithPairStateMachines(CuTest *testCase) { int64_t spanningTrees = st_randomInt(1, 5); bool progressive = st_random() > 0.5; MultipleAlignment *mA = makeAlignment(stateMachine, randomSeqFrags, spanningTrees, 10000000, progressive, 0.5, pabp); - // makeAlignment can give different alignments of the same sequences from one call to the next (seen - // without progressive merging, on the spanning-tree path), so equality below is only asked for where it - // agrees with itself - MultipleAlignment *again = makeAlignment(stateMachine, randomSeqFrags, spanningTrees, 10000000, progressive, 0.5, pabp); - bool isRepeatable = sameAlignedPairs(mA->alignedPairs, again->alignedPairs); - repeatable += isRepeatable; - multipleAlignment_destruct(again); // a chooser that returns NULL, or the default machine itself, changes nothing, and is asked once per // pairwise alignment about a real pair of sequences @@ -291,7 +306,7 @@ static void test_makeAlignmentWithPairStateMachines(CuTest *testCase) { PairChoice choice = { i == 0 ? NULL : stateMachine, seqNo, 0, 0 }; MultipleAlignment *mA2 = makeAlignmentWithPairStateMachines(stateMachine, choosePairStateMachine, &choice, randomSeqFrags, spanningTrees, 10000000, progressive, 0.5, pabp); - CuAssertTrue(testCase, !isRepeatable || sameAlignedPairs(mA->alignedPairs, mA2->alignedPairs)); + CuAssertTrue(testCase, sameMultipleAlignment(mA, mA2)); CuAssertIntEquals(testCase, stList_length(mA2->chosenPairwiseAlignments), choice.calls); CuAssertTrue(testCase, !choice.badPair); multipleAlignment_destruct(mA2); @@ -306,7 +321,7 @@ static void test_makeAlignmentWithPairStateMachines(CuTest *testCase) { randomSeqFrags, spanningTrees, 10000000, progressive, 0.5, pabp); checkAlignment(testCase, randomSeqFrags, mA3->alignedPairs); CuAssertIntEquals(testCase, stList_length(mA3->chosenPairwiseAlignments), choice.calls); - differed += isRepeatable && !sameAlignedPairs(mA->alignedPairs, mA3->alignedPairs); + differed += !sameTuples(mA->alignedPairs, mA3->alignedPairs); multipleAlignment_destruct(mA3); stateMachine_destruct(other); hmm_destruct(hmm); @@ -315,10 +330,39 @@ static void test_makeAlignmentWithPairStateMachines(CuTest *testCase) { stList_destruct(randomSeqFrags); teardown(); } - CuAssertTrue(testCase, repeatable >= 40); CuAssertTrue(testCase, differed > 0); } +/* + * The same sequences give the same alignment, whatever ran before: neither the state of the random number + * generator nor where memory is allocated makes a difference. + */ +static void test_multipleAlignerRepeatable(CuTest *testCase) { + for (int64_t test = 0; test < 100; test++) { + setup(); + // as in test_makeAlignmentWithPairStateMachines, and sometimes with few enough pairs considered for the + // distances that which columns give them matters + stList *randomSeqFrags = getRandomSeqFrags(st_randomInt(2, 10), st_randomInt(20, 100)); + int64_t spanningTrees = st_randomInt(1, 5); + int64_t maxPairsToConsider = st_random() > 0.5 ? 10000000 : st_randomInt(1, 200); + bool progressive = st_random() > 0.5; + MultipleAlignment *mA = makeAlignment(stateMachine, randomSeqFrags, spanningTrees, maxPairsToConsider, progressive, 0.5, pabp); + // move the generator on, and hold some memory so that the second alignment's allocations land elsewhere + stList *junk = stList_construct3(0, free); + for (int64_t i = st_randomInt(1, 1000); i > 0; i--) { + stList_append(junk, st_malloc(st_randomInt(1, 200))); + } + MultipleAlignment *mA2 = makeAlignment(stateMachine, randomSeqFrags, spanningTrees, maxPairsToConsider, progressive, 0.5, pabp); + checkAlignment(testCase, randomSeqFrags, mA2->alignedPairs); + CuAssertTrue(testCase, sameMultipleAlignment(mA, mA2)); + stList_destruct(junk); + multipleAlignment_destruct(mA2); + multipleAlignment_destruct(mA); + stList_destruct(randomSeqFrags); + teardown(); + } +} + CuSuite* multipleAlignerTestSuite(void) { CuSuite* suite = CuSuiteNew(); SUITE_ADD_TEST(suite, test_getDistanceMatrix); @@ -330,6 +374,7 @@ CuSuite* multipleAlignerTestSuite(void) { SUITE_ADD_TEST(suite, test_multipleAlignerAllPairsRandom); SUITE_ADD_TEST(suite, test_multipleAlignerRandom); SUITE_ADD_TEST(suite, test_makeAlignmentWithPairStateMachines); + SUITE_ADD_TEST(suite, test_multipleAlignerRepeatable); return suite; }