From 44c96a75e7efb7b88c75389e1201f98ba73690e5 Mon Sep 17 00:00:00 2001 From: Codex Date: Wed, 20 May 2026 15:20:14 +0000 Subject: [PATCH 1/3] Add genomeGenerate benchmark harness Add a local benchmark script for genomeGenerate that records wall time, CPU and I/O samples, selected index-build settings, and per-stage timings from Log.out. This keeps performance validation reproducible while leaving STAR runtime behavior unchanged. --- .../tests/scripts/benchmarkGenomeGenerate.sh | 436 ++++++++++++++++++ 1 file changed, 436 insertions(+) create mode 100755 extras/tests/scripts/benchmarkGenomeGenerate.sh diff --git a/extras/tests/scripts/benchmarkGenomeGenerate.sh b/extras/tests/scripts/benchmarkGenomeGenerate.sh new file mode 100755 index 00000000..ba9885ec --- /dev/null +++ b/extras/tests/scripts/benchmarkGenomeGenerate.sh @@ -0,0 +1,436 @@ +#!/usr/bin/env bash +set -euo pipefail + +# Genome-generation benchmark and idempotence harness. +# +# Typical use after a genomeGenerate optimization: +# make -C source STAR +# BASELINE_STAR_BIN=/path/to/baseline/STAR THREADS="1 4 8 16" \ +# extras/tests/scripts/benchmarkGenomeGenerate.sh +# +# To compare against both a recent optimized baseline and a raw/base STAR build: +# BASELINE_STAR_BIN=/path/to/optimized/STAR BASE_STAR_BIN=/path/to/base/STAR \ +# THREADS="64 128" extras/tests/scripts/benchmarkGenomeGenerate.sh +# +# For performance measurements, prefer a realistic FASTA: +# BENCH_FASTA=/path/to/genome.fa LIMIT_GENOME_GENERATE_RAM=31000000000 \ +# THREADS="1 4 8 16 32" extras/tests/scripts/benchmarkGenomeGenerate.sh +# +# To benchmark a splice-junction-augmented genome index: +# BENCH_FASTA=/path/to/genome.fa BENCH_GTF=/path/to/genes.gtf SJDB_OVERHANG=93 \ +# THREADS="32 64" extras/tests/scripts/benchmarkGenomeGenerate.sh +# +# To pass multiple FASTA files, use a space-separated list: +# BENCH_FASTA_FILES="/path/to/genome.fa /path/to/spikeins.fa" ... +# +# The built-in FASTA is intentionally tiny. It is for fast byte-identity checks, +# not for meaningful wall-clock speedup measurements. +# Use SYNTHETIC_REPEAT=N to scale it up deterministically when BENCH_FASTA is +# not set. + +script_dir="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +repo_root="$(cd "${script_dir}/../../.." && pwd)" + +star_bin="${STAR_BIN:-${repo_root}/source/STAR}" +baseline_star_bin="${BASELINE_STAR_BIN:-}" +base_star_bin="${BASE_STAR_BIN:-}" +out_root="${OUT_DIR:-/tmp/star-genome-generate-bench.$(date +%Y%m%d_%H%M%S).$$}" +threads_list="${THREADS:-1 4}" +runs_per_thread="${RUNS_PER_THREAD:-2}" +baseline_runs_per_thread="${BASELINE_RUNS_PER_THREAD:-${runs_per_thread}}" +base_runs_per_thread="${BASE_RUNS_PER_THREAD:-${baseline_runs_per_thread}}" +genome_saindex_nbases="${GENOME_SAINDEX_NBASES:-2}" +genome_chrbin_nbits="${GENOME_CHRBIN_NBITS:-4}" +limit_genome_ram="${LIMIT_GENOME_GENERATE_RAM:-20000}" +bench_gtf="${BENCH_GTF:-}" +sjdb_overhang="${SJDB_OVERHANG:-}" +time_bin="${TIME_BIN:-/usr/bin/time}" +synthetic_repeat="${SYNTHETIC_REPEAT:-1}" +iostat_bin="${IOSTAT_BIN:-$(command -v iostat || true)}" +iostat_interval="${IOSTAT_INTERVAL:-5}" +cpu_sample_interval="${CPU_SAMPLE_INTERVAL:-1}" + +if [[ ! -x "${star_bin}" ]]; then + echo "ERROR: STAR_BIN is not executable: ${star_bin}" >&2 + exit 1 +fi + +if [[ -n "${baseline_star_bin}" && ! -x "${baseline_star_bin}" ]]; then + echo "ERROR: BASELINE_STAR_BIN is not executable: ${baseline_star_bin}" >&2 + exit 1 +fi + +if [[ -n "${base_star_bin}" && ! -x "${base_star_bin}" ]]; then + echo "ERROR: BASE_STAR_BIN is not executable: ${base_star_bin}" >&2 + exit 1 +fi + +if [[ ! -x "${time_bin}" ]]; then + echo "ERROR: TIME_BIN is not executable: ${time_bin}" >&2 + exit 1 +fi + +if [[ -n "${bench_gtf}" && ! -r "${bench_gtf}" ]]; then + echo "ERROR: BENCH_GTF is not readable: ${bench_gtf}" >&2 + exit 1 +fi + +if [[ -n "${sjdb_overhang}" && -z "${bench_gtf}" ]]; then + echo "ERROR: SJDB_OVERHANG requires BENCH_GTF" >&2 + exit 1 +fi + +mkdir -p "${out_root}" +out_root="$(cd "${out_root}" && pwd)" + +fasta="${BENCH_FASTA:-${out_root}/synthetic.fa}" +bench_fasta_files="${BENCH_FASTA_FILES:-${BENCH_FASTA:-}}" +fasta_files=() + +make_synthetic_fasta() { + local out_fasta="$1" + local repeat_i + { + echo ">chrA" + for ((repeat_i=0; repeat_ichrB" + for ((repeat_i=0; repeat_i "${out_fasta}" +} + +if [[ -z "${bench_fasta_files}" ]]; then + make_synthetic_fasta "${fasta}" + fasta_files=("${fasta}") +else + read -r -a fasta_files <<< "${bench_fasta_files}" +fi + +for fasta_file in "${fasta_files[@]}"; do + if [[ ! -r "${fasta_file}" ]]; then + echo "ERROR: FASTA is not readable: ${fasta_file}" >&2 + exit 1 + fi +done + +if [[ "${#fasta_files[@]}" -eq 0 ]]; then + echo "ERROR: no FASTA files were provided" >&2 + exit 1 +fi + +extract_stage_times() { + local log_path="$1" + local stage_path="$2" + local log_year + log_year="$(date +%Y)" + + awk -v year="${log_year}" ' + BEGIN { + month["Jan"]=1; month["Feb"]=2; month["Mar"]=3; month["Apr"]=4; + month["May"]=5; month["Jun"]=6; month["Jul"]=7; month["Aug"]=8; + month["Sep"]=9; month["Oct"]=10; month["Nov"]=11; month["Dec"]=12; + print "stage\tstart_time\tend_time\tseconds"; + } + + function timestamp(line, mon, day, hour, minute, second) { + mon=substr(line,1,3); + if (!(mon in month)) { + return ""; + } + day=substr(line,5,2)+0; + hour=substr(line,8,2)+0; + minute=substr(line,11,2)+0; + second=substr(line,14,2)+0; + return mktime(year " " month[mon] " " day " " hour " " minute " " second); + } + + function mark(key,line, parsed_time) { + if (!(key in times)) { + parsed_time=timestamp(line); + if (parsed_time!="") { + times[key]=parsed_time; + labels[key]=substr(line,1,15); + } + } + } + + function emit(stage,start_key,end_key) { + if ((start_key in times) && (end_key in times)) { + printf "%s\t%s\t%s\t%d\n", stage, labels[start_key], labels[end_key], times[end_key]-times[start_key]; + } else { + printf "%s\tNA\tNA\tNA\n", stage; + } + } + + /starting to sort Suffix Array/ { mark("sa_sort_start",$0); next } + /sorting Suffix Array chunks/ { mark("sa_chunk_sort",$0); next } + /packing SA from RAM chunks|loading chunks from disk, packing SA/ { mark("sa_pack",$0); next } + /finished generating suffix array/ { mark("sa_done",$0); next } + /generating Suffix Array index/ { mark("saindex_start",$0); next } + /completed Suffix Array index/ { mark("saindex_done",$0); next } + /Finished preparing junctions/ { mark("sj_prepare_done",$0); next } + /inserting junctions into the genome indices/ { mark("sj_insert_start",$0); next } + /Finished SA search:/ { mark("sj_sa_search_done",$0); next } + /Finished sorting SA indicesL/ { mark("sj_sort_done",$0); next } + /Finished inserting junction indices/ { mark("sj_insert_indices_done",$0); next } + /Finished SAi$/ { mark("sj_sai_done",$0); next } + /writing Genome to disk/ { mark("write_genome_start",$0); next } + /writing Suffix Array to disk/ { mark("write_sa_start",$0); next } + /writing SAindex to disk/ { mark("write_saindex_start",$0); next } + /finished successfully/ { mark("finish_success",$0); next } + + END { + emit("sa_prefix_plan_seconds","sa_sort_start","sa_chunk_sort"); + emit("sa_chunk_sort_seconds","sa_chunk_sort","sa_pack"); + emit("sa_pack_seconds","sa_pack","sa_done"); + emit("saindex_seconds","saindex_start","saindex_done"); + emit("sj_prepare_seconds","saindex_done","sj_prepare_done"); + emit("sj_sa_search_seconds","sj_insert_start","sj_sa_search_done"); + emit("sj_sort_seconds","sj_sa_search_done","sj_sort_done"); + emit("sj_insert_indices_seconds","sj_sort_done","sj_insert_indices_done"); + emit("sj_sai_seconds","sj_insert_indices_done","sj_sai_done"); + emit("write_genome_seconds","write_genome_start","write_sa_start"); + emit("write_sa_seconds","write_sa_start","write_saindex_start"); + emit("write_saindex_seconds","write_saindex_start","finish_success"); + } + ' "${log_path}" > "${stage_path}" +} + +stage_metric() { + local stage_path="$1" + local stage_name="$2" + awk -F '\t' -v stage="${stage_name}" '$1==stage {print $4; found=1; exit} END {if (!found) print "NA"}' "${stage_path}" +} + +summary="${out_root}/summary.tsv" +printf "label\tthreads\trun\treal_seconds\tuser_seconds\tsys_seconds\tmax_rss_kb\tsa_chunk_fill_strategy\tsa_chunk_sort_prefix_length\tsa_chunk_sort_granularity\tsa_chunk_storage_strategy\tsa_chunk_system_available_bytes\tsa_chunk_ram_peak_bytes\tsa_all_scatter_available_bytes\tsa_all_scatter_required_bytes\tsa_batched_fill_batches\tsaindex_traversal_strategy\tsaindex_event_count\tsa_prefix_plan_seconds\tsa_chunk_sort_seconds\tsa_pack_seconds\tsaindex_seconds\tsj_prepare_seconds\tsj_sa_search_seconds\tsj_sort_seconds\tsj_insert_indices_seconds\tsj_sai_seconds\twrite_genome_seconds\twrite_sa_seconds\twrite_saindex_seconds\tstage_times_log\tgenome_dir\tio_log\tcpu_log\n" > "${summary}" + +run_dir="" + +run_generate() { + local label="$1" + local bin="$2" + local threads="$3" + local run_id="$4" + run_dir="${out_root}/${label}_t${threads}_r${run_id}" + local log_file="${out_root}/${label}_t${threads}_r${run_id}.log" + local time_file="${out_root}/${label}_t${threads}_r${run_id}.time" + local io_file="${out_root}/${label}_t${threads}_r${run_id}.iostat" + local cpu_file="${out_root}/${label}_t${threads}_r${run_id}.threads" + local stage_file="${out_root}/${label}_t${threads}_r${run_id}.stage_times.tsv" + local iostat_pid="" + local cpu_sampler_pid="" + local star_wrapper_pid="" + + if [[ -e "${run_dir}" ]]; then + echo "ERROR: output directory already exists: ${run_dir}" >&2 + exit 1 + fi + mkdir -p "${run_dir}" + + if [[ -n "${iostat_bin}" && -x "${iostat_bin}" ]]; then + "${iostat_bin}" -xz "${iostat_interval}" > "${io_file}" 2>&1 & + iostat_pid="$!" + else + printf "iostat unavailable\n" > "${io_file}" + fi + + printf "timestamp\tpid\ttid\tpsr\tpcpu\tstat\tcomm\targs\n" > "${cpu_file}" + + local star_args=( + --runMode genomeGenerate + --runThreadN "${threads}" + --genomeDir "${run_dir}" + --genomeFastaFiles "${fasta_files[@]}" + --genomeSAindexNbases "${genome_saindex_nbases}" + --genomeChrBinNbits "${genome_chrbin_nbits}" + --limitGenomeGenerateRAM "${limit_genome_ram}" + ) + + if [[ -n "${bench_gtf}" ]]; then + star_args+=(--sjdbGTFfile "${bench_gtf}") + fi + + if [[ -n "${sjdb_overhang}" ]]; then + star_args+=(--sjdbOverhang "${sjdb_overhang}") + fi + + set +e + ( + cd "${out_root}" + "${time_bin}" -f "%e\t%U\t%S\t%M" -o "${time_file}" \ + "${bin}" "${star_args[@]}" \ + > "${log_file}" 2>&1 + ) & + star_wrapper_pid="$!" + + ( + while kill -0 "${star_wrapper_pid}" 2>/dev/null; do + local_now="$(date -Ins)" + ps -eLo pid,tid,psr,pcpu,stat,comm,args \ + | awk -v now="${local_now}" -v dir="${run_dir}" \ + 'NR>1 && $6=="STAR" && index($0,dir)>0 {printf "%s\t%s\t%s\t%s\t%s\t%s\t%s\t", now, $1, $2, $3, $4, $5, $6; for (i=7; i<=NF; i++) printf "%s%s", $i, (i> "${cpu_file}" + sleep "${cpu_sample_interval}" + done + ) & + cpu_sampler_pid="$!" + + wait "${star_wrapper_pid}" + local star_status="$?" + set -e + + if [[ -n "${cpu_sampler_pid}" ]]; then + kill "${cpu_sampler_pid}" 2>/dev/null || true + wait "${cpu_sampler_pid}" 2>/dev/null || true + fi + + if [[ -n "${iostat_pid}" ]]; then + kill "${iostat_pid}" 2>/dev/null || true + wait "${iostat_pid}" 2>/dev/null || true + fi + + if [[ "${star_status}" -ne 0 ]]; then + echo "ERROR: STAR genomeGenerate failed for ${label} t${threads} run ${run_id}; see ${log_file}" >&2 + exit "${star_status}" + fi + + extract_stage_times "${run_dir}/Log.out" "${stage_file}" + + local real_s user_s sys_s max_rss sa_chunk_fill_strategy sa_chunk_sort_prefix_length sa_chunk_sort_granularity sa_chunk_storage_strategy sa_chunk_system_available_bytes sa_chunk_ram_peak_bytes sa_all_scatter_available_bytes sa_all_scatter_required_bytes sa_batched_fill_batches saindex_traversal_strategy saindex_event_count + local sa_prefix_plan_seconds sa_chunk_sort_seconds sa_pack_seconds saindex_seconds sj_prepare_seconds sj_sa_search_seconds sj_sort_seconds sj_insert_indices_seconds sj_sai_seconds write_genome_seconds write_sa_seconds write_saindex_seconds + read -r real_s user_s sys_s max_rss < "${time_file}" + sa_chunk_fill_strategy="$(awk -F': ' '/SA chunk fill strategy:/ {print $2; found=1; exit} END {if (!found) print "NA"}' "${run_dir}/Log.out")" + sa_chunk_sort_prefix_length="$(awk -F': ' '/SA chunk sort prefix length:/ {print $2; found=1; exit} END {if (!found) print "NA"}' "${run_dir}/Log.out")" + sa_chunk_sort_granularity="$(awk -F': ' '/SA chunk sort granularity:/ {print $2; found=1; exit} END {if (!found) print "NA"}' "${run_dir}/Log.out")" + sa_chunk_storage_strategy="$(awk -F': ' '/SA chunk storage strategy:/ {print $2; found=1; exit} END {if (!found) print "NA"}' "${run_dir}/Log.out")" + sa_chunk_system_available_bytes="$(awk '/SA chunk system available bytes:/ {gsub(/;/,""); print $6; found=1; exit} END {if (!found) print "NA"}' "${run_dir}/Log.out")" + sa_chunk_ram_peak_bytes="$(awk '/SA chunk retained bytes:/ {gsub(/;/,""); print $9; found=1; exit} END {if (!found) print "NA"}' "${run_dir}/Log.out")" + read -r sa_all_scatter_available_bytes sa_all_scatter_required_bytes < <(awk '/SA chunk all-scatter available bytes:/ {gsub(/;/,""); print $6, $10; found=1; exit} END {if (!found) print "NA NA"}' "${run_dir}/Log.out") + sa_batched_fill_batches="$(awk -F': ' '/SA chunk estimated batched-fill batches:/ {print $2; found=1; exit} END {if (!found) print "NA"}' "${run_dir}/Log.out")" + saindex_traversal_strategy="$(awk -F': ' '/SAindex traversal strategy:/ {print $2; found=1; exit} END {if (!found) print "NA"}' "${run_dir}/Log.out")" + saindex_event_count="$(awk -F': ' '/SAindex event count:/ {print $2; found=1; exit} END {if (!found) print "NA"}' "${run_dir}/Log.out")" + sa_prefix_plan_seconds="$(stage_metric "${stage_file}" "sa_prefix_plan_seconds")" + sa_chunk_sort_seconds="$(stage_metric "${stage_file}" "sa_chunk_sort_seconds")" + sa_pack_seconds="$(stage_metric "${stage_file}" "sa_pack_seconds")" + saindex_seconds="$(stage_metric "${stage_file}" "saindex_seconds")" + sj_prepare_seconds="$(stage_metric "${stage_file}" "sj_prepare_seconds")" + sj_sa_search_seconds="$(stage_metric "${stage_file}" "sj_sa_search_seconds")" + sj_sort_seconds="$(stage_metric "${stage_file}" "sj_sort_seconds")" + sj_insert_indices_seconds="$(stage_metric "${stage_file}" "sj_insert_indices_seconds")" + sj_sai_seconds="$(stage_metric "${stage_file}" "sj_sai_seconds")" + write_genome_seconds="$(stage_metric "${stage_file}" "write_genome_seconds")" + write_sa_seconds="$(stage_metric "${stage_file}" "write_sa_seconds")" + write_saindex_seconds="$(stage_metric "${stage_file}" "write_saindex_seconds")" + printf "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n" \ + "${label}" "${threads}" "${run_id}" "${real_s}" "${user_s}" "${sys_s}" "${max_rss}" \ + "${sa_chunk_fill_strategy}" "${sa_chunk_sort_prefix_length}" "${sa_chunk_sort_granularity}" \ + "${sa_chunk_storage_strategy}" "${sa_chunk_system_available_bytes}" "${sa_chunk_ram_peak_bytes}" \ + "${sa_all_scatter_available_bytes}" "${sa_all_scatter_required_bytes}" "${sa_batched_fill_batches}" "${saindex_traversal_strategy}" "${saindex_event_count}" \ + "${sa_prefix_plan_seconds}" "${sa_chunk_sort_seconds}" "${sa_pack_seconds}" "${saindex_seconds}" \ + "${sj_prepare_seconds}" "${sj_sa_search_seconds}" "${sj_sort_seconds}" "${sj_insert_indices_seconds}" "${sj_sai_seconds}" \ + "${write_genome_seconds}" "${write_sa_seconds}" "${write_saindex_seconds}" \ + "${stage_file}" "${run_dir}" "${io_file}" "${cpu_file}" >> "${summary}" +} + +compare_core_outputs() { + local dir_a="$1" + local dir_b="$2" + local label="$3" + local file + + for file in Genome SA SAindex chrName.txt chrStart.txt chrLength.txt chrNameLength.txt sjdbList.out.tab; do + if [[ -e "${dir_a}/${file}" || -e "${dir_b}/${file}" ]]; then + if ! cmp -s "${dir_a}/${file}" "${dir_b}/${file}"; then + echo "ERROR: ${label}: ${file} differs between ${dir_a} and ${dir_b}" >&2 + exit 1 + fi + fi + done +} + +first_candidate_dir="" + +for threads in ${threads_list}; do + candidate_reference_dir="" + + run_id=1 + while [[ "${run_id}" -le "${runs_per_thread}" ]]; do + run_generate candidate "${star_bin}" "${threads}" "${run_id}" + candidate_dir="${run_dir}" + + if [[ -z "${candidate_reference_dir}" ]]; then + candidate_reference_dir="${candidate_dir}" + else + compare_core_outputs "${candidate_reference_dir}" "${candidate_dir}" "candidate repeat t${threads}" + fi + + run_id=$((run_id+1)) + done + + if [[ -z "${first_candidate_dir}" ]]; then + first_candidate_dir="${candidate_reference_dir}" + else + compare_core_outputs "${first_candidate_dir}" "${candidate_reference_dir}" "candidate thread-count t${threads}" + fi + + if [[ -n "${baseline_star_bin}" ]]; then + baseline_reference_dir="" + run_id=1 + while [[ "${run_id}" -le "${baseline_runs_per_thread}" ]]; do + run_generate baseline "${baseline_star_bin}" "${threads}" "${run_id}" + baseline_dir="${run_dir}" + + if [[ -z "${baseline_reference_dir}" ]]; then + baseline_reference_dir="${baseline_dir}" + else + compare_core_outputs "${baseline_reference_dir}" "${baseline_dir}" "baseline repeat t${threads}" + fi + + run_id=$((run_id+1)) + done + + compare_core_outputs "${baseline_reference_dir}" "${candidate_reference_dir}" "baseline vs candidate t${threads}" + fi + + if [[ -n "${base_star_bin}" ]]; then + base_reference_dir="" + run_id=1 + while [[ "${run_id}" -le "${base_runs_per_thread}" ]]; do + run_generate base "${base_star_bin}" "${threads}" "${run_id}" + base_dir="${run_dir}" + + if [[ -z "${base_reference_dir}" ]]; then + base_reference_dir="${base_dir}" + else + compare_core_outputs "${base_reference_dir}" "${base_dir}" "base repeat t${threads}" + fi + + run_id=$((run_id+1)) + done + + compare_core_outputs "${base_reference_dir}" "${candidate_reference_dir}" "base vs candidate t${threads}" + fi +done + +echo "Benchmark complete: ${out_root}" +echo "Summary: ${summary}" From 2d12a0c62a8dca4e3fdefbae3369447202b94403 Mon Sep 17 00:00:00 2001 From: Codex Date: Wed, 20 May 2026 15:20:19 +0000 Subject: [PATCH 2/3] Parallelize SAindex and junction indexing Use parallel traversal for SAindex generation and parallel bucket sorting for junction insertion indices. Keep deterministic output ordering while reducing the remaining annotation-heavy genomeGenerate stages. --- source/genomeSAindex.cpp | 146 +++++++++++++++++++++++++++---------- source/sjdbBuildIndex.cpp | 149 +++++++++++++++++++++++++++++++++----- 2 files changed, 239 insertions(+), 56 deletions(-) diff --git a/source/genomeSAindex.cpp b/source/genomeSAindex.cpp index 140b171e..7e180cb3 100644 --- a/source/genomeSAindex.cpp +++ b/source/genomeSAindex.cpp @@ -3,6 +3,93 @@ #include "SuffixArrayFuns.h" #include "ErrorWarning.h" +struct SAindexEvent { + uint isa; + uint indFull; + int iL4; +}; + +static bool SAindexEventEqual(uint indFull1, int iL41, uint indFull2, int iL42) +{ + return indFull1==indFull2 && iL41==iL42; +}; + +static void SAindexProcessEvent(PackedArray &SAi, Genome &mapGen, uint isa, uint indFull, int iL4, uint* ind0) +{ + for (uint iL=0; iL < mapGen.pGe.gSAindexNbases; iL++) {//calculate index + + uint indPref = indFull >> (2*(mapGen.pGe.gSAindexNbases-1-iL)); + + if ( (int)iL==iL4 ) {//this suffix contains N and does not belong in SAi + for (uint iL1=iL; iL1 < mapGen.pGe.gSAindexNbases; iL1++) { + SAi.writePacked(mapGen.genomeSAindexStart[iL1]+ind0[iL1],SAi[mapGen.genomeSAindexStart[iL1]+ind0[iL1]] | mapGen.SAiMarkNmaskC); + }; + break;//break the iL cycle + }; + + if ( indPref > ind0[iL] || isa==0 ) {//new && good index, record it + + SAi.writePacked(mapGen.genomeSAindexStart[iL]+indPref, isa); + + for (uint ii=ind0[iL]+1; iilogMain, EXIT_CODE_INPUT_FILES, mapGen.P); + }; + }; +}; + +static void SAindexScanEventsParallel(char *G, PackedArray &SA, Parameters &P, Genome &mapGen, uint iSA1, uint iSA2, vector &events) +{ + const uint nSA=iSA2-iSA1+1; + const uint chunkN=min(nSA, max(1, P.runThreadN*4)); + const int threadN=(int) min(P.runThreadN, chunkN); + vector > eventsByChunk(chunkN); + + #pragma omp parallel for num_threads(threadN) schedule(static) + for (uint iChunk=0; iChunk &chunkEvents=eventsByChunk[iChunk]; + chunkEvents.reserve((iEnd-iStart)/32+1); + + int iL4prev=-1; + uint indFullPrev=0; + if (iStart>iSA1) { + indFullPrev=funCalcSAiFromSA(G,SA,mapGen,iStart-1,mapGen.pGe.gSAindexNbases,iL4prev); + }; + + for (uint isa=iStart; isa=16 && iSA1==0 && iSA2+1==mapGen.nSA; + + P.inOut->logMain << "SAindex traversal strategy: " << (parallelEvents ? "parallel-events" : "skip-search") << "\n" << flush; + uint* ind0=new uint[mapGen.pGe.gSAindexNbases]; for (uint ii=0; ii events; + SAindexScanEventsParallel(G, SA, P, mapGen, iSA1, iSA2, events); - uint isaStep=mapGen.nSA/(1llu<<(2*mapGen.pGe.gSAindexNbases))+1; + P.inOut->logMain << "SAindex event count: " << events.size() << "\n" << flush; + for (uint iEvent=0; iEvent> (2*(mapGen.pGe.gSAindexNbases-1-iL)); - - if ( (int)iL==iL4 ) {//this suffix contains N and does not belong in SAi - for (uint iL1=iL; iL1 < mapGen.pGe.gSAindexNbases; iL1++) { - SAi.writePacked(mapGen.genomeSAindexStart[iL1]+ind0[iL1],SAi[mapGen.genomeSAindexStart[iL1]+ind0[iL1]] | mapGen.SAiMarkNmaskC); - }; - break;//break the iL cycle - }; - - if ( indPref > ind0[iL] || isa==0 ) {//new && good index, record it - - SAi.writePacked(mapGen.genomeSAindexStart[iL]+indPref, isa); - - for (uint ii=ind0[iL]+1; iilogMain, EXIT_CODE_INPUT_FILES, P); - }; - }; + uint isa=iSA1; + int iL4; + uint indFull=funCalcSAiFromSA(G,SA,mapGen,isa,mapGen.pGe.gSAindexNbases,iL4); + while (isa<=iSA2) {//for all suffixes + SAindexProcessEvent(SAi, mapGen, isa, indFull, iL4, ind0); - //find next index not equal to the current one - funSAiFindNextIndex(G, SA, isaStep, isa, indFull, iL4, mapGen);//indFull and iL4 have been already defined at the previous step + //find next index not equal to the current one + funSAiFindNextIndex(G, SA, isaStep, isa, indFull, iL4, mapGen);//indFull and iL4 have been already defined at the previous step - };//isa cycle + };//isa cycle + }; for (uint iL=0; iL < mapGen.pGe.gSAindexNbases; iL++) {//fill up unfilled indexes for (uint ii=mapGen.genomeSAindexStart[iL]+ind0[iL]+1; iinInd) { + bucketN=1; + while ((bucketN<<1)<=nInd) bucketN <<= 1; + }; + if (bucketN<2) { + qsort((void*) indArray, nInd, 2*sizeof(uint64), funCompareUintAndSuffixes); + return false; + }; + + const uint64 indBytes=2*nInd*sizeof(uint64); + const uint64 auxBytes=indBytes + (uint64) P.runThreadN*bucketN*2*sizeof(uint64) + (bucketN+1)*sizeof(uint64); + if (P.limitGenomeGenerateRAM>0 && auxBytes>P.limitGenomeGenerateRAM/4) { + qsort((void*) indArray, nInd, 2*sizeof(uint64), funCompareUintAndSuffixes); + return false; + }; + + vector bucketCount((uint) P.runThreadN*bucketN,0); + #pragma omp parallel num_threads(P.runThreadN) + { + uint tid=(uint) omp_get_thread_num(); + uint64* threadBucketCount=bucketCount.data()+tid*bucketN; + #pragma omp for schedule(static) + for (uint ii=0; ii bucketStart(bucketN+1,0); + for (uint iBucket=0; iBucket bucketThreadStart((uint) P.runThreadN*bucketN,0); + for (uint iBucket=0; iBucketlogMain << timeMonthDayTime(rawtime) << " Finished sorting SA indicesL nInd="<logMain << timeMonthDayTime(rawtime) << " Finished sorting SA indicesL nInd="< sjdbFirstN(2*nGsj+1,mapGen.pGe.gSAindexNbases); + uint nextN=2*nGsj; + for (int64 ii=2*(int64)nGsj; ii>=0; ii--) { + if (Gsj[ii]>3) { + nextN=(uint) ii; + }; + uint nDist=nextN-(uint) ii; + if (nDist=mapGen.pGe.gSAindexNbases) { + continue; + }; + sjdbSAiNmarkN++; + int64 ind1=0; - for (uint iL=0; iL < mapGen.pGe.gSAindexNbases; iL++) { - uint g=(uint) Gsj[indArray[2*isj+1]+iL]; + for (uint iL0=0; iL03) {//this iSA contains N, need to mark the previous - for (uint iL1=iL; iL1 < mapGen.pGe.gSAindexNbases; iL1++) { - ind1+=3; - int64 ind2=mapGen.genomeSAindexStart[iL1]+ind1; - for (; ind2>=0; ind2--) {//find previous index that is not absent - if ( (SAi[ind2] & mapGen.SAiMarkAbsentMaskC)==0 ) { - break; - }; - }; - SAi.writePacked(ind2,SAi[ind2] | mapGen.SAiMarkNmaskC); - ind1 <<= 2; + ind1 += (uint) Gsj[isjG+iL0]; + }; + + ind1 <<= 2; + for (uint iL1=iL; iL1 < mapGen.pGe.gSAindexNbases; iL1++) { + ind1+=3; + int64 ind2=mapGen.genomeSAindexStart[iL1]+ind1; + for (; ind2>=0; ind2--) {//find previous index that is not absent + if ( (SAi[ind2] & mapGen.SAiMarkAbsentMaskC)==0 ) { + break; }; - break; - } else { - ind1 += g; }; + SAi.writePacked(ind2,SAi[ind2] | mapGen.SAiMarkNmaskC); + ind1 <<= 2; }; }; time ( &rawtime ); + P.inOut->logMain << timeMonthDayTime(rawtime) << " Finished SAi N marking: candidates=" << sjdbSAiNmarkN + << " strategy=first-N-cache" <logMain << timeMonthDayTime(rawtime) << " Finished SAi" < Date: Wed, 20 May 2026 15:20:25 +0000 Subject: [PATCH 3/3] Improve genomeGenerate suffix-array construction Reduce genomeGenerate suffix-array build time by batching prefix-bin fills, retaining sorted chunks in RAM when memory allows, splitting very large prefix bins into ordered sub-bins, and using a comparator that can skip already-known prefix words. Optional SA sort profiling remains gated behind STAR_PROFILE_SA_SORT=1. --- source/Genome_genomeGenerate.cpp | 569 ++++++++++++++++++++++++++++--- 1 file changed, 516 insertions(+), 53 deletions(-) diff --git a/source/Genome_genomeGenerate.cpp b/source/Genome_genomeGenerate.cpp index 213edf7f..49051f45 100755 --- a/source/Genome_genomeGenerate.cpp +++ b/source/Genome_genomeGenerate.cpp @@ -1,4 +1,6 @@ #include +#include +#include #include "Genome.h" @@ -25,18 +27,34 @@ char* globalG; uint globalL; +static uint systemAvailableMemoryBytes() +{ +#ifdef __linux__ + ifstream memInfo("/proc/meminfo"); + string field, unit; + uint value=0; -inline int funCompareSuffixes ( const void *a, const void *b){ + while (memInfo >> field >> value >> unit) { + if (field=="MemAvailable:") { + return value*1024LLU; + }; + }; +#endif + return 0; +}; + + +inline int funCompareSuffixesFromWord ( const void *a, const void *b, uint wordStart){ - uint *ga=(uint*)((globalG-7LLU)+(*((uint*)a))); - uint *gb=(uint*)((globalG-7LLU)+(*((uint*)b))); + uint *ga=(uint*)((globalG-7LLU)+(*((uint*)a)))-wordStart; + uint *gb=(uint*)((globalG-7LLU)+(*((uint*)b)))-wordStart; uint jj=0; int ii=0; uint va=0,vb=0; uint8 *va1, *vb1; - while (jj < globalL) { + while (jj+wordStart < globalL) { va=*(ga-jj); vb=*(gb-jj); @@ -88,6 +106,14 @@ inline int funCompareSuffixes ( const void *a, const void *b){ }; }; +inline int funCompareSuffixes ( const void *a, const void *b){ + return funCompareSuffixesFromWord(a,b,0); +}; + +inline int funCompareSuffixesSkipFirstWord ( const void *a, const void *b){ + return funCompareSuffixesFromWord(a,b,1); +}; + inline uint funG2strLocus (uint SAstr, uint const N, char const GstrandBit, uint const GstrandMask) { bool strandG = (SAstr>>GstrandBit) == 0; SAstr &= GstrandMask; @@ -95,6 +121,43 @@ inline uint funG2strLocus (uint SAstr, uint const N, char const GstrandBit, uint return SAstr; }; +inline uint funSAsortPrefix(char* G, uint ii, uint prefixLength) { + uint prefix=0; + for (uint jj=0; jj5) g=5; + prefix=prefix*6+g; + }; + return prefix; +}; + +inline uint funSAsortPrefixAtOffset(char* G, uint ii, uint offset, uint prefixLength) { + uint prefix=0; + for (uint jj=0; jj5) g=5; + prefix=prefix*6+g; + }; + return prefix; +}; + +inline bool funSAsortPrefixFirstWordHas5(uint prefix, uint prefixLength) { + uint wordBases=(uint) sizeof(uint); + if (prefixLength>wordBases) { + for (uint jj=0; jjnG1alloc ? P.limitGenomeGenerateRAM-nG1alloc : 0; + if (P.runThreadN>=16 && nGenome>=50000000) { + uint indPrefN8=1; + for (uint ii=0;ii<8;ii++) indPrefN8 *= 6; + const uint64 indPrefTempBytes8=(uint64) P.runThreadN*indPrefN8*2*sizeof(uint) + (indPrefN8+1)*sizeof(uint); + if (indPrefTempBytes8 < saAvailableBytesForPrefix/4) { + indPrefLen=8; + }; + }; + uint indPrefN=1; + for (uint ii=0;iiindPrefMaxCount) { + indPrefMaxCount=indPrefCount[ii]; + indPrefMaxIndex=ii; + }; + }; + }; + }; + uint saChunkSize=(P.limitGenomeGenerateRAM-nG1alloc)/8/P.runThreadN; //number of SA indexes per chunk saChunkSize=saChunkSize*6/10; //allow extra space for qsort //uint saChunkN=((nSA/saChunkSize+1)/P.runThreadN+1)*P.runThreadN;//ensure saChunkN is divisible by P.runThreadN @@ -254,44 +354,393 @@ void Genome::genomeGenerate() { indPrefStart[saChunkN]=indPrefN+1; indPrefChunkCount[saChunkN-1]=chunkSize1; + uint* saChunkStart = new uint [saChunkN+1]; + saChunkStart[0]=0; + for (uint iChunk=0; iChunklogMain << "Number of chunks: " << saChunkN <<"; chunks size limit: " << saChunkSize*8 <<" bytes\n" <logMain << timeMonthDayTime(rawTime) <<" ... sorting Suffix Array chunks and saving them to disk...\n" <logStdOut << timeMonthDayTime(rawTime) <<" ... sorting Suffix Array chunks and saving them to disk...\n" <=indPrefStart[iChunk] && p1logMain << timeMonthDayTime(rawTime) <<" ... sorting Suffix Array chunks...\n" <logStdOut << timeMonthDayTime(rawTime) <<" ... sorting Suffix Array chunks...\n" <nG1alloc ? P.limitGenomeGenerateRAM-nG1alloc : 0; + const uint64 saAllChunkScatterBytes=nSA*sizeof(uint) + (uint64) P.runThreadN*indPrefN*2*sizeof(uint) + (indPrefN+1)*sizeof(uint) + (saChunkN+1)*sizeof(uint); + const uint saChunkBatchSize=max(saChunkSize*(uint) P.runThreadN,saChunkSize); + uint saChunkBatchN=0; + if (P.runThreadN>=16 && saChunkN>1) { + uint iChunkBatch=0; + while (iChunkBatch=16 && saChunkN>1; + const uint64 saRetainedChunkBytes=nSA*sizeof(uint); + const uint64 saRamPackPeakBytes=saRetainedChunkBytes+SApass1.lengthByte; + const uint64 saRamPeakBytes=max(saAllChunkScatterBytes,saRamPackPeakBytes); + const uint64 saRamHeadroomBytes=max(saRamPeakBytes/10,(uint64) 2000000000LLU); + const uint64 saRamRequiredBytes=saRamPeakBytes+saRamHeadroomBytes; + const uint64 systemAvailableBytes=systemAvailableMemoryBytes(); + const bool saRamLimitOK=saAvailableBytes>=saRamRequiredBytes; + const bool saRamSystemOK=systemAvailableBytes>0 && systemAvailableBytes>=saRamRequiredBytes; + bool saChunksInMemoryActive=saRamLimitOK && saRamSystemOK; + P.inOut->logMain << "SA chunk all-scatter available bytes: " << saAvailableBytes << "; required temporary bytes: " << saAllChunkScatterBytes << "\n" <logMain << "SA chunk estimated batched-fill batches: " << saChunkBatchN << "\n" <logMain << "SA chunk fill strategy: " << (saChunkBatchFill ? "batched" : "per-chunk") << "\n" <logMain << "SA chunk sort prefix length: " << indPrefLen << "\n" <logMain << "SA chunk sort granularity: " << (saChunkBatchFill ? "prefix-bin" : "chunk") << "\n" <logMain << "SA chunk retained bytes: " << saRetainedChunkBytes << "; RAM peak bytes: " << saRamPeakBytes << "; RAM headroom bytes: " << saRamHeadroomBytes << "\n" <logMain << "SA chunk system available bytes: " << systemAvailableBytes << "; limit available bytes: " << saAvailableBytes << "\n" <saChunkMaxCount) { + saChunkMaxCount=indPrefChunkCount[iChunk]; + saChunkMaxIndex=iChunk; + }; + }; + P.inOut->logMain << "SA sort profile: enabled\n" <logMain << "SA sort profile prefix bins: total=" << indPrefN + << "; empty=" << indPrefEmptyN + << "; nonempty=" << indPrefNonEmptyN + << "; singleton=" << indPrefSingletonN + << "; sortable=" << indPrefSortableN + << "; max_bin=" << indPrefMaxCount + << "; max_prefix=" << indPrefMaxIndex << "\n" <logMain << "SA sort profile chunks: count=" << saChunkN + << "; min_indices=" << saChunkMinCount + << "; max_indices=" << saChunkMaxCount + << "; max_chunk=" << saChunkMaxIndex + << "; mean_indices=" << (saChunkN>0 ? nSA/saChunkN : 0) << "\n" <logMain << "SA chunk storage fallback: disk-backed after RAM allocation failure\n" <logMain << "SA chunk storage strategy: " << (saChunksInMemoryActive ? "RAM-retained" : "disk-backed") << "\n" <logMain, EXIT_CODE_INPUT_FILES, P); + }; + + uint* saThreadPrefCount=new uint [(uint) P.runThreadN*batchPrefN]; + memset(saThreadPrefCount,0,sizeof(saThreadPrefCount[0])*(uint) P.runThreadN*batchPrefN); + + double saProfileTime0=omp_get_wtime(); + #pragma omp parallel num_threads(P.runThreadN) + { + int tid=omp_get_thread_num(); + uint* threadPrefCount=saThreadPrefCount+(uint) tid*batchPrefN; + #pragma omp for schedule(static) + for (uint iScan=0;iScan=batchPrefStart && p1logMain, EXIT_CODE_INPUT_FILES, P); + }; + }; + + saProfileTime0=omp_get_wtime(); + #pragma omp parallel num_threads(P.runThreadN) + { + int tid=omp_get_thread_num(); + uint* threadPrefStart=saThreadPrefStart+(uint) tid*batchPrefN; + #pragma omp for schedule(static) + for (uint iScan=0;iScan=batchPrefStart && p1=(uint) sizeof(uint) && globalL>1; + const uint saSubBinExtraLen=3; + const uint saSubBinN=216; // 6^3 + const uint saSubBinMinN=1000000; + vector saSortRange; + vector saSortRangeN; + vector saSortRangeSkipFirstWord; + saSortRange.reserve(batchPrefN); + saSortRangeN.reserve(batchPrefN); + saSortRangeSkipFirstWord.reserve(batchPrefN); + for (uint iPref=0; iPref < batchPrefN; iPref++) { + uint* saPref=saBatch+saPrefStart[iPref]; + uint saPrefN=saPrefStart[iPref+1]-saPrefStart[iPref]; + if (saPrefN>1) { + bool skipFirstWord=saCanSkipFirstWord && !funSAsortPrefixFirstWordHas5(batchPrefStart+iPref,indPrefLen); + bool subBinned=false; + if (skipFirstWord && saPrefN>=saSubBinMinN) { + uint saSubBinCount[saSubBinN]; + memset(saSubBinCount,0,sizeof(saSubBinCount)); + bool subBinOK=true; + for (uint ii=0; ii0) { + saSubBinNonEmptyN++; + }; + saSubBinStart[ii+1]=saSubBinStart[ii]+saSubBinCount[ii]; + }; + }; + + uint* saPrefTmp=NULL; + if (saSubBinNonEmptyN>1) { + saPrefTmp=new (nothrow) uint [saPrefN]; + }; + if (saPrefTmp!=NULL) { + uint saSubBinCursor[saSubBinN]; + memcpy(saSubBinCursor,saSubBinStart,sizeof(saSubBinCursor)); + for (uint ii=0; ii1) { + saSortRange.push_back(saPref+saSubBinStart[iSub]); + saSortRangeN.push_back(saSubN); + saSortRangeSkipFirstWord.push_back(1); + saSubBinSortRangeN++; + }; + }; + if (saSortProfile) { + saProfileSubBinPrefixN++; + saProfileSubBinRangeN += saSubBinSortRangeN; + saProfileSubBinIndexN += saPrefN; + }; + subBinned=true; + }; + }; + if (!subBinned) { + saSortRange.push_back(saPref); + saSortRangeN.push_back(saPrefN); + saSortRangeSkipFirstWord.push_back(skipFirstWord ? 1 : 0); + }; + }; + }; + + #pragma omp parallel for num_threads(P.runThreadN) schedule(dynamic,16) reduction(+:saProfileSkipFirstWordRangeN,saProfileSkipFirstWordIndexN) + for (int iRange=0; iRange < (int) saSortRangeN.size(); iRange++) { + uint* saRange=saSortRange[iRange]; + uint saRangeN=saSortRangeN[iRange]; + if (saSortRangeSkipFirstWord[iRange]) { + qsort(saRange,saRangeN,sizeof(saRange[0]),funCompareSuffixesSkipFirstWord); + if (saSortProfile) { + saProfileSkipFirstWordRangeN++; + saProfileSkipFirstWordIndexN += saRangeN; + }; + } else { + qsort(saRange,saRangeN,sizeof(saRange[0]),funCompareSuffixes); + }; + }; + if (saSortProfile) { + saProfileSortSeconds += omp_get_wtime()-saProfileTime0; + }; + + saProfileTime0=omp_get_wtime(); + #pragma omp parallel for num_threads(P.runThreadN) schedule(dynamic,1) + for (int iChunk=0; iChunk < (int) batchChunkN; iChunk++) { + uint iChunkAbs=iChunkBatch+(uint) iChunk; + uint* saChunk=saBatch+saBatchStart[iChunk]; + for (uint ii=0;ii=indPrefStart[iChunk] && p11) { + qsort(saChunk,indPrefChunkCount[iChunk],sizeof(saChunk[0]),funCompareSuffixes); + }; + if (saSortProfile) { + saProfileSortSeconds += omp_get_wtime()-saProfileTime0; + }; + saProfileTime0=omp_get_wtime(); + for (uint ii=0;iilogMain << "SA sort profile comparator: skip_first_word_ranges=" << saProfileSkipFirstWordRangeN + << "; skip_first_word_indices=" << saProfileSkipFirstWordIndexN << "\n" <logMain << "SA sort profile sub-bins: prefixes=" << saProfileSubBinPrefixN + << "; ranges=" << saProfileSubBinRangeN + << "; indices=" << saProfileSubBinIndexN << "\n" <logMain << timeMonthDayTime(rawTime) <<" ... loading chunks from disk, packing SA...\n" <logStdOut << timeMonthDayTime(rawTime) <<" ... loading chunks from disk, packing SA...\n" <logMain << timeMonthDayTime(rawTime) << (saChunksInMemoryActive ? " ... packing SA from RAM chunks...\n" : " ... loading chunks from disk, packing SA...\n") <logStdOut << timeMonthDayTime(rawTime) << (saChunksInMemoryActive ? " ... packing SA from RAM chunks...\n" : " ... loading chunks from disk, packing SA...\n") <