From 6f74aac1971a96a364880cf40c2a84cee0bac279 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Tue, 25 Aug 2026 07:42:43 -0400 Subject: [PATCH 1/5] Make cactus_consolidated's output reproducible The same input and command produced a different alignment on every run, at any thread count. Four causes lived here; a fifth was in sonLib, whose submodule pointer this bumps to the two commits that fix it. Names came off one shared counter inside the OpenMP loops in bar and in reference building, so which object got which Name depended on how the threads interleaved -- and Names are written into the .c2h. A loop now reserves one block of names up front and hands each iteration a private sub-interval chosen from the iteration index, so naming depends on a flower's position in the list rather than on scheduling. The intervals are sized from an upper bound on what an iteration can build; overrunning one falls back to the shared counter and says so, rather than failing. st_random drew from one process-global sequence. Both places cactus breaks ties randomly -- the reference ordering's simulated annealing, and ancestral base calling -- do so inside parallel loops, so the draws were a different slice of that sequence each run. With sonLib's generator now thread-local, each iteration seeds from its own flower, and the root flower is seeded before the serial call that follows the loop. Two comparators could not give a stable order: sortByEvent ranked segments by the address of their Event, and several others left ties to qsort. They now order by name, with names as tie-breaks, so each is a total order. getFirstSegmentMatchingEvent moves with sortByEvent, since it binary-searches the list that sort produces. Thread components were attached to the dead-end component in the order an stSortedSet of pointers yielded them, which decided which threads got attached to the root. They are now sorted by their lowest thread name, and each component's threads by name. This changes cactus's alignments, not just their serialisation: different tie-breaks and a different random sequence give slightly different -- not worse -- output than 3.3.x. Stored baselines will need refreshing. What this establishes is same input + same command + same binary. It is not bit-identity across platforms: several comparators still hand ties to qsort_r, whose behaviour on equal elements is unspecified and changed in glibc 2.37. Verified: cactus_consolidated byte-identical across runs at 1, 2, 4, 8, 16 and 32 threads on two jobs; the progressive and pangenome pipelines identical in every cactus-produced output; alignment accuracy against the evolver truth maf unchanged, inside the old build's own run-to-run spread; peak RSS and wall time unchanged on a 266 Mb job, where hash entries account for 1.3% of peak and this adds 0.43% of it. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_017E9w7KrbFEfqNRtkam4u1F --- api/impl/cactusDisk.c | 58 +++++++++++++++++++++++++++++++++ api/inc/cactusDisk.h | 30 +++++++++++++++++ bar/impl/bar.c | 35 ++++++++++++++++++++ bar/impl/poaBarAligner.c | 7 +++- caf/impl/pinchToCactus.c | 38 +++++++++++++++++---- pipeline/cactus_consolidated.c | 17 ++++++++-- reference/impl/blockMLString.c | 13 ++++++-- reference/impl/buildReference.c | 31 +++++++++++++++++- submodules/sonLib | 2 +- 9 files changed, 218 insertions(+), 13 deletions(-) diff --git a/api/impl/cactusDisk.c b/api/impl/cactusDisk.c index 69f06f8f4..39aa5d089 100644 --- a/api/impl/cactusDisk.c +++ b/api/impl/cactusDisk.c @@ -207,7 +207,35 @@ void cactusDisk_setEventTree(CactusDisk *cactusDisk, EventTree *eventTree) { * Function to get unique ID. */ +/* + * The interval a parallel loop has handed to this thread, if any. See + * cactusDisk_pushNameInterval in cactusDisk.h. + */ +#if defined(__GNUC__) || defined(__clang__) +#define CACTUS_THREAD_LOCAL __thread +#else +#define CACTUS_THREAD_LOCAL +#endif + +static CACTUS_THREAD_LOCAL Name intervalNext = 0; +static CACTUS_THREAD_LOCAL int64_t intervalRemaining = 0; +static CACTUS_THREAD_LOCAL bool intervalActive = 0; + int64_t cactusDisk_getUniqueIDInterval(CactusDisk *cactusDisk, int64_t intervalSize) { + if (intervalRemaining >= intervalSize) { // Thread has names of its own left, so no lock and no shared state + Name n = intervalNext; + intervalNext += intervalSize; + intervalRemaining -= intervalSize; + return n; + } + if (intervalActive) { // Ran out: correct, but the names stop being reproducible from here on + static bool warned = 0; + if (!warned) { + warned = 1; + st_logCritical("Warning: ran out of reserved names in a parallel region, so names, and any " + "output containing them, will not be identical from run to run\n"); + } + } #if defined(_OPENMP) omp_set_lock(&(cactusDisk->writelock)); #endif @@ -223,6 +251,36 @@ int64_t cactusDisk_getUniqueID(CactusDisk *cactusDisk) { return cactusDisk_getUniqueIDInterval(cactusDisk, 1); } +Name cactusDisk_reserveNames(CactusDisk *cactusDisk, int64_t nameNumber) { + assert(nameNumber >= 0); + assert(!intervalActive); // Reserving is done outside the parallel loop that uses the names +#if defined(_OPENMP) + omp_set_lock(&(cactusDisk->writelock)); +#endif + Name n = cactusDisk->currentName; + cactusDisk->currentName += nameNumber; + assert(cactusDisk->currentName > n || nameNumber == 0); // Overflow of the name space +#if defined(_OPENMP) + omp_unset_lock(&(cactusDisk->writelock)); +#endif + return n; +} + +void cactusDisk_pushNameInterval(CactusDisk *cactusDisk, Name start, int64_t nameNumber) { + assert(nameNumber >= 0); + assert(!intervalActive); // Intervals do not nest + intervalNext = start; + intervalRemaining = nameNumber; + intervalActive = 1; +} + +void cactusDisk_popNameInterval(CactusDisk *cactusDisk) { + assert(intervalActive); + intervalNext = 0; + intervalRemaining = 0; + intervalActive = 0; +} + EventTree *cactusDisk_getEventTree(CactusDisk *cactusDisk) { return cactusDisk->eventTree; } diff --git a/api/inc/cactusDisk.h b/api/inc/cactusDisk.h index fdbeb0e0c..bc72b0ce8 100644 --- a/api/inc/cactusDisk.h +++ b/api/inc/cactusDisk.h @@ -41,6 +41,36 @@ int64_t cactusDisk_getUniqueID(CactusDisk *cactusDisk); */ int64_t cactusDisk_getUniqueIDInterval(CactusDisk *cactusDisk, int64_t intervalSize); +/* + * Deterministic naming inside parallel loops. + * + * Names normally come off one shared counter, so an object's name depends on + * how the threads happened to interleave and differs from run to run. To avoid + * that, a parallel loop hands each iteration its own interval of names, chosen + * from the iteration index rather than from the thread: + * + * Name base = cactusDisk_reserveNames(cactusDisk, total); // before the loop + * #pragma omp parallel for + * for (i ...) { + * cactusDisk_pushNameInterval(cactusDisk, base + offset[i], size[i]); + * ... construct objects ... + * cactusDisk_popNameInterval(cactusDisk); + * } + * + * The intervals must not overlap, and the caller has to size them: sizing is a + * per-loop matter, since only the caller knows what an iteration can build. + * Anything an iteration allocates beyond its interval falls back to the shared + * counter -- correct, but not reproducible, so it is logged as a warning. + * + * The interval is held in thread-local state, so it applies to the thread that + * pushed it and only until it pops. + */ +Name cactusDisk_reserveNames(CactusDisk *cactusDisk, int64_t nameNumber); + +void cactusDisk_pushNameInterval(CactusDisk *cactusDisk, Name start, int64_t nameNumber); + +void cactusDisk_popNameInterval(CactusDisk *cactusDisk); + /* * Gets a flower the cactusDisk contains. If the flower is not in memory it will be loaded. If not in memory or on disk, returns NULL. */ diff --git a/bar/impl/bar.c b/bar/impl/bar.c index cece3c849..bbc87754e 100644 --- a/bar/impl/bar.c +++ b/bar/impl/bar.c @@ -82,6 +82,35 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * st_errAbort("We have precomputed alignments but %" PRIi64 " flowers to align.\n", stList_length(flowers)); } + /* + * Give each flower its own interval of names, so that the names of the blocks, + * segments and so on that get built below depend on the flower's position in + * the list and not on the order the threads happen to reach it. + * + * The bound: every base in the flower can end up in a pinch segment of its own, + * each segment costs three names (itself and its two caps), and there is at most + * one block per segment, at three names again. Doubling that covers the ends, + * groups, chains and nested flowers, which are all bounded by the number of + * blocks. A flower's bases are its unaligned length plus two per adjacency. + */ + int64_t flowerNumber = stList_length(flowers); + int64_t *nameOffsets = st_malloc(sizeof(int64_t) * (flowerNumber + 1)); +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1) +#endif + for (int64_t j = 0; j < flowerNumber; j++) { + Flower *flower = stList_get(flowers, j); + nameOffsets[j] = 12 * (flower_getTotalBaseLength(flower) + flower_getCapNumber(flower)) + 4096; + } + int64_t nameTotal = 0; // Turn the sizes into offsets + for (int64_t j = 0; j < flowerNumber; j++) { + int64_t size = nameOffsets[j]; + nameOffsets[j] = nameTotal; + nameTotal += size; + } + nameOffsets[flowerNumber] = nameTotal; + Name nameBase = cactusDisk_reserveNames(cactusDisk, nameTotal); + #if defined(_OPENMP) // Enable nested parallelism for large flowers with many ends // Level 0: serial (main thread) @@ -94,6 +123,9 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * for (int64_t j = 0; jminimumIngroupDegree = cactusParams_get_int(params, 2, "bar", "minimumIngroupDegree"); @@ -160,8 +192,11 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * } free(fa); + cactusDisk_popNameInterval(cactusDisk); + st_logDebug("Finished filling in the alignments for the flower\n"); } + free(nameOffsets); ////////////////////////////////////////////// //Clean up diff --git a/bar/impl/poaBarAligner.c b/bar/impl/poaBarAligner.c index 843ad3476..d4f162d9b 100644 --- a/bar/impl/poaBarAligner.c +++ b/bar/impl/poaBarAligner.c @@ -1055,7 +1055,12 @@ int caps_comp_by_adjacency_length(const void *a, const void *b) { int length1, length2; get_adjacency_string((Cap *)a, &length1, 0); get_adjacency_string((Cap *)b, &length2, 0); - return length1 > length2 ? -1 : (length1 < length2 ? 1 : 0); // sort in descending order of length + if (length1 != length2) { + return length1 > length2 ? -1 : 1; // sort in descending order of length + } + // Equal lengths are common, and the row order decides the alignment, so + // separate them by name instead of leaving it to whatever qsort does + return cactusMisc_nameCompare(cap_getName((Cap *)a), cap_getName((Cap *)b)); } /* diff --git a/caf/impl/pinchToCactus.c b/caf/impl/pinchToCactus.c index 38dcbb7b8..408e0b7fb 100644 --- a/caf/impl/pinchToCactus.c +++ b/caf/impl/pinchToCactus.c @@ -105,10 +105,26 @@ static void attachThreadToDeadEndComponent(stPinchThread *thread, stList *deadEn } } +static int comparePinchThreadsByName(const void *a, const void *b) { + int64_t i = stPinchThread_getName((stPinchThread *) a); + int64_t j = stPinchThread_getName((stPinchThread *) b); + return i > j ? 1 : (i < j ? -1 : 0); +} + int comparePinchThreadsByLength(const void *a, const void *b) { int64_t i = stPinchThread_getLength((stPinchThread *) a); int64_t j = stPinchThread_getLength((stPinchThread *) b); - return i > j ? 1 : (i < j ? -1 : 0); + if (i != j) { + return i > j ? 1 : -1; + } + return comparePinchThreadsByName(a, b); // Break ties by name: which of two equally long + // threads gets attached first has to be decided the same way on every run +} + +static int compareThreadComponentsByFirstThread(const void *a, const void *b) { + // The components are disjoint and thread names are unique, so comparing their + // lowest-named thread is a total order over the components + return comparePinchThreadsByName(stList_get((stList *) a, 0), stList_get((stList *) b, 0)); } static void attachThreadComponentToDeadEndComponent(stList *threadComponent, stList *deadEndComponent, @@ -220,13 +236,23 @@ static void stCaf_attachUnattachedThreadComponents(Flower *flower, stPinchThread } stSortedSet *threadComponents = stPinchThreadSet_getThreadComponents(threadSet); assert(stSortedSet_size(threadComponents) > 0); - stSortedSetIterator *threadIt = stSortedSet_getIterator(threadComponents); - stList *threadComponent; - while ((threadComponent = stSortedSet_getNext(threadIt)) != NULL) { - attachThreadComponentToDeadEndComponent(threadComponent, deadEndComponent, pinchEndsToAdjacencyComponents, markEndsAttached, + /* + * The set holds the components in address order, and each component holds its + * threads in whatever order the union-find produced them. Which threads get + * attached to the root depends on the order they are visited in, so put both + * into name order first, which is the same on every run. + */ + stList *threadComponents2 = stSortedSet_getList(threadComponents); + for (int64_t i = 0; i < stList_length(threadComponents2); i++) { + stList_sort(stList_get(threadComponents2, i), comparePinchThreadsByName); + } + stList_sort(threadComponents2, compareThreadComponentsByFirstThread); + for (int64_t i = 0; i < stList_length(threadComponents2); i++) { + attachThreadComponentToDeadEndComponent(stList_get(threadComponents2, i), deadEndComponent, + pinchEndsToAdjacencyComponents, markEndsAttached, minLengthForChromosome, proportionOfUnalignedBasesForNewChromosome, flower); } - stSortedSet_destructIterator(threadIt); + stList_destruct(threadComponents2); stSortedSet_destruct(threadComponents); } diff --git a/pipeline/cactus_consolidated.c b/pipeline/cactus_consolidated.c index b0e92fe57..7a60fe32c 100644 --- a/pipeline/cactus_consolidated.c +++ b/pipeline/cactus_consolidated.c @@ -126,6 +126,9 @@ static RecordHolder *doBottomUpTraversal(stList *flowerLayers, #pragma omp parallel for schedule(dynamic) #endif for (int64_t j = 0; j < stList_length(flowers); j++) { + // Ancestral base calling breaks ties randomly: seed from the flower so the + // draws are the flower's own and not a slice of one sequence shared by the threads + st_randomSeed(flower_getName(stList_get(flowers, j))); stList_set(recordHoldersForFlowers, j, getMergedRecordHolders(recordHolders, stList_get(flowers, j))); bottomUpFn(stList_get(flowers, j), stList_get(recordHoldersForFlowers, j), extraArgs); } @@ -176,17 +179,24 @@ stHash *compute_flower_length_hash(stList *flowers) { return flower_to_length; } +// Flowers of equal size are ordered by name. The sort decides which flower gets +// which interval of names, so leaving equal-sized flowers to be separated by +// whatever qsort does with them would tie the output to the C library's sort. +static int flower_nameCmpFn(const void *a, const void *b) { + return cactusMisc_nameCompare(flower_getName((Flower *)a), flower_getName((Flower *)b)); +} + int flower_lengthCmpFn(const void *a, const void *b, void *flower_to_length_hash) { // Sort by hashed length value of the flowers int64_t i = (int64_t)stHash_search((stHash*)flower_to_length_hash, (void*)a); int64_t j = (int64_t)stHash_search((stHash*)flower_to_length_hash, (void*)b); - return i < j ? 1 : (i > j ? -1 : 0); // Sort in descending order + return i < j ? 1 : (i > j ? -1 : flower_nameCmpFn(a, b)); // Sort in descending order } int flower_sizeCmpFn(const void *a, const void *b) { // Sort by number of caps the flowers contains int64_t i = flower_getCapNumber((Flower *)a), j = flower_getCapNumber((Flower *)b); - return i < j ? 1 : (i > j ? -1 : 0); // Sort in descending order + return i < j ? 1 : (i > j ? -1 : flower_nameCmpFn(a, b)); // Sort in descending order } int main(int argc, char *argv[]) { @@ -486,6 +496,9 @@ int main(int argc, char *argv[]) { // Bottom-up reference coordinates phase RecordHolder *rh = doBottomUpTraversal(flowerLayers, callBottomUp, (void *)referenceEventName); + // The traversal above left this thread's generator wherever the flowers it + // happened to be handed took it, so seed the root flower like any other + st_randomSeed(flower_getName(flower)); bottomUpNoDb(flower, rh, referenceEventName, 1, generateJukesCantorMatrix); assert(recordHolder_size(rh) == 0); recordHolder_destruct(rh); diff --git a/reference/impl/blockMLString.c b/reference/impl/blockMLString.c index 77657d0e5..9d82e02a2 100644 --- a/reference/impl/blockMLString.c +++ b/reference/impl/blockMLString.c @@ -225,7 +225,9 @@ static void multiply(double *baseProbs1, double *baseProbs2, int64_t blockLength static int getFirstSegmentMatchingEvent(const void *a, const void *b) { Event *e1 = (Event *)a, *e2 = segment_getEvent((Segment *)b); assert(e1 != NULL && e2 != NULL); - return e1 < e2 ? -1 : (e1 > e2 ? 1 : 0); + // Must order events the same way sortByEvent does, as this searches the list + // it produced + return cactusMisc_nameCompare(event_getName(e1), event_getName(e2)); } static double *computeBaseProbs(stTree *tree, stList *eventSortedSegments, int64_t blockLength) { @@ -315,7 +317,14 @@ void maskAncestralRepeatBases(Block *block, stList *segments, char *mlString) { static int sortByEvent(const void *a, const void *b) { Event *e1 = segment_getEvent((Segment *)a), *e2 = segment_getEvent((Segment *)b); assert(e1 != NULL && e2 != NULL); - return e1 < e2 ? -1 : (e1 > e2 ? 1 : 0); + // By name rather than by address: the addresses order the events differently + // from one run to the next, and the order the segments are summed in decides + // the ancestral bases. Segment name breaks ties, so the order is total. + int i = cactusMisc_nameCompare(event_getName(e1), event_getName(e2)); + if (i != 0) { + return i; + } + return cactusMisc_nameCompare(segment_getName((Segment *)a), segment_getName((Segment *)b)); } static stList *segmentsSortedByEvent(Block *block) { diff --git a/reference/impl/buildReference.c b/reference/impl/buildReference.c index 1345d92d8..b6263ec24 100644 --- a/reference/impl/buildReference.c +++ b/reference/impl/buildReference.c @@ -1320,13 +1320,42 @@ void cactus_make_reference(stList *flowers, char *referenceEventString, double (*temperatureFn)(double) = useSimulatedAnnealing ? exponentiallyDecreasingTemperatureFn : constantTemperatureFn; + /* + * Give each flower its own interval of names, so that the reference caps and + * segments built below are named from the flower's position in the list rather + * than from the order the threads happen to reach it. Per flower this adds at + * most a cap per end, a segment (three names) per block, and a block plus an end + * per group for the scaffold gaps. + */ + int64_t flowerNumber = stList_length(flowers); + int64_t *nameOffsets = st_malloc(sizeof(int64_t) * (flowerNumber + 1)); +#pragma omp parallel for schedule(dynamic, 1) + for (int64_t i = 0; i < flowerNumber; i++) { + Flower *flower = stList_get(flowers, i); + nameOffsets[i] = 8 * (flower_getEndNumber(flower) + flower_getBlockNumber(flower) + + flower_getGroupNumber(flower) + flower_getCapNumber(flower)) + 4096; + } + int64_t nameTotal = 0; // Turn the sizes into offsets + for (int64_t i = 0; i < flowerNumber; i++) { + int64_t size = nameOffsets[i]; + nameOffsets[i] = nameTotal; + nameTotal += size; + } + nameOffsets[flowerNumber] = nameTotal; + Name nameBase = cactusDisk_reserveNames(cactusDisk, nameTotal); + #pragma omp parallel for schedule(dynamic, 1) - for(int64_t i=0; i Date: Tue, 25 Aug 2026 07:50:41 -0400 Subject: [PATCH 2/5] Pin hal to the external-array buffer zeroing fix hal wrote uninitialised heap bytes into every file: Hdf5ExternalArray's buffer comes from new char[], and not every byte of it is set before the buffer is flushed. Two cactus runs producing the same alignment therefore wrote .hal files that differed in bytes nothing ever reads, which meant h5diff and checksums could not be used to compare them. This is the last piece of making cactus's output reproducible; the rest is in the previous commit. Only that fix comes in: the merge commit in between changes no files, so 77075929..7b418b56 is exactly the 12 added lines. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_017E9w7KrbFEfqNRtkam4u1F --- submodules/hal | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/submodules/hal b/submodules/hal index 770759293..7b418b561 160000 --- a/submodules/hal +++ b/submodules/hal @@ -1 +1 @@ -Subproject commit 770759293decff32fa036696d10590cac55f0724 +Subproject commit 7b418b56102dc92c21b8dd4aac1df56e088ba26d From ee4251ac9b600aab317c617f38a9410f41ac25ae Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Tue, 25 Aug 2026 08:16:28 -0400 Subject: [PATCH 3/5] Pin paffy to its sonLib determinism bump The build forces paffy's nested sonLib to whatever SHA cactus's own sonLib is at, so the two have to agree or a fresh clone checks out a commit paffy has never fetched. 0d26e841..67ca7364 changes one file: paffy's sonLib pointer, now cbb2285 like ours. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_017E9w7KrbFEfqNRtkam4u1F --- submodules/paffy | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/submodules/paffy b/submodules/paffy index 0d26e841f..67ca7364f 160000 --- a/submodules/paffy +++ b/submodules/paffy @@ -1 +1 @@ -Subproject commit 0d26e841f6ed493ce703fe331c369af3615119a1 +Subproject commit 67ca7364fed880e395a667da42f2347f991ce0d6 From baa72d159d01885f9e8c5c2c3626624fab184ed4 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Tue, 25 Aug 2026 08:20:11 -0400 Subject: [PATCH 4/5] Update the pinned taffy to its sonLib determinism bump downloadMafTools builds taffy from a pinned commit and initialises its submodules recursively, so taffy brings its own sonLib. Move the pin so that sonLib is the same one cactus, hal and paffy now use. 0b6e42fd..2f964e8c is a fast-forward whose whole file diff is taffy's sonLib pointer; the merge commit in between changes nothing. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_017E9w7KrbFEfqNRtkam4u1F --- build-tools/downloadMafTools | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build-tools/downloadMafTools b/build-tools/downloadMafTools index 3eacacd55..e9f152509 100755 --- a/build-tools/downloadMafTools +++ b/build-tools/downloadMafTools @@ -63,7 +63,7 @@ export HTSLIB_LIBS="$(pwd)/libhts.a -lbz2 -ldeflate -lm -lpthread -lz -llzma -pt cd ${mafBuildDir} git clone https://github.com/ComparativeGenomicsToolkit/taffy.git cd taffy -git checkout 0b6e42fd8cf3267ae84fe84ea855f7883c2ee75d +git checkout 2f964e8c3e8f61a134b7c4dd93ba1a5056116f65 git submodule update --init --recursive export HALDIR=${CWD}/submodules/hal LIBS="${jemallocLib} -llz4 -lzstd" make -j ${numcpu} From a3e7176a9db4a7d10085ff83b2947606c0b6f5b9 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Wed, 26 Aug 2026 09:50:16 -0400 Subject: [PATCH 5/5] Stop taking workflow ordering from os.listdir Every place the workflow collected files from a directory did so in os.listdir order, which is whatever the filesystem hands back. The one that matters is the lastz chunking. faffy names its chunks 0.fa, 1.fa, ... and cactus took them in listdir order, so chunk i -- and hence the lastz job named _i -- was an arbitrary piece of the genome rather than the i'th piece of it. faffy already prints each chunk's path to stdout as it completes it, in the order the chunks tile the input, so take the list from there. Sorting the names would not have worked: they sort lexicographically, 0, 1, 10, 11, 2. The remaining sites take sorted(). Most concatenate a genome supplied as a directory of fastas, where the order sets sequence order in the merged file. This is not a run-to-run fix. listdir is stable for a given directory on a given filesystem, so a machine's own results were already reproducible; the exposure is across filesystems. The same 53 chunk names list in three different orders on tmpfs and on two ext4 volumes, because ext4 hashes directory entries against a seed chosen at mkfs time -- so the same job run on another disk, or in a container, chunked the genome differently. Behaviour is otherwise unchanged: a 120-lastz-job chunked blast produces a byte-identical paf before and after. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_017E9w7KrbFEfqNRtkam4u1F --- src/cactus/blast/cactus_blast.py | 2 +- src/cactus/paf/local_alignment.py | 10 +++++++--- src/cactus/preprocessor/cactus_preprocessor.py | 6 ++++-- src/cactus/progressive/cactus_prepare.py | 2 +- src/cactus/progressive/cactus_progressive.py | 2 +- src/cactus/refmap/cactus_graphmap.py | 2 +- src/cactus/refmap/cactus_graphmap_split.py | 4 ++-- src/cactus/refmap/cactus_minigraph.py | 2 +- src/cactus/refmap/cactus_pangenome.py | 2 +- src/cactus/refmap/cactus_refmap.py | 2 +- src/cactus/setup/cactus_align.py | 2 +- 11 files changed, 21 insertions(+), 15 deletions(-) diff --git a/src/cactus/blast/cactus_blast.py b/src/cactus/blast/cactus_blast.py index 275ed4f23..1a1047174 100644 --- a/src/cactus/blast/cactus_blast.py +++ b/src/cactus/blast/cactus_blast.py @@ -150,7 +150,7 @@ def runCactusBlastOnly(options): if genome in event_set: if os.path.isdir(seq): tmpSeq = getTempFile() - catFiles([os.path.join(seq, subSeq) for subSeq in os.listdir(seq)], tmpSeq) + catFiles([os.path.join(seq, subSeq) for subSeq in sorted(os.listdir(seq))], tmpSeq) seq = tmpSeq seq = makeURL(seq) input_seq_id_map[genome] = toil.importFile(seq) diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py index 193df4bd3..e06db799e 100755 --- a/src/cactus/paf/local_alignment.py +++ b/src/cactus/paf/local_alignment.py @@ -382,9 +382,13 @@ def make_chunks(genome): '-o', params.find("blast").attrib["overlapSize"], '--dir', output_chunks_dir, job.fileStore.readGlobalFile(genome)] - cactus_call(parameters=fasta_chunk_cmd) - return [job.fileStore.writeGlobalFile(os.path.join(output_chunks_dir, chunk), cleanup=True) - for chunk in os.listdir(output_chunks_dir)] + # faffy prints the path of each chunk as it finishes it, in the order the chunks + # tile the input. Use that rather than os.listdir, whose order is arbitrary: + # it decides which piece of sequence is chunk i, and so what the lastz job named + # _i actually aligned. Sorting the names would not fix it either, since + # they are 0.fa, 1.fa, ... 10.fa and sort lexicographically. + chunk_paths = cactus_call(parameters=fasta_chunk_cmd, check_output=True).split() + return [job.fileStore.writeGlobalFile(chunk_path, cleanup=True) for chunk_path in chunk_paths] # Chunk each input genome chunks_a = make_chunks(genome_a) chunks_b = make_chunks(genome_b) diff --git a/src/cactus/preprocessor/cactus_preprocessor.py b/src/cactus/preprocessor/cactus_preprocessor.py index f97e08fb2..74e53cf65 100644 --- a/src/cactus/preprocessor/cactus_preprocessor.py +++ b/src/cactus/preprocessor/cactus_preprocessor.py @@ -597,8 +597,10 @@ def main(): except: pass assert os.path.isdir(inPath) == os.path.isdir(outPath) - inSeqPaths += [os.path.join(inPath, seqPath) for seqPath in os.listdir(inPath)] - outSeqPaths += [os.path.join(outPath, seqPath) for seqPath in os.listdir(inPath)] + # sorted so the two lists pair up the same way on every run, and so a + # directory of fastas is preprocessed in a stable order + inSeqPaths += [os.path.join(inPath, seqPath) for seqPath in sorted(os.listdir(inPath))] + outSeqPaths += [os.path.join(outPath, seqPath) for seqPath in sorted(os.listdir(inPath))] else: inSeqPaths += [inPath] outSeqPaths += [outPath] diff --git a/src/cactus/progressive/cactus_prepare.py b/src/cactus/progressive/cactus_prepare.py index 36d793816..f9f257798 100644 --- a/src/cactus/progressive/cactus_prepare.py +++ b/src/cactus/progressive/cactus_prepare.py @@ -1045,7 +1045,7 @@ def toil_call_blast(job, options, seq_file, mc_tree, og_map, event, cigar_name, # scrape the output files out of the workdir out_nameids = [] - for out_file in [f for f in os.listdir(work_dir) if os.path.isfile(os.path.join(work_dir, f))]: + for out_file in sorted(f for f in os.listdir(work_dir) if os.path.isfile(os.path.join(work_dir, f))): if out_file.startswith(os.path.basename(cigar_name)): out_nameids.append((os.path.basename(out_file), job.fileStore.writeGlobalFile(os.path.join(work_dir, out_file)))) diff --git a/src/cactus/progressive/cactus_progressive.py b/src/cactus/progressive/cactus_progressive.py index 61ff4d963..1baaa9607 100755 --- a/src/cactus/progressive/cactus_progressive.py +++ b/src/cactus/progressive/cactus_progressive.py @@ -484,7 +484,7 @@ def main(): if genome in event_set: if os.path.isdir(seq): tmpSeq = getTempFile() - catFiles([os.path.join(seq, subSeq) for subSeq in os.listdir(seq)], tmpSeq) + catFiles([os.path.join(seq, subSeq) for subSeq in sorted(os.listdir(seq))], tmpSeq) seq = tmpSeq seq = makeURL(seq) input_seq_id_map[genome] = toil.importFile(seq) diff --git a/src/cactus/refmap/cactus_graphmap.py b/src/cactus/refmap/cactus_graphmap.py index 710f90a0c..5d0bb4b9d 100644 --- a/src/cactus/refmap/cactus_graphmap.py +++ b/src/cactus/refmap/cactus_graphmap.py @@ -223,7 +223,7 @@ def graph_map(options): if genome != graph_event: if os.path.isdir(seq): tmpSeq = getTempFile() - catFiles([os.path.join(seq, subSeq) for subSeq in os.listdir(seq)], tmpSeq) + catFiles([os.path.join(seq, subSeq) for subSeq in sorted(os.listdir(seq))], tmpSeq) seq = tmpSeq seq = makeURL(seq) seq_id_map[genome] = toil.importFile(seq) diff --git a/src/cactus/refmap/cactus_graphmap_split.py b/src/cactus/refmap/cactus_graphmap_split.py index 6a4e8b30c..434359d1b 100644 --- a/src/cactus/refmap/cactus_graphmap_split.py +++ b/src/cactus/refmap/cactus_graphmap_split.py @@ -166,7 +166,7 @@ def cactus_graphmap_split(options): if genome in leaves: if os.path.isdir(seq): tmpSeq = getTempFile() - catFiles([os.path.join(seq, subSeq) for subSeq in os.listdir(seq)], tmpSeq) + catFiles([os.path.join(seq, subSeq) for subSeq in sorted(os.listdir(seq))], tmpSeq) seq = tmpSeq seq = makeURL(seq) input_seq_id_map[genome] = toil.importFile(seq) @@ -446,7 +446,7 @@ def split_gfa(job, config, gfa_id, paf_ids, ref_contigs, other_contig, reference cactus_call(parameters=cmd, work_dir=work_dir, job_memory=job.memory) output_id_map = {} - for out_name in os.listdir(work_dir): + for out_name in sorted(os.listdir(work_dir)): file_name, ext = os.path.splitext(out_name) if file_name.startswith(os.path.basename(out_prefix)) and ext in [".gfa", ".paf", ".fa_contigs"] and \ os.path.isfile(os.path.join(work_dir, file_name + ".fa_contigs")): diff --git a/src/cactus/refmap/cactus_minigraph.py b/src/cactus/refmap/cactus_minigraph.py index 037ef6f80..3492eecb4 100644 --- a/src/cactus/refmap/cactus_minigraph.py +++ b/src/cactus/refmap/cactus_minigraph.py @@ -184,7 +184,7 @@ def minigraph_construct_import_sequences(options, config_wrapper, input_seqfiles if genome != graph_event and genome in leaves: if os.path.isdir(seq): tmpSeq = getTempFile() - catFiles([os.path.join(seq, subSeq) for subSeq in os.listdir(seq)], tmpSeq) + catFiles([os.path.join(seq, subSeq) for subSeq in sorted(os.listdir(seq))], tmpSeq) seq = tmpSeq seq = makeURL(seq) input_seq_id_map[genome] = file_store.importFile(seq) diff --git a/src/cactus/refmap/cactus_pangenome.py b/src/cactus/refmap/cactus_pangenome.py index 92fb440ef..17aeaba5e 100644 --- a/src/cactus/refmap/cactus_pangenome.py +++ b/src/cactus/refmap/cactus_pangenome.py @@ -287,7 +287,7 @@ def main(): if genome != graph_event and genome in leaves: if os.path.isdir(seq): tmpSeq = getTempFile() - catFiles([os.path.join(seq, subSeq) for subSeq in os.listdir(seq)], tmpSeq) + catFiles([os.path.join(seq, subSeq) for subSeq in sorted(os.listdir(seq))], tmpSeq) seq = tmpSeq seq = makeURL(seq) input_seq_id_map[genome] = toil.importFile(seq) diff --git a/src/cactus/refmap/cactus_refmap.py b/src/cactus/refmap/cactus_refmap.py index 4c2f2c867..ed032daf0 100644 --- a/src/cactus/refmap/cactus_refmap.py +++ b/src/cactus/refmap/cactus_refmap.py @@ -320,7 +320,7 @@ def main(): if genome in event_set: if os.path.isdir(seq): tmpSeq = getTempFile() - catFiles([os.path.join(seq, subSeq) for subSeq in os.listdir(seq)], tmpSeq) + catFiles([os.path.join(seq, subSeq) for subSeq in sorted(os.listdir(seq))], tmpSeq) seq = tmpSeq seq = makeURL(seq) input_seq_id_map[genome] = toil.importFile(seq) diff --git a/src/cactus/setup/cactus_align.py b/src/cactus/setup/cactus_align.py index 68ab06976..d8cee3532 100644 --- a/src/cactus/setup/cactus_align.py +++ b/src/cactus/setup/cactus_align.py @@ -373,7 +373,7 @@ def make_align_job(options, toil, config_wrapper=None, chrom_name=None): if genome in event_set: if os.path.isdir(seq): tmpSeq = getTempFile() - catFiles([os.path.join(seq, subSeq) for subSeq in os.listdir(seq)], tmpSeq) + catFiles([os.path.join(seq, subSeq) for subSeq in sorted(os.listdir(seq))], tmpSeq) seq = tmpSeq seq = makeURL(seq) input_seq_id_map[genome] = toil.importFile(seq)