diff --git a/impl/multipleAligner.c b/impl/multipleAligner.c index 7baacb7..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 @@ -1025,14 +1070,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 +1100,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 +1109,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 +1138,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); @@ -1244,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; @@ -1262,12 +1336,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 +1363,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 +1388,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..01d05ed 100644 --- a/tests/multipleAlignerTest.c +++ b/tests/multipleAlignerTest.c @@ -231,6 +231,138 @@ 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 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; + } + bool same = 1; + 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; + } + 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; + 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); + + // 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, sameMultipleAlignment(mA, mA2)); + 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 += !sameTuples(mA->alignedPairs, mA3->alignedPairs); + multipleAlignment_destruct(mA3); + stateMachine_destruct(other); + hmm_destruct(hmm); + + multipleAlignment_destruct(mA); + stList_destruct(randomSeqFrags); + teardown(); + } + 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); @@ -241,6 +373,8 @@ 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); + SUITE_ADD_TEST(suite, test_multipleAlignerRepeatable); return suite; }