Make cactus reproducible

api/impl/cactusDisk.c

@@ -253,7 +253,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
@@ -269,6 +297,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;
 }

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.
  */

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; j<stList_length(flowers); j++) {
         Flower *flower = stList_get(flowers, j);
 
+        cactusDisk_pushNameInterval(cactusDisk, nameBase + nameOffsets[j], nameOffsets[j+1] - nameOffsets[j]);
+        st_randomSeed(flower_getName(flower)); // Any random choices below are the flower's own, not the thread's
+
         // These are all variables used by the filter fns
         FilterArgs *fa = st_calloc(1, sizeof(FilterArgs));
         fa->minimumIngroupDegree = 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

bar/impl/poaBarAligner.c

@@ -1090,7 +1090,12 @@ int caps_comp_by_adjacency_length(const void *a, const void *b) {
     int64_t 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));
 }
 
 /*

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}

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);
 }
 

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);

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) {

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<stList_length(flowers); i++) {
+    for(int64_t i=0; i<flowerNumber; i++) {
         Flower *flower = stList_get(flowers, i);
         st_logDebug("Processing flower %" PRIi64 "\n", flower_getName(flower));
+        cactusDisk_pushNameInterval(cactusDisk, nameBase + nameOffsets[i], nameOffsets[i+1] - nameOffsets[i]);
+        st_randomSeed(flower_getName(flower)); // The reference ordering is randomised: make it the flower's own
+        // sequence of random choices, so it does not depend on the thread it runs on
         buildReferenceTopDown(flower, referenceEventString, permutations, matchingAlgorithm, temperatureFn, theta,
                               phi, maxWalkForCalculatingZ, ignoreUnalignedGaps, wiggle, numberOfNsForScaffoldGap,
                               minNumberOfSequencesToSupportAdjacency, makeScaffolds,
                               minimumNestedBasesToBreakAdjacency, maximumChainBasesToBreakAdjacency);
+        cactusDisk_popNameInterval(cactusDisk);
     }
+    free(nameOffsets);
 }

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)

src/cactus/paf/local_alignment.py

@@ -408,9 +408,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
+        # <event>_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)