diff --git a/.gitignore b/.gitignore index 3094469..51855af 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,2 @@ .DS_Store -build +build* diff --git a/CMakeLists.txt b/CMakeLists.txt index a4c624c..3c4ed9c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -106,6 +106,11 @@ if(SSHASH_BUILD_EXECUTABLES) sshash_static ) + add_executable(test_minimizer test/test_minimizer.cpp) + target_link_libraries(test_minimizer + sshash_static + ) + add_executable(check test/check.cpp) target_link_libraries(check sshash_static diff --git a/README.md b/README.md index 77101b1..d1ab97c 100644 --- a/README.md +++ b/README.md @@ -204,20 +204,14 @@ if your queries are to be read from a (multi-line) FASTA file. ### Example 3 - ./sshash build -i ../data/unitigs_stitched/salmonella_100_k31_ust.fa.gz -k 31 -m 13 --canonical -o salmonella_100.canon.sshash - -This example builds a dictionary from the input file `../data/unitigs_stitched/salmonella_100_k31_ust.fa.gz` (same used in Example 2), with k = 31, m = 13, and with the canonical parsing modality (option `--canonical`). The dictionary is serialized on disk to the file `salmonella_100.canon.sshash`. - -The "canonical" version of the dictionary offers more speed for only a little space increase, especially under low-hit workloads -- when the majority of k-mers are not found in the dictionary. (For all details, refer to the paper.) - -Below a comparison between the dictionary built in Example 2 (not canonical) -and the one just built (Example 3, canonical). - ./sshash query -i salmonella_100.sshash -q ../data/queries/SRR5833294.10K.fastq.gz - ./sshash query -i salmonella_100.canon.sshash -q ../data/queries/SRR5833294.10K.fastq.gz +A k-mer and its reverse complement always share a bucket, so a lookup costs a +single bucket probe and reports, in `kmer_orientation`, which of the two was +found. This is the only indexing modality: there is no `--canonical` flag to +choose, and no space-versus-speed trade-off to make. -Both queries should originate the following report (reported here for reference): +The query above should originate the following report (reported here for reference): ==== query report: num_kmers = 460000 @@ -225,9 +219,6 @@ Both queries should originate the following report (reported here for reference) num_searches = 42/46 (91.3043%) num_extensions = 4/46 (8.69565%) -The canonical dictionary can be twice as fast as the regular dictionary -for low-hit workloads, even on this tiny example, for only +0.3 bits/k-mer. - ### Example 4 ./sshash permute -i ../data/unitigs_stitched/with_weights/ecoli_sakai.ust.k31.fa.gz -k 31 -o ecoli_sakai.permuted.fa diff --git a/benchmarks/README.md b/benchmarks/README.md index 5e65e01..c17b87e 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -14,16 +14,42 @@ To run the benchmarks, from within the `build` directory, run where `` should be replaced by a suitable basename, e.g., the current date. -These are the results obtained on 21/01/26 (see logs [here](results-21-01-26)) +These are the results obtained on 28/08/26 (see logs [here](results-28-08-26)) on a machine equipped with an AMD Ryzen Threadripper PRO 7985WX processor clocked at 5.40GHz. The code was compiled with `gcc` 13.3.0. - The indexes were build with a max RAM usage of 16 GB and 64 threads. Queries were run using one thread, instead. -![](results-21-01-26/results.png) - -The results can be exported to CSV format with - - python3 ../script/print_csv.py ../benchmarks/results-21-01-26/k31 - python3 ../script/print_csv.py ../benchmarks/results-21-01-26/k63 +| k | Collection | m | Space (bits/kmer) | Space (total GB) | Building time (m:ss) | Positive random lookup (µs/kmer) | Negative random lookup (µs/kmer) | Random Access (µs/kmer) | Streaming Lookup high-hit (ns/kmer) | +|----|------------|:--:|:-----------------:|:----------------:|:--------------------:|:--------------------------------:|:--------------------------------:|:-----------------------:|:-----------------------------------:| +|| +| 31 | Cod | 20 | 7.94 | 0.50 | 0:27 | 0.43 | 0.36 | 0.27 | 24 | +| | Kestrel | 20 | 7.52 | 1.08 | 1:00 | 0.44 | 0.39 | 0.28 | 39 | +| | Human | 21 | 8.77 | 2.75 | 2:35 | 0.56 | 0.43 | 0.35 | 76 | +| | NCBI-v | 19 | 7.43 | 0.35 | 0:14 | 0.41 | 0.35 | 0.26 | 26 | +| | SE | 21 | 10.24 | 1.14 | 0:54 | 0.62 | 0.40 | 0.36 | 163 | +| | HPRC | 21 | 10.64 | 4.94 | 4:36 | 0.70 | 0.46 | 0.53 | 84 | +|| +| 63 | Cod | 24 | 4.62 | 0.32 | 0:24 | 0.57 | 0.45 | 0.29 | 37 | +| | Kestrel | 24 | 3.84 | 0.55 | 0:19 | 0.54 | 0.48 | 0.32 | 44 | +| | Human | 25 | 4.92 | 1.70 | 1:13 | 0.66 | 0.51 | 0.36 | 110 | +| | NCBI-v | 23 | 4.06 | 0.21 | 0:07 | 0.53 | 0.44 | 0.28 | 30 | +| | SE | 31 | 7.42 | 1.41 | 1:20 | 0.97 | 0.50 | 0.40 | 278 | +| | HPRC | 31 | 7.63 | 5.65 | 4:46 | 0.94 | 0.58 | 0.62 | 137 | + +The results of a run, for both k=31 and k=63 at once, can be exported to CSV +format (every collected quantity) or to a markdown table (the main quantities) +with + + python3 ../benchmarks/print_table.py + python3 ../benchmarks/print_table.py --md + +where `` is the directory written by the scripts above (e.g., +`results-21-01-26`), containing the `k31` and `k63` subdirectories. + +Note that the scripts now produce a single set of result files per `k` +(`build.json`, `bench.json`, and `streaming-queries.json`), since the +regular/canonical distinction is gone: there is one indexing modality only. +Result directories archived before this change instead hold a `regular-` and a +`canon-` file for each of those; to re-read them, use the `print_csv.py` +script from the corresponding commit. diff --git a/benchmarks/print_csv.py b/benchmarks/print_csv.py deleted file mode 100644 index 6e60816..0000000 --- a/benchmarks/print_csv.py +++ /dev/null @@ -1,175 +0,0 @@ -#!/usr/bin/env python3 - -import sys -import json -import os -from statistics import mean, StatisticsError -import math - -def format_time(microseconds): - seconds = microseconds / 1_000_000 - minutes = int(seconds // 60) - seconds = int(seconds % 60) - return f"{minutes}:{seconds:02d}" - -def parse_build_file(path, canonical_flag): - """Parse build JSONL file.""" - results = [] - with open(path) as f: - for line in f: - line = line.strip() - if not line: - continue - try: - d = json.loads(line) - except json.JSONDecodeError: - print(f"Skipping invalid JSON line in {path}", file=sys.stderr) - continue - - num_kmers = int(d["num_kmers"]) - index_bytes = int(d["index_size_in_bytes"]) - build_time_us = int(d["total_build_time_in_microsec"]) - - bits_per_kmer = (index_bytes * 8) / num_kmers - gb = index_bytes / 1e9 - build_time_fmt = format_time(build_time_us) - - fname = os.path.basename(d["input_filename"]) - collection = fname.split(".")[0].capitalize() - k = d["k"] - - results.append({ - "k": k, - "Collection": collection, - "m": d["m"], - "canonical": "yes" if canonical_flag else "no", - "bits_per_kmer": f"{bits_per_kmer:.2f}", - "total_GB": f"{gb:.2f}", - "build_time": build_time_fmt - }) - return results - -def parse_bench_file(path, canonical_flag): - """Parse benchmark JSONL file and average per collection.""" - lookup_data = {} - with open(path) as f: - for line in f: - line = line.strip() - if not line: - continue - try: - d = json.loads(line) - except json.JSONDecodeError: - print(f"Skipping invalid JSON line in {path}", file=sys.stderr) - continue - - fname = os.path.basename(d["index_filename"]) - collection = fname.split(".")[0].capitalize() - m = d["m"] - k = d["k"] - canonical = "yes" if canonical_flag else "no" - - key = (collection, m, canonical) - entry = lookup_data.setdefault(key, { - "k": k, - "pos": [], "neg": [], "access": [], "iter": [] - }) - entry["pos"].append(float(d["positive lookup (avg_nanosec_per_kmer)"])) - entry["neg"].append(float(d["negative lookup (avg_nanosec_per_kmer)"])) - entry["access"].append(float(d["access (avg_nanosec_per_kmer)"])) - entry["iter"].append(float(d["iterator (avg_nanosec_per_kmer)"])) - - # average the results - for k, v in lookup_data.items(): - try: - lookup_data[k] = { - "k": v["k"], - "pos": f"{mean(v['pos'])/1000:.2f}", - "neg": f"{mean(v['neg'])/1000:.2f}", - "access": f"{mean(v['access'])/1000:.2f}", - "iter": f"{mean(v['iter']):.2f}", - } - except StatisticsError: - lookup_data[k] = {"k": v["k"], "pos": "NA", "neg": "NA", "access": "NA", "iter": "NA"} - return lookup_data - - -def parse_streaming_file(path, canonical_flag): - """Parse streaming queries JSON file.""" - stream_data = {} - if not os.path.exists(path): - return stream_data - - with open(path) as f: - for line in f: - line = line.strip() - if not line: - continue - try: - d = json.loads(line) - except json.JSONDecodeError: - print(f"Skipping invalid JSON line in {path}", file=sys.stderr) - continue - - fname = os.path.basename(d["index_filename"]) - collection = fname.split(".")[0].capitalize() - canonical = "yes" if canonical_flag else "no" - - key = (collection, canonical) - num_kmers = int(d["num_kmers"]) - num_pos = int(d["num_positive_kmers"]) - num_ext = int(d["num_extensions"]) - elapsed_ms = int(d["elapsed_millisec"]) - - ns_per_kmer = int(math.ceil(elapsed_ms * 1e6 / num_kmers)) - hit_rate = (num_pos / num_kmers) * 100 if num_kmers else 0 - extension_rate = (num_ext / num_pos) * 100 if num_pos else 0 - - stream_data[key] = { - "ns_per_kmer": f"{ns_per_kmer}", - "hit_rate": f"{hit_rate:.2f}", - "extension_rate": f"{extension_rate:.2f}" - } - return stream_data - - -def main(): - if len(sys.argv) != 2: - print("Usage: print.py input_dir", file=sys.stderr) - sys.exit(1) - - input_dir = sys.argv[1] - reg_build_path = input_dir + "/regular-build.json" - canon_build_path = input_dir + "/canon-build.json" - reg_bench_path = input_dir + "/regular-bench.json" - canon_bench_path = input_dir + "/canon-bench.json" - reg_stream_path = input_dir + "/regular-streaming-queries-high-hit.json" - canon_stream_path = input_dir + "/canon-streaming-queries-high-hit.json" - - reg_build = parse_build_file(reg_build_path, False) - canon_build = parse_build_file(canon_build_path, True) - reg_bench = parse_bench_file(reg_bench_path, False) - canon_bench = parse_bench_file(canon_bench_path, True) - reg_stream = parse_streaming_file(reg_stream_path, False) - canon_stream = parse_streaming_file(canon_stream_path, True) - - # merge everything - all_builds = reg_build + canon_build - lookup_all = {**reg_bench, **canon_bench} - stream_all = {**reg_stream, **canon_stream} - - # CSV header - print("k,Collection,m,canonical,bits_per_kmer,total_GB,build_time,positive_lookup_ns,negative_lookup_ns,access_ns,iteration_ns,ns_per_kmer,hit_rate,extension_rate") - - for r in sorted(all_builds, key=lambda x: (int(x["k"]), x["Collection"], x["canonical"])): - lookup = lookup_all.get( - (r["Collection"], r["m"], r["canonical"]), # key - {"pos": "NA", "neg": "NA", "access": "NA", "iter": "NA", "k": r["k"]}) - stream = stream_all.get( - (r["Collection"], r["canonical"]), # key - {"ns_per_kmer": "NA", "hit_rate": "NA", "extension_rate": "NA"}) - - print(f"{r['k']},{r['Collection']},{r['m']},{r['canonical']},{r['bits_per_kmer']},{r['total_GB']},{r['build_time']},{lookup['pos']},{lookup['neg']},{lookup['access']},{lookup['iter']},{stream['ns_per_kmer']},{stream['hit_rate']},{stream['extension_rate']}") - -if __name__ == "__main__": - main() diff --git a/benchmarks/print_table.py b/benchmarks/print_table.py new file mode 100644 index 0000000..b79fda8 --- /dev/null +++ b/benchmarks/print_table.py @@ -0,0 +1,214 @@ +#!/usr/bin/env python3 + +""" +Print the benchmark results of both k=31 and k=63 as one table. + +The input directory must contain the two subdirectories `k31` and `k63`, +each with the files written by the benchmark scripts: +`build.json`, `bench.json`, and `streaming-queries.json`. + +The default output is CSV, with every collected quantity; with `--md`, +a markdown table with the main quantities is printed instead. +""" + +import argparse +import json +import os +import sys +import math +from statistics import mean, StatisticsError + +K_VALUES = [31, 63] + +# fixed dataset order and display names +DATASETS = ["cod", "kestrel", "human", "ncbi-virus", "se", "hprc"] +DISPLAY_NAME = { + "cod": "Cod", + "kestrel": "Kestrel", + "human": "Human", + "ncbi-virus": "NCBI-v", + "se": "SE", + "hprc": "HPRC", +} + +def format_time(microseconds): + seconds = microseconds / 1_000_000 + minutes = int(seconds // 60) + seconds = int(seconds % 60) + return f"{minutes}:{seconds:02d}" + +def dataset_of(filename): + return os.path.basename(filename).split(".")[0] + +def parse_build_file(path): + """Parse build JSONL file: one record per dataset.""" + results = {} + with open(path) as f: + for line in f: + line = line.strip() + if not line: + continue + try: + d = json.loads(line) + except json.JSONDecodeError: + print(f"Skipping invalid JSON line in {path}", file=sys.stderr) + continue + + num_kmers = int(d["num_kmers"]) + index_bytes = int(d["index_size_in_bytes"]) + build_time_us = int(d["total_build_time_in_microsec"]) + + results[dataset_of(d["input_filename"])] = { + "m": d["m"], + "bits_per_kmer": f"{(index_bytes * 8) / num_kmers:.2f}", + "total_GB": f"{index_bytes / 1e9:.2f}", + "build_time": format_time(build_time_us), + } + return results + +def parse_bench_file(path): + """Parse benchmark JSONL file and average the runs per dataset.""" + runs = {} + with open(path) as f: + for line in f: + line = line.strip() + if not line: + continue + try: + d = json.loads(line) + except json.JSONDecodeError: + print(f"Skipping invalid JSON line in {path}", file=sys.stderr) + continue + + entry = runs.setdefault(dataset_of(d["index_filename"]), + {"pos": [], "neg": [], "access": [], "iter": []}) + entry["pos"].append(float(d["positive lookup (avg_nanosec_per_kmer)"])) + entry["neg"].append(float(d["negative lookup (avg_nanosec_per_kmer)"])) + entry["access"].append(float(d["access (avg_nanosec_per_kmer)"])) + entry["iter"].append(float(d["iterator (avg_nanosec_per_kmer)"])) + + results = {} + for ds, v in runs.items(): + try: + results[ds] = { + "pos_us": f"{mean(v['pos']) / 1000:.2f}", + "neg_us": f"{mean(v['neg']) / 1000:.2f}", + "access_us": f"{mean(v['access']) / 1000:.2f}", + "iter_ns": f"{mean(v['iter']):.2f}", + } + except StatisticsError: + results[ds] = {"pos_us": "NA", "neg_us": "NA", "access_us": "NA", "iter_ns": "NA"} + return results + +def parse_streaming_file(path): + """Parse streaming queries JSONL file: one record per dataset.""" + results = {} + if not os.path.exists(path): + return results + + with open(path) as f: + for line in f: + line = line.strip() + if not line: + continue + try: + d = json.loads(line) + except json.JSONDecodeError: + print(f"Skipping invalid JSON line in {path}", file=sys.stderr) + continue + + num_kmers = int(d["num_kmers"]) + num_pos = int(d["num_positive_kmers"]) + num_ext = int(d["num_extensions"]) + elapsed_ms = int(d["elapsed_millisec"]) + + results[dataset_of(d["index_filename"])] = { + "ns_per_kmer": f"{int(math.ceil(elapsed_ms * 1e6 / num_kmers))}", + "hit_rate": f"{(num_pos / num_kmers) * 100 if num_kmers else 0:.2f}", + "extension_rate": f"{(num_ext / num_pos) * 100 if num_pos else 0:.2f}", + } + return results + +def collect(input_dir): + """One row per (k, dataset), in fixed order.""" + rows = [] + for k in K_VALUES: + k_dir = os.path.join(input_dir, f"k{k}") + if not os.path.isdir(k_dir): + print(f"Error: expected subdirectory '{k_dir}'", file=sys.stderr) + sys.exit(1) + build = parse_build_file(os.path.join(k_dir, "build.json")) + bench = parse_bench_file(os.path.join(k_dir, "bench.json")) + stream = parse_streaming_file(os.path.join(k_dir, "streaming-queries.json")) + first_row_of_k = len(rows) + for ds in DATASETS: + if ds not in build: + continue + row = {"k": str(k), "Collection": DISPLAY_NAME[ds], "first_of_k": False} + row.update(build[ds]) + row.update(bench.get(ds, {"pos_us": "NA", "neg_us": "NA", + "access_us": "NA", "iter_ns": "NA"})) + row.update(stream.get(ds, {"ns_per_kmer": "NA", "hit_rate": "NA", + "extension_rate": "NA"})) + rows.append(row) + if len(rows) > first_row_of_k: + rows[first_row_of_k]["first_of_k"] = True + return rows + +def print_csv(rows): + print("k,Collection,m,bits_per_kmer,total_GB,build_time," + "positive_lookup_us,negative_lookup_us,access_us,iteration_ns," + "streaming_ns_per_kmer,hit_rate,extension_rate") + for r in rows: + print(f"{r['k']},{r['Collection']},{r['m']},{r['bits_per_kmer']},{r['total_GB']}," + f"{r['build_time']},{r['pos_us']},{r['neg_us']},{r['access_us']},{r['iter_ns']}," + f"{r['ns_per_kmer']},{r['hit_rate']},{r['extension_rate']}") + +MD_COLUMNS = [ + # (header, row key, alignment: 'left' or 'center') + ("k", "k", "left"), + ("Collection", "Collection", "left"), + ("m", "m", "center"), + ("Space (bits/kmer)", "bits_per_kmer", "center"), + ("Space (total GB)", "total_GB", "center"), + ("Building time (m:ss)", "build_time", "center"), + ("Positive random lookup (µs/kmer)", "pos_us", "center"), + ("Negative random lookup (µs/kmer)", "neg_us", "center"), + ("Random Access (µs/kmer)", "access_us", "center"), + ("Streaming Lookup high-hit (ns/kmer)", "ns_per_kmer", "center"), +] + +def print_md(rows): + """The k value is shown on the first row of its block only, + and each block is preceded by an empty '||' row.""" + cells = [[str(r[key]) if key != "k" or r["first_of_k"] else "" + for _, key, _ in MD_COLUMNS] for r in rows] + widths = [max(len(header), *(len(row[i]) for row in cells)) + for i, (header, _, _) in enumerate(MD_COLUMNS)] + + def line(values): + return "| " + " | ".join(v.ljust(w) for v, w in zip(values, widths)) + " |" + + print(line([header for header, _, _ in MD_COLUMNS])) + print("|" + "|".join(f":{'-' * w}:" if align == "center" else "-" * (w + 2) + for w, (_, _, align) in zip(widths, MD_COLUMNS)) + "|") + for r, row_cells in zip(rows, cells): + if r["first_of_k"]: + print("||") + print(line(row_cells)) + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("input_dir", help="directory containing the k31 and k63 subdirectories") + parser.add_argument("--md", action="store_true", + help="print a markdown table instead of CSV") + args = parser.parse_args() + + rows = collect(args.input_dir) + if args.md: + print_md(rows) + else: + print_csv(rows) + +if __name__ == "__main__": + main() diff --git a/benchmarks/results-28-08-26/k31/bench.json b/benchmarks/results-28-08-26/k31/bench.json new file mode 100644 index 0000000..ffdc5ec --- /dev/null +++ b/benchmarks/results-28-08-26/k31/bench.json @@ -0,0 +1,18 @@ +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "436.318757", "negative lookup (avg_nanosec_per_kmer)": "362.749596", "access (avg_nanosec_per_kmer)": "272.585695", "iterator (avg_nanosec_per_kmer)": "2.482751"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "433.819361", "negative lookup (avg_nanosec_per_kmer)": "362.259503", "access (avg_nanosec_per_kmer)": "274.144289", "iterator (avg_nanosec_per_kmer)": "2.456011"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "433.500928", "negative lookup (avg_nanosec_per_kmer)": "364.396456", "access (avg_nanosec_per_kmer)": "273.884328", "iterator (avg_nanosec_per_kmer)": "2.569083"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "439.680768", "negative lookup (avg_nanosec_per_kmer)": "395.603409", "access (avg_nanosec_per_kmer)": "276.121400", "iterator (avg_nanosec_per_kmer)": "2.499332"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "444.361067", "negative lookup (avg_nanosec_per_kmer)": "394.031758", "access (avg_nanosec_per_kmer)": "274.707766", "iterator (avg_nanosec_per_kmer)": "2.524888"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "440.037087", "negative lookup (avg_nanosec_per_kmer)": "393.136266", "access (avg_nanosec_per_kmer)": "275.896044", "iterator (avg_nanosec_per_kmer)": "2.538501"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "562.807039", "negative lookup (avg_nanosec_per_kmer)": "431.088550", "access (avg_nanosec_per_kmer)": "347.198096", "iterator (avg_nanosec_per_kmer)": "2.465738"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "555.812279", "negative lookup (avg_nanosec_per_kmer)": "422.575776", "access (avg_nanosec_per_kmer)": "348.303889", "iterator (avg_nanosec_per_kmer)": "2.496474"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "561.630325", "negative lookup (avg_nanosec_per_kmer)": "424.718263", "access (avg_nanosec_per_kmer)": "347.852823", "iterator (avg_nanosec_per_kmer)": "2.547386"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash", "k": "31", "m": "19", "positive lookup (avg_nanosec_per_kmer)": "412.510594", "negative lookup (avg_nanosec_per_kmer)": "348.734417", "access (avg_nanosec_per_kmer)": "259.806846", "iterator (avg_nanosec_per_kmer)": "2.500924"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash", "k": "31", "m": "19", "positive lookup (avg_nanosec_per_kmer)": "407.386907", "negative lookup (avg_nanosec_per_kmer)": "353.408523", "access (avg_nanosec_per_kmer)": "260.417569", "iterator (avg_nanosec_per_kmer)": "2.459740"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash", "k": "31", "m": "19", "positive lookup (avg_nanosec_per_kmer)": "408.124222", "negative lookup (avg_nanosec_per_kmer)": "346.954378", "access (avg_nanosec_per_kmer)": "260.970479", "iterator (avg_nanosec_per_kmer)": "2.504823"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "625.169634", "negative lookup (avg_nanosec_per_kmer)": "396.027242", "access (avg_nanosec_per_kmer)": "357.092978", "iterator (avg_nanosec_per_kmer)": "2.475710"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "614.495341", "negative lookup (avg_nanosec_per_kmer)": "396.575080", "access (avg_nanosec_per_kmer)": "354.582941", "iterator (avg_nanosec_per_kmer)": "2.414229"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "623.312167", "negative lookup (avg_nanosec_per_kmer)": "400.475688", "access (avg_nanosec_per_kmer)": "357.971520", "iterator (avg_nanosec_per_kmer)": "2.508970"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "696.768429", "negative lookup (avg_nanosec_per_kmer)": "459.656809", "access (avg_nanosec_per_kmer)": "528.457965", "iterator (avg_nanosec_per_kmer)": "2.471353"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "699.498620", "negative lookup (avg_nanosec_per_kmer)": "461.383063", "access (avg_nanosec_per_kmer)": "527.811735", "iterator (avg_nanosec_per_kmer)": "2.509292"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "696.245840", "negative lookup (avg_nanosec_per_kmer)": "458.446995", "access (avg_nanosec_per_kmer)": "521.370161", "iterator (avg_nanosec_per_kmer)": "2.440193"} diff --git a/benchmarks/results-28-08-26/k31/bench.log b/benchmarks/results-28-08-26/k31/bench.log new file mode 100644 index 0000000..5259139 --- /dev/null +++ b/benchmarks/results-28-08-26/k31/bench.log @@ -0,0 +1,90 @@ +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 436.319 +negative lookup (avg_nanosec_per_kmer) 362.75 +access (avg_nanosec_per_kmer) = 272.586 +iterator (avg_nanosec_per_kmer) = 2.48275 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 433.819 +negative lookup (avg_nanosec_per_kmer) 362.26 +access (avg_nanosec_per_kmer) = 274.144 +iterator (avg_nanosec_per_kmer) = 2.45601 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 433.501 +negative lookup (avg_nanosec_per_kmer) 364.396 +access (avg_nanosec_per_kmer) = 273.884 +iterator (avg_nanosec_per_kmer) = 2.56908 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 439.681 +negative lookup (avg_nanosec_per_kmer) 395.603 +access (avg_nanosec_per_kmer) = 276.121 +iterator (avg_nanosec_per_kmer) = 2.49933 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 444.361 +negative lookup (avg_nanosec_per_kmer) 394.032 +access (avg_nanosec_per_kmer) = 274.708 +iterator (avg_nanosec_per_kmer) = 2.52489 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 440.037 +negative lookup (avg_nanosec_per_kmer) 393.136 +access (avg_nanosec_per_kmer) = 275.896 +iterator (avg_nanosec_per_kmer) = 2.5385 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 562.807 +negative lookup (avg_nanosec_per_kmer) 431.089 +access (avg_nanosec_per_kmer) = 347.198 +iterator (avg_nanosec_per_kmer) = 2.46574 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 555.812 +negative lookup (avg_nanosec_per_kmer) 422.576 +access (avg_nanosec_per_kmer) = 348.304 +iterator (avg_nanosec_per_kmer) = 2.49647 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 561.63 +negative lookup (avg_nanosec_per_kmer) 424.718 +access (avg_nanosec_per_kmer) = 347.853 +iterator (avg_nanosec_per_kmer) = 2.54739 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 412.511 +negative lookup (avg_nanosec_per_kmer) 348.734 +access (avg_nanosec_per_kmer) = 259.807 +iterator (avg_nanosec_per_kmer) = 2.50092 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 407.387 +negative lookup (avg_nanosec_per_kmer) 353.409 +access (avg_nanosec_per_kmer) = 260.418 +iterator (avg_nanosec_per_kmer) = 2.45974 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 408.124 +negative lookup (avg_nanosec_per_kmer) 346.954 +access (avg_nanosec_per_kmer) = 260.97 +iterator (avg_nanosec_per_kmer) = 2.50482 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 625.17 +negative lookup (avg_nanosec_per_kmer) 396.027 +access (avg_nanosec_per_kmer) = 357.093 +iterator (avg_nanosec_per_kmer) = 2.47571 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 614.495 +negative lookup (avg_nanosec_per_kmer) 396.575 +access (avg_nanosec_per_kmer) = 354.583 +iterator (avg_nanosec_per_kmer) = 2.41423 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 623.312 +negative lookup (avg_nanosec_per_kmer) 400.476 +access (avg_nanosec_per_kmer) = 357.972 +iterator (avg_nanosec_per_kmer) = 2.50897 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 696.768 +negative lookup (avg_nanosec_per_kmer) 459.657 +access (avg_nanosec_per_kmer) = 528.458 +iterator (avg_nanosec_per_kmer) = 2.47135 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 699.499 +negative lookup (avg_nanosec_per_kmer) 461.383 +access (avg_nanosec_per_kmer) = 527.812 +iterator (avg_nanosec_per_kmer) = 2.50929 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 696.246 +negative lookup (avg_nanosec_per_kmer) 458.447 +access (avg_nanosec_per_kmer) = 521.37 +iterator (avg_nanosec_per_kmer) = 2.44019 diff --git a/benchmarks/results-28-08-26/k31/streaming-queries.json b/benchmarks/results-28-08-26/k31/streaming-queries.json new file mode 100644 index 0000000..d7cbe12 --- /dev/null +++ b/benchmarks/results-28-08-26/k31/streaming-queries.json @@ -0,0 +1,6 @@ +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz", "num_kmers": "163287360", "num_positive_kmers": "132860997", "num_negative_kmers": "30426363", "num_invalid_kmers": "0", "num_searches": "6576340", "num_extensions": "126284657", "num_memo_singleton": "193441", "num_memo_light": "6812848", "elapsed_millisec": "3770"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz", "num_kmers": "695737535", "num_positive_kmers": "525542891", "num_negative_kmers": "170183654", "num_invalid_kmers": "10990", "num_searches": "12437476", "num_extensions": "513105415", "num_memo_singleton": "997817", "num_memo_light": "34475911", "elapsed_millisec": "26537"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz", "num_kmers": "1569974986", "num_positive_kmers": "1437949378", "num_negative_kmers": "130996597", "num_invalid_kmers": "1029011", "num_searches": "100222623", "num_extensions": "1337726755", "num_memo_singleton": "1072506", "num_memo_light": "27289697", "elapsed_millisec": "119021"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz", "num_kmers": "14092875", "num_positive_kmers": "13983775", "num_negative_kmers": "108161", "num_invalid_kmers": "939", "num_searches": "590894", "num_extensions": "13392881", "num_memo_singleton": "433", "num_memo_light": "107455", "elapsed_millisec": "356"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz", "num_kmers": "789838196", "num_positive_kmers": "764882549", "num_negative_kmers": "24935381", "num_invalid_kmers": "20266", "num_searches": "218875709", "num_extensions": "546006840", "num_memo_singleton": "85669", "num_memo_light": "7940520", "elapsed_millisec": "128190"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz", "num_kmers": "1569974986", "num_positive_kmers": "1485223278", "num_negative_kmers": "83722697", "num_invalid_kmers": "1029011", "num_searches": "135823240", "num_extensions": "1349400038", "num_memo_singleton": "565997", "num_memo_light": "17514165", "elapsed_millisec": "131342"} diff --git a/benchmarks/results-28-08-26/k31/streaming-queries.log b/benchmarks/results-28-08-26/k31/streaming-queries.log new file mode 100644 index 0000000..b2b5e62 --- /dev/null +++ b/benchmarks/results-28-08-26/k31/streaming-queries.log @@ -0,0 +1,78 @@ +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz +2026-08-28 20:34:21: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz'... +2026-08-28 20:34:24: DONE +==== query report: +num_kmers = 163287360 +num_positive_kmers = 132860997 (81.3664%) +num_negative_kmers = 30426363 (18.6336%) +num_invalid_kmers = 0 (0%) +num_searches = 6576340/132860997 (4.94979%) +num_extensions = 126284657/132860997 (95.0502%) +num_memo_singleton = 193441 +num_memo_light = 6812848 +elapsed = 3.77 sec / 0.0628333 min / 23.0881 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz +2026-08-28 20:34:25: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz'... +2026-08-28 20:34:52: DONE +==== query report: +num_kmers = 695737535 +num_positive_kmers = 525542891 (75.5375%) +num_negative_kmers = 170183654 (24.4609%) +num_invalid_kmers = 10990 (0.00157962%) +num_searches = 12437476/525542891 (2.3666%) +num_extensions = 513105415/525542891 (97.6334%) +num_memo_singleton = 997817 +num_memo_light = 34475911 +elapsed = 26.537 sec / 0.442283 min / 38.1423 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz +2026-08-28 20:34:53: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz'... +2026-08-28 20:36:52: DONE +==== query report: +num_kmers = 1569974986 +num_positive_kmers = 1437949378 (91.5906%) +num_negative_kmers = 130996597 (8.34387%) +num_invalid_kmers = 1029011 (0.0655431%) +num_searches = 100222623/1437949378 (6.96983%) +num_extensions = 1337726755/1437949378 (93.0302%) +num_memo_singleton = 1072506 +num_memo_light = 27289697 +elapsed = 119.021 sec / 1.98368 min / 75.8108 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz +2026-08-28 20:36:52: performing queries from file '/mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz'... +2026-08-28 20:36:53: DONE +==== query report: +num_kmers = 14092875 +num_positive_kmers = 13983775 (99.2258%) +num_negative_kmers = 108161 (0.767487%) +num_invalid_kmers = 939 (0.00666294%) +num_searches = 590894/13983775 (4.22557%) +num_extensions = 13392881/13983775 (95.7744%) +num_memo_singleton = 433 +num_memo_light = 107455 +elapsed = 0.356 sec / 0.00593333 min / 25.261 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz +2026-08-28 20:36:53: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz'... +2026-08-28 20:39:01: DONE +==== query report: +num_kmers = 789838196 +num_positive_kmers = 764882549 (96.8404%) +num_negative_kmers = 24935381 (3.15702%) +num_invalid_kmers = 20266 (0.00256584%) +num_searches = 218875709/764882549 (28.6156%) +num_extensions = 546006840/764882549 (71.3844%) +num_memo_singleton = 85669 +num_memo_light = 7940520 +elapsed = 128.19 sec / 2.1365 min / 162.299 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz +2026-08-28 20:39:04: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz'... +2026-08-28 20:41:15: DONE +==== query report: +num_kmers = 1569974986 +num_positive_kmers = 1485223278 (94.6017%) +num_negative_kmers = 83722697 (5.33274%) +num_invalid_kmers = 1029011 (0.0655431%) +num_searches = 135823240/1485223278 (9.14497%) +num_extensions = 1349400038/1485223278 (90.855%) +num_memo_singleton = 565997 +num_memo_light = 17514165 +elapsed = 131.342 sec / 2.18903 min / 83.6587 ns/kmer diff --git a/benchmarks/results-28-08-26/k63/bench.json b/benchmarks/results-28-08-26/k63/bench.json new file mode 100644 index 0000000..49edb22 --- /dev/null +++ b/benchmarks/results-28-08-26/k63/bench.json @@ -0,0 +1,18 @@ +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "567.271179", "negative lookup (avg_nanosec_per_kmer)": "447.896700", "access (avg_nanosec_per_kmer)": "292.035089", "iterator (avg_nanosec_per_kmer)": "2.723525"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "567.158980", "negative lookup (avg_nanosec_per_kmer)": "442.781347", "access (avg_nanosec_per_kmer)": "290.059058", "iterator (avg_nanosec_per_kmer)": "2.716768"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "577.632402", "negative lookup (avg_nanosec_per_kmer)": "445.099827", "access (avg_nanosec_per_kmer)": "289.093225", "iterator (avg_nanosec_per_kmer)": "2.713838"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "536.517146", "negative lookup (avg_nanosec_per_kmer)": "479.223758", "access (avg_nanosec_per_kmer)": "324.438066", "iterator (avg_nanosec_per_kmer)": "2.757624"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "537.308816", "negative lookup (avg_nanosec_per_kmer)": "482.519123", "access (avg_nanosec_per_kmer)": "327.085464", "iterator (avg_nanosec_per_kmer)": "2.705321"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "543.875553", "negative lookup (avg_nanosec_per_kmer)": "484.284402", "access (avg_nanosec_per_kmer)": "321.921078", "iterator (avg_nanosec_per_kmer)": "2.708895"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash", "k": "63", "m": "25", "positive lookup (avg_nanosec_per_kmer)": "655.823771", "negative lookup (avg_nanosec_per_kmer)": "512.844700", "access (avg_nanosec_per_kmer)": "355.444918", "iterator (avg_nanosec_per_kmer)": "2.719615"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash", "k": "63", "m": "25", "positive lookup (avg_nanosec_per_kmer)": "657.564050", "negative lookup (avg_nanosec_per_kmer)": "514.922087", "access (avg_nanosec_per_kmer)": "356.297916", "iterator (avg_nanosec_per_kmer)": "2.773051"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash", "k": "63", "m": "25", "positive lookup (avg_nanosec_per_kmer)": "659.735498", "negative lookup (avg_nanosec_per_kmer)": "510.912399", "access (avg_nanosec_per_kmer)": "355.614307", "iterator (avg_nanosec_per_kmer)": "2.715916"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash", "k": "63", "m": "23", "positive lookup (avg_nanosec_per_kmer)": "531.164346", "negative lookup (avg_nanosec_per_kmer)": "434.532611", "access (avg_nanosec_per_kmer)": "281.460643", "iterator (avg_nanosec_per_kmer)": "2.733554"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash", "k": "63", "m": "23", "positive lookup (avg_nanosec_per_kmer)": "525.985610", "negative lookup (avg_nanosec_per_kmer)": "437.954528", "access (avg_nanosec_per_kmer)": "277.253213", "iterator (avg_nanosec_per_kmer)": "2.727797"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash", "k": "63", "m": "23", "positive lookup (avg_nanosec_per_kmer)": "531.091279", "negative lookup (avg_nanosec_per_kmer)": "434.230904", "access (avg_nanosec_per_kmer)": "283.331560", "iterator (avg_nanosec_per_kmer)": "2.727683"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "965.438332", "negative lookup (avg_nanosec_per_kmer)": "495.346946", "access (avg_nanosec_per_kmer)": "402.682535", "iterator (avg_nanosec_per_kmer)": "2.755412"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "969.421164", "negative lookup (avg_nanosec_per_kmer)": "494.920548", "access (avg_nanosec_per_kmer)": "404.439983", "iterator (avg_nanosec_per_kmer)": "2.749083"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "961.154818", "negative lookup (avg_nanosec_per_kmer)": "498.388122", "access (avg_nanosec_per_kmer)": "400.903357", "iterator (avg_nanosec_per_kmer)": "2.820549"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "941.367280", "negative lookup (avg_nanosec_per_kmer)": "581.844150", "access (avg_nanosec_per_kmer)": "622.161489", "iterator (avg_nanosec_per_kmer)": "2.827105"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "944.656877", "negative lookup (avg_nanosec_per_kmer)": "588.830562", "access (avg_nanosec_per_kmer)": "621.531063", "iterator (avg_nanosec_per_kmer)": "2.838949"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "942.756666", "negative lookup (avg_nanosec_per_kmer)": "582.045724", "access (avg_nanosec_per_kmer)": "624.022485", "iterator (avg_nanosec_per_kmer)": "2.915508"} diff --git a/benchmarks/results-28-08-26/k63/bench.log b/benchmarks/results-28-08-26/k63/bench.log new file mode 100644 index 0000000..dc29412 --- /dev/null +++ b/benchmarks/results-28-08-26/k63/bench.log @@ -0,0 +1,90 @@ +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 567.271 +negative lookup (avg_nanosec_per_kmer) 447.897 +access (avg_nanosec_per_kmer) = 292.035 +iterator (avg_nanosec_per_kmer) = 2.72352 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 567.159 +negative lookup (avg_nanosec_per_kmer) 442.781 +access (avg_nanosec_per_kmer) = 290.059 +iterator (avg_nanosec_per_kmer) = 2.71677 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 577.632 +negative lookup (avg_nanosec_per_kmer) 445.1 +access (avg_nanosec_per_kmer) = 289.093 +iterator (avg_nanosec_per_kmer) = 2.71384 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 536.517 +negative lookup (avg_nanosec_per_kmer) 479.224 +access (avg_nanosec_per_kmer) = 324.438 +iterator (avg_nanosec_per_kmer) = 2.75762 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 537.309 +negative lookup (avg_nanosec_per_kmer) 482.519 +access (avg_nanosec_per_kmer) = 327.085 +iterator (avg_nanosec_per_kmer) = 2.70532 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 543.876 +negative lookup (avg_nanosec_per_kmer) 484.284 +access (avg_nanosec_per_kmer) = 321.921 +iterator (avg_nanosec_per_kmer) = 2.70889 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 655.824 +negative lookup (avg_nanosec_per_kmer) 512.845 +access (avg_nanosec_per_kmer) = 355.445 +iterator (avg_nanosec_per_kmer) = 2.71961 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 657.564 +negative lookup (avg_nanosec_per_kmer) 514.922 +access (avg_nanosec_per_kmer) = 356.298 +iterator (avg_nanosec_per_kmer) = 2.77305 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 659.735 +negative lookup (avg_nanosec_per_kmer) 510.912 +access (avg_nanosec_per_kmer) = 355.614 +iterator (avg_nanosec_per_kmer) = 2.71592 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 531.164 +negative lookup (avg_nanosec_per_kmer) 434.533 +access (avg_nanosec_per_kmer) = 281.461 +iterator (avg_nanosec_per_kmer) = 2.73355 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 525.986 +negative lookup (avg_nanosec_per_kmer) 437.955 +access (avg_nanosec_per_kmer) = 277.253 +iterator (avg_nanosec_per_kmer) = 2.7278 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 531.091 +negative lookup (avg_nanosec_per_kmer) 434.231 +access (avg_nanosec_per_kmer) = 283.332 +iterator (avg_nanosec_per_kmer) = 2.72768 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 965.438 +negative lookup (avg_nanosec_per_kmer) 495.347 +access (avg_nanosec_per_kmer) = 402.683 +iterator (avg_nanosec_per_kmer) = 2.75541 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 969.421 +negative lookup (avg_nanosec_per_kmer) 494.921 +access (avg_nanosec_per_kmer) = 404.44 +iterator (avg_nanosec_per_kmer) = 2.74908 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 961.155 +negative lookup (avg_nanosec_per_kmer) 498.388 +access (avg_nanosec_per_kmer) = 400.903 +iterator (avg_nanosec_per_kmer) = 2.82055 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 941.367 +negative lookup (avg_nanosec_per_kmer) 581.844 +access (avg_nanosec_per_kmer) = 622.161 +iterator (avg_nanosec_per_kmer) = 2.8271 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 944.657 +negative lookup (avg_nanosec_per_kmer) 588.831 +access (avg_nanosec_per_kmer) = 621.531 +iterator (avg_nanosec_per_kmer) = 2.83895 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 942.757 +negative lookup (avg_nanosec_per_kmer) 582.046 +access (avg_nanosec_per_kmer) = 624.022 +iterator (avg_nanosec_per_kmer) = 2.91551 diff --git a/benchmarks/results-28-08-26/k63/streaming-queries.json b/benchmarks/results-28-08-26/k63/streaming-queries.json new file mode 100644 index 0000000..decc135 --- /dev/null +++ b/benchmarks/results-28-08-26/k63/streaming-queries.json @@ -0,0 +1,6 @@ +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz", "num_kmers": "97972416", "num_positive_kmers": "67275966", "num_negative_kmers": "30696450", "num_invalid_kmers": "0", "num_searches": "2280719", "num_extensions": "64995247", "num_memo_singleton": "1309608", "num_memo_light": "11103991", "elapsed_millisec": "3569"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz", "num_kmers": "461383839", "num_positive_kmers": "293470517", "num_negative_kmers": "167902332", "num_invalid_kmers": "10990", "num_searches": "6395458", "num_extensions": "287075059", "num_memo_singleton": "4872310", "num_memo_light": "56564763", "elapsed_millisec": "20048"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz", "num_kmers": "477818474", "num_positive_kmers": "406529529", "num_negative_kmers": "70615167", "num_invalid_kmers": "673778", "num_searches": "33432030", "num_extensions": "373097499", "num_memo_singleton": "3056220", "num_memo_light": "23887860", "elapsed_millisec": "52305"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz", "num_kmers": "10330949", "num_positive_kmers": "10230224", "num_negative_kmers": "99451", "num_invalid_kmers": "1274", "num_searches": "279697", "num_extensions": "9950527", "num_memo_singleton": "1727", "num_memo_light": "53202", "elapsed_millisec": "304"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz", "num_kmers": "541466405", "num_positive_kmers": "507202856", "num_negative_kmers": "34238416", "num_invalid_kmers": "25133", "num_searches": "131137192", "num_extensions": "376065664", "num_memo_singleton": "126573", "num_memo_light": "3425578", "elapsed_millisec": "150386"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz", "num_kmers": "477818474", "num_positive_kmers": "434532302", "num_negative_kmers": "42612394", "num_invalid_kmers": "673778", "num_searches": "43895917", "num_extensions": "390636385", "num_memo_singleton": "877256", "num_memo_light": "13835358", "elapsed_millisec": "65020"} diff --git a/benchmarks/results-28-08-26/k63/streaming-queries.log b/benchmarks/results-28-08-26/k63/streaming-queries.log new file mode 100644 index 0000000..89312c6 --- /dev/null +++ b/benchmarks/results-28-08-26/k63/streaming-queries.log @@ -0,0 +1,78 @@ +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz +2026-08-28 20:41:29: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz'... +2026-08-28 20:41:32: DONE +==== query report: +num_kmers = 97972416 +num_positive_kmers = 67275966 (68.6683%) +num_negative_kmers = 30696450 (31.3317%) +num_invalid_kmers = 0 (0%) +num_searches = 2280719/67275966 (3.39009%) +num_extensions = 64995247/67275966 (96.6099%) +num_memo_singleton = 1309608 +num_memo_light = 11103991 +elapsed = 3.569 sec / 0.0594833 min / 36.4286 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz +2026-08-28 20:41:33: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz'... +2026-08-28 20:41:53: DONE +==== query report: +num_kmers = 461383839 +num_positive_kmers = 293470517 (63.6066%) +num_negative_kmers = 167902332 (36.391%) +num_invalid_kmers = 10990 (0.00238196%) +num_searches = 6395458/293470517 (2.17925%) +num_extensions = 287075059/293470517 (97.8207%) +num_memo_singleton = 4872310 +num_memo_light = 56564763 +elapsed = 20.048 sec / 0.334133 min / 43.4519 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz +2026-08-28 20:41:53: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz'... +2026-08-28 20:42:46: DONE +==== query report: +num_kmers = 477818474 +num_positive_kmers = 406529529 (85.0803%) +num_negative_kmers = 70615167 (14.7787%) +num_invalid_kmers = 673778 (0.141011%) +num_searches = 33432030/406529529 (8.22376%) +num_extensions = 373097499/406529529 (91.7762%) +num_memo_singleton = 3056220 +num_memo_light = 23887860 +elapsed = 52.305 sec / 0.87175 min / 109.466 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz +2026-08-28 20:42:46: performing queries from file '/mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz'... +2026-08-28 20:42:46: DONE +==== query report: +num_kmers = 10330949 +num_positive_kmers = 10230224 (99.025%) +num_negative_kmers = 99451 (0.962651%) +num_invalid_kmers = 1274 (0.0123319%) +num_searches = 279697/10230224 (2.73403%) +num_extensions = 9950527/10230224 (97.266%) +num_memo_singleton = 1727 +num_memo_light = 53202 +elapsed = 0.304 sec / 0.00506667 min / 29.4261 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz +2026-08-28 20:42:47: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz'... +2026-08-28 20:45:17: DONE +==== query report: +num_kmers = 541466405 +num_positive_kmers = 507202856 (93.6721%) +num_negative_kmers = 34238416 (6.32328%) +num_invalid_kmers = 25133 (0.00464165%) +num_searches = 131137192/507202856 (25.855%) +num_extensions = 376065664/507202856 (74.145%) +num_memo_singleton = 126573 +num_memo_light = 3425578 +elapsed = 150.386 sec / 2.50643 min / 277.738 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz +2026-08-28 20:45:20: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz'... +2026-08-28 20:46:25: DONE +==== query report: +num_kmers = 477818474 +num_positive_kmers = 434532302 (90.9409%) +num_negative_kmers = 42612394 (8.91811%) +num_invalid_kmers = 673778 (0.141011%) +num_searches = 43895917/434532302 (10.1019%) +num_extensions = 390636385/434532302 (89.8981%) +num_memo_singleton = 877256 +num_memo_light = 13835358 +elapsed = 65.02 sec / 1.08367 min / 136.077 ns/kmer diff --git a/include/builder/dictionary_builder.hpp b/include/builder/dictionary_builder.hpp index dea89a5..5f63d54 100644 --- a/include/builder/dictionary_builder.hpp +++ b/include/builder/dictionary_builder.hpp @@ -20,13 +20,11 @@ struct dictionary_builder // d.m_m = build_config.m; d.m_spss.k = build_config.k; d.m_spss.m = build_config.m; - d.m_canonical = build_config.canonical; d.m_hasher.seed(build_config.seed); build_stats.add("input_filename", filename.c_str()); build_stats.add("k", d.m_k); build_stats.add("m", d.m_m); - build_stats.add("canonical", d.m_canonical ? "true" : "false"); build_stats.add("seed", build_config.seed); build_stats.add("num_threads", build_config.num_threads); @@ -50,7 +48,6 @@ struct dictionary_builder // std::cout << "num_minimizers = " << minimizers.num_minimizers() << std::endl; std::cout << "num_minimizer_positions = " << minimizers.num_minimizer_positions() << std::endl; - std::cout << "num_super_kmers = " << minimizers.num_super_kmers() << std::endl; } do_step("step 4 (build mphf)", [&]() { build_mphf(d); }); @@ -143,21 +140,22 @@ struct dictionary_builder // } } - const uint64_t num_super_kmers = minimizers.num_super_kmers(); + /* one tuple per minimizer position (the scheme is forward), so this is + the number of records in the tuples file */ + const uint64_t num_tuples = minimizers.num_minimizer_positions(); const uint64_t buffer_size = num_files_to_merge == 1 - ? num_super_kmers + ? num_tuples : (RAM_available_in_bytes / (3 * sizeof(minimizer_tuple))); - const uint64_t num_blocks = (num_super_kmers + buffer_size - 1) / buffer_size; - assert(num_super_kmers > (num_blocks - 1) * buffer_size); + const uint64_t num_blocks = (num_tuples + buffer_size - 1) / buffer_size; + assert(num_tuples > (num_blocks - 1) * buffer_size); std::vector threads; threads.reserve(num_threads); std::vector buffer; for (uint64_t i = 0; i != num_blocks; ++i) { - const uint64_t n = (i == num_blocks - 1) - ? num_super_kmers - (num_blocks - 1) * buffer_size - : buffer_size; + const uint64_t n = + (i == num_blocks - 1) ? num_tuples - (num_blocks - 1) * buffer_size : buffer_size; buffer.resize(n); input.read(reinterpret_cast(buffer.data()), buffer.size() * sizeof(minimizer_tuple)); diff --git a/include/builder/util.hpp b/include/builder/util.hpp index b390436..add5e00 100644 --- a/include/builder/util.hpp +++ b/include/builder/util.hpp @@ -58,23 +58,24 @@ inline std::ostream& operator<<(std::ostream& os, minimizer_tuple const& mt) { return os; } +/* + The bucket of a minimizer: its tuples, sorted by pos_in_seq. The minimizer + is forward (see `util::compute_minimizer`), so a locus is never abandoned + and later re-selected: each tuple carries a distinct pos_in_seq, and the + number of tuples (= super-kmers) equals the number of minimizer positions, + which is what size() returns. `minimizers_tuples::merge` checks this. +*/ struct bucket_type { bucket_type(minimizer_tuple const* begin, minimizer_tuple const* end) : m_begin(begin) - , m_end(end) - , m_num_super_kmers(std::distance(begin, end)) - , m_num_minimizer_positions(0) // + , m_end(end) // { - uint64_t prev_pos_in_seq = constants::invalid_uint64; - while (begin != end) { - uint64_t pos_in_seq = (*begin).pos_in_seq; - if (pos_in_seq != prev_pos_in_seq) { - ++m_num_minimizer_positions; - prev_pos_in_seq = pos_in_seq; - } - ++begin; - } - assert(m_num_minimizer_positions <= m_num_super_kmers); + assert([&] { /* one tuple per minimizer position: the scheme is forward */ + for (auto it = begin; it + 1 < end; ++it) { + if (it->pos_in_seq == (it + 1)->pos_in_seq) return false; + } + return true; + }()); } struct iterator { @@ -92,19 +93,7 @@ struct bucket_type { iterator begin() const { return iterator(m_begin); } iterator end() const { return iterator(m_end); } - /* - When a canonical index is built (option `--canonical`), - a minimizer offset can correspond to more than one super-kmer. - A super-kmer is uniquely identified by the couple - (minimizer offset, position of minimizer in the first kmer of the super-kmer). - These two components, together, give the - starting position of a super-kmer in the sequence. - - So the method size() returns the number of minimizer - positions which is <= the number of superkmers. - */ - uint64_t num_super_kmers() const { return m_num_super_kmers; } - uint64_t size() const { return m_num_minimizer_positions; } + uint64_t size() const { return std::distance(m_begin, m_end); } minimizer_tuple const* begin_ptr() const { return m_begin; } minimizer_tuple const* end_ptr() const { return m_end; } @@ -112,8 +101,6 @@ struct bucket_type { private: minimizer_tuple const* m_begin; minimizer_tuple const* m_end; - uint64_t m_num_super_kmers; - uint64_t m_num_minimizer_positions; }; /* @@ -159,7 +146,6 @@ struct minimizers_tuples { minimizers_tuples(build_configuration const& build_config) : m_num_minimizers(0) , m_num_minimizer_positions(0) - , m_num_super_kmers(0) , m_run_identifier(pthash::clock_type::now().time_since_epoch().count()) , m_build_config(build_config) // { @@ -217,17 +203,9 @@ struct minimizers_tuples { assert(m_num_minimizers == 0); assert(m_num_minimizer_positions == 0); - assert(m_num_super_kmers == 0); mm::file_source input(get_minimizers_filename(), mm::advice::sequential); - for (minimizers_tuples_iterator it(input.data(), input.data() + input.size()); - it.has_next(); it.next()) // - { - auto bucket = it.bucket(); - m_num_minimizers += 1; - m_num_minimizer_positions += bucket.size(); - m_num_super_kmers += bucket.num_super_kmers(); - } + count_and_check_forward(input.data(), input.data() + input.size()); input.close(); return; } @@ -245,23 +223,15 @@ struct minimizers_tuples { m_num_minimizers = 0; m_num_minimizer_positions = 0; - m_num_super_kmers = 0; uint64_t prev_minimizer = constants::invalid_uint64; uint64_t prev_pos_in_seq = constants::invalid_uint64; while (fm_iterator.has_next()) { minimizer_tuple mt = *fm_iterator; - if (mt.minimizer != prev_minimizer) { - prev_minimizer = mt.minimizer; - ++m_num_minimizers; - ++m_num_minimizer_positions; - } else { - if (mt.pos_in_seq != prev_pos_in_seq) ++m_num_minimizer_positions; - } + count_and_check_forward_one(mt, prev_minimizer, prev_pos_in_seq); out.write(reinterpret_cast(&mt), sizeof(minimizer_tuple)); - prev_pos_in_seq = mt.pos_in_seq; - ++m_num_super_kmers; - if (m_build_config.verbose and m_num_super_kmers % 100'000'000 == 0) { - std::cout << "processed " << m_num_super_kmers << " minimizer tuples" << std::endl; + if (m_build_config.verbose and m_num_minimizer_positions % 100'000'000 == 0) { + std::cout << "processed " << m_num_minimizer_positions << " minimizer tuples" + << std::endl; } fm_iterator.next(); } @@ -278,8 +248,11 @@ struct minimizers_tuples { uint64_t num_files_to_merge() const { return m_num_files_to_merge; } uint64_t num_minimizers() const { return m_num_minimizers; } + + /* One tuple per minimizer position, the scheme being forward (checked by + `merge`), so this is also the number of super-kmers and the number of + records in the merged tuples file. */ uint64_t num_minimizer_positions() const { return m_num_minimizer_positions; } - uint64_t num_super_kmers() const { return m_num_super_kmers; } void remove_tmp_file() { std::remove(get_minimizers_filename().c_str()); } @@ -287,10 +260,39 @@ struct minimizers_tuples { std::atomic m_num_files_to_merge; uint64_t m_num_minimizers; uint64_t m_num_minimizer_positions; - uint64_t m_num_super_kmers; uint64_t m_run_identifier; build_configuration m_build_config; + /* + Count minimizers and minimizer positions over tuples sorted by + (minimizer, pos_in_seq), checking the forwardness requirement: the + scheme never re-selects an abandoned locus, so no (minimizer, + pos_in_seq) pair may appear twice. Everything downstream sizes the + index by the tuple count, so a violation must stop the build. + */ + void count_and_check_forward_one(minimizer_tuple const& mt, uint64_t& prev_minimizer, + uint64_t& prev_pos_in_seq) // + { + if (mt.minimizer != prev_minimizer) { + prev_minimizer = mt.minimizer; + ++m_num_minimizers; + } else if (mt.pos_in_seq == prev_pos_in_seq) { + throw std::runtime_error( + "the minimizer scheme is not forward: " + "a (minimizer, position) pair occurs in more than one super-kmer"); + } + prev_pos_in_seq = mt.pos_in_seq; + ++m_num_minimizer_positions; + } + + void count_and_check_forward(minimizer_tuple const* begin, minimizer_tuple const* end) { + uint64_t prev_minimizer = constants::invalid_uint64; + uint64_t prev_pos_in_seq = constants::invalid_uint64; + for (; begin != end; ++begin) { + count_and_check_forward_one(*begin, prev_minimizer, prev_pos_in_seq); + } + } + std::string get_tmp_output_filename(uint64_t id) const { std::stringstream filename; filename << m_build_config.tmp_dirname << "/sshash.tmp.run_" << m_run_identifier diff --git a/include/constants.hpp b/include/constants.hpp index b252f6e..215901a 100644 --- a/include/constants.hpp +++ b/include/constants.hpp @@ -20,9 +20,9 @@ constexpr int forward_orientation = 1; constexpr int backward_orientation = -1; namespace current_version_number { -constexpr uint8_t x = 5; -constexpr uint8_t y = 1; -constexpr uint8_t z = 1; +constexpr uint8_t x = 6; +constexpr uint8_t y = 0; +constexpr uint8_t z = 0; } // namespace current_version_number } // namespace sshash::constants diff --git a/include/dictionary.hpp b/include/dictionary.hpp index a30b8c4..b5fb618 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -11,6 +11,7 @@ template struct dictionary // { using kmer_type = Kmer; + using spss_type = spectrum_preserving_string_set; template friend struct dictionary_builder; @@ -22,8 +23,7 @@ struct dictionary // , m_num_kmers(0) , m_num_strings(0) , m_k(0) - , m_m(0) - , m_canonical(false) {} + , m_m(0) {} /* Build from input file. */ void build(std::string const& input_filename, build_configuration const& build_config); @@ -33,11 +33,16 @@ struct dictionary // uint64_t num_strings() const { return m_num_strings; } uint64_t k() const { return m_k; } uint64_t m() const { return m_m; } - bool canonical() const { return m_canonical; } bool weighted() const { return !m_weights.empty(); } hasher_type const& hasher() const { return m_hasher; } - /* Lookup queries. */ + /* + Lookup queries. A kmer and its reverse complement share a bucket, so a + lookup always costs a single probe and reports, via + `lookup_result::kmer_orientation`, which of the two was found. Pass + `check_reverse_complement = false` to restrict the answer to kmers + occurring in forward orientation. + */ lookup_result lookup(char const* string_kmer, bool check_reverse_complement = true) const; lookup_result lookup(Kmer uint_kmer, bool check_reverse_complement = true) const; @@ -75,7 +80,7 @@ struct dictionary // bool is_member(char const* string_kmer, bool check_reverse_complement = true) const; bool is_member(Kmer uint_kmer, bool check_reverse_complement = true) const; - template + template friend struct streaming_query; streaming_query_report // @@ -144,7 +149,6 @@ struct dictionary // visitor.visit(t.m_num_strings); visitor.visit(t.m_k); visitor.visit(t.m_m); - visitor.visit(t.m_canonical); visitor.visit(t.m_hasher); visitor.visit(t.m_spss); visitor.visit(t.m_ssi); @@ -156,20 +160,15 @@ struct dictionary // uint64_t m_num_strings; uint16_t m_k; uint16_t m_m; - bool m_canonical; hasher_type m_hasher; - spectrum_preserving_string_set m_spss; + spss_type m_spss; sparse_and_skew_index m_ssi; weights m_weights; - lookup_result lookup_regular(Kmer uint_kmer) const; - lookup_result lookup_regular(Kmer uint_kmer, minimizer_info mini_info) const; - - lookup_result lookup_canonical(Kmer uint_kmer) const; - lookup_result lookup_canonical(Kmer uint_kmer, Kmer uint_kmer_rc, - minimizer_info mini_info) const; + lookup_result lookup(Kmer uint_kmer, Kmer uint_kmer_rc, minimizer_info mini_info, + typename spss_type::bucket_cache* cache = nullptr) const; void forward_neighbours(Kmer suffix, neighbourhood& res, bool check_reverse_complement) const; diff --git a/include/kmer.hpp b/include/kmer.hpp index 6fe1520..3dae600 100644 --- a/include/kmer.hpp +++ b/include/kmer.hpp @@ -105,7 +105,15 @@ struct alpha_kmer_t : uint_kmer_t { static char uint64_to_char(uint64_t x) { return alphabet[x]; } // Revcompl only makes sense for DNA, fallback to noop otherwise + static constexpr bool has_reverse_complement = false; [[maybe_unused]] virtual void reverse_complement_inplace(uint64_t) {} + + /* Reverse complement of an m-mer packed in the low m*bits_per_char bits of a + word, for m <= max_m. Same fallback as above: the identity. */ + [[maybe_unused]] static uint64_t reverse_complement_mmer(uint64_t mmer, uint64_t) { + return mmer; + } + [[maybe_unused]] static void compute_reverse_complement(char const* input, char* output, uint64_t size) { for (uint64_t i = 0; i != size; ++i) output[i] = input[i]; @@ -156,6 +164,8 @@ struct dna_uint_kmer_t : alpha_kmer_t { return res; } + static constexpr bool has_reverse_complement = true; + [[maybe_unused]] void reverse_complement_inplace(uint64_t k) override { assert(k <= max_k); dna_uint_kmer_t rev(0); @@ -164,6 +174,20 @@ struct dna_uint_kmer_t : alpha_kmer_t { *this = rev; } + /* + Reverse complement of an m-mer packed in the low m*bits_per_char bits of a + single word. This is `reverse_complement_inplace` specialized to a value + that is known to fit in 64 bits, hence a single crc64 regardless of how + wide the kmer type is. Used to canonicalize the m-mers of a kmer, which is + done once per character of the input during construction and once per + character of the query at lookup time. + */ + [[maybe_unused]] static uint64_t reverse_complement_mmer(uint64_t mmer, uint64_t m) { + assert(m <= max_m); + assert(m * bits_per_char < 64); + return crc64(mmer) >> (64 - m * bits_per_char); + } + #ifdef SSHASH_USE_TRADITIONAL_NUCLEOTIDE_ENCODING /* char decimal binary diff --git a/include/kmer_iterator.hpp b/include/kmer_iterator.hpp index 8672ab1..356db0e 100644 --- a/include/kmer_iterator.hpp +++ b/include/kmer_iterator.hpp @@ -72,7 +72,10 @@ struct kmer_iterator // inline void fill_buff_reverse() { static_assert(Kmer::uint_kmer_bits % 64 == 0); - for (int i = Kmer::uint_kmer_bits; i > 0; i -= 64) { + /* fill the buffer with the bits at [m_pos - uint_kmer_bits, m_pos), + in order: `append64` appends at the low end, so words must be + appended from the highest address down to the lowest */ + for (int i = 64; i <= Kmer::uint_kmer_bits; i += 64) { m_buff.append64(m_bv->get_word64(std::max(m_pos, Kmer::uint_kmer_bits) - i)); } m_avail = std::min(m_pos, Kmer::uint_kmer_bits); diff --git a/include/minimizer_iterator.hpp b/include/minimizer_iterator.hpp index b53adc8..56f26f0 100644 --- a/include/minimizer_iterator.hpp +++ b/include/minimizer_iterator.hpp @@ -1,11 +1,26 @@ #pragma once #include "kmer.hpp" +#include "util.hpp" namespace sshash { /* - "Re-scan" method. + "Re-scan" method: computes the minimizer of each kmer of a sequence, + sliding one character at a time. + + See `util::compute_minimizer` for the definition of the minimizer; this + iterator computes exactly the same thing incrementally, which the assertion + at the end of `next` checks. + + The only extra state compared to a plain forward minimizer is `m_num_mins`, + the number of loci of the current window attaining the minimum hash: a tie + cannot be broken by position alone without breaking mirror-equivariance, so + when there is one we have to look at the kmer's own orientation. + + The kmer is reverse-complemented once per call and that single value serves + every locus, via rc(x_i) = rc(x)_{k-m-i}, so no m-mer is ever + reverse-complemented on its own. */ template struct minimizer_iterator { @@ -30,141 +45,125 @@ struct minimizer_iterator { void reset() { m_min_pos_in_kmer = 0; m_min_position = m_position - 1; + m_num_mins = 0; } minimizer_info next(Kmer kmer) { + Kmer kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(m_k); + if (m_min_pos_in_kmer == 0) { /* min leaves the window: re-scan to compute the new min */ m_position = m_min_position + 1; - rescan(kmer); + rescan(kmer, kmer_rc); } else { m_position += 1; - Kmer mmer = kmer; - mmer.drop_chars(m_k - m_m); - uint64_t hash = m_hasher.hash(uint64_t(mmer)); + Kmer window = kmer; + window.drop_chars(m_k - m_m); + uint64_t value = util::canonical_mmer_at(window, kmer_rc, m_k, m_m, m_k - m_m); + uint64_t hash = m_hasher.hash(value); if (hash < m_min_hash) { m_min_hash = hash; - m_min_value = uint64_t(mmer); + m_min_value = value; m_min_position = m_position; m_min_pos_in_kmer = m_k - m_m; + m_num_mins = 1; } else { + /* Only the leftmost minimum can leave the window without a + re-scan, so the count stays valid across the slide. */ assert(m_min_pos_in_kmer > 0); m_min_pos_in_kmer -= 1; + if (hash == m_min_hash) m_num_mins += 1; } } - assert(minimizer_info(m_min_value, m_min_pos_in_kmer) == - util::compute_minimizer(kmer, m_k, m_m, m_hasher)); + minimizer_info mini_info(m_min_value, m_min_position, m_min_pos_in_kmer); + if (m_num_mins > 1) break_tie(kmer, kmer_rc, mini_info); - return {m_min_value, m_min_position, m_min_pos_in_kmer}; + assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == + util::compute_minimizer(kmer, kmer_rc, m_k, m_m, m_hasher)); + + return mini_info; } private: uint64_t m_k, m_m; uint64_t m_position, m_min_pos_in_kmer; uint64_t m_min_value, m_min_position, m_min_hash; + uint64_t m_num_mins; hasher_type m_hasher; - void rescan(Kmer kmer) { - m_min_hash = constants::invalid_uint64; - m_min_value = constants::invalid_uint64; - m_min_pos_in_kmer = 0; - uint64_t begin = m_position; - for (uint64_t i = 0; i != m_k - m_m + 1; ++i, ++m_position) { - Kmer mmer = kmer; + void rescan(Kmer kmer, Kmer const& kmer_rc) { + const uint64_t begin = m_position; + + /* first locus, peeled off the loop: see `util::compute_minimizer` */ + { + m_min_value = util::canonical_mmer_at(kmer, kmer_rc, m_k, m_m, 0); kmer.drop_char(); - mmer.take_chars(m_m); - uint64_t hash = m_hasher.hash(uint64_t(mmer)); + m_min_hash = m_hasher.hash(m_min_value); + m_min_pos_in_kmer = 0; + m_num_mins = 1; + ++m_position; + } + + for (uint64_t i = 1; i != m_k - m_m + 1; ++i, ++m_position) { + uint64_t value = util::canonical_mmer_at(kmer, kmer_rc, m_k, m_m, i); + kmer.drop_char(); + uint64_t hash = m_hasher.hash(value); if (hash < m_min_hash) { // leftmost m_min_hash = hash; - m_min_value = uint64_t(mmer); + m_min_value = value; m_min_pos_in_kmer = i; + m_num_mins = 1; + } else if (hash == m_min_hash) { + m_num_mins += 1; } } + m_position -= 1; m_min_position = begin + m_min_pos_in_kmer; } -}; - -/* - "Re-scan" method. -*/ -template -struct minimizer_iterator_rc { - minimizer_iterator_rc() {} - - minimizer_iterator_rc(uint64_t k, uint64_t m, hasher_type const& hasher, uint64_t position = 0) - : m_k(k) - , m_m(m) - , m_min_value(constants::invalid_uint64) - , m_min_hash(constants::invalid_uint64) - , m_hasher(hasher) // - { - assert(k > 0 and m <= k); - set_position(position); - } - - void set_position(uint64_t position) { - m_position = position; - reset(); - } - - void reset() { - m_min_pos_in_kmer = m_k - m_m; - m_min_position = m_position - 1; - } - minimizer_info next(Kmer kmer) { - if (m_min_pos_in_kmer == m_k - m_m) { - /* min leaves the window: re-scan to compute the new min */ - m_position = m_min_position + 1; - rescan(kmer); - } else { - m_position += 1; - Kmer mmer = kmer; - mmer.take_chars(m_m); - uint64_t hash = m_hasher.hash(uint64_t(mmer)); - if (hash <= m_min_hash) { - m_min_hash = hash; - m_min_value = uint64_t(mmer); - m_min_position = m_position; - m_min_pos_in_kmer = 0; - } else { - m_min_pos_in_kmer += 1; - assert(m_min_pos_in_kmer <= m_k - m_m); + /* + Two or more loci of the window attain the minimum hash (all carrying the + same class, so only the sampled position is at stake). Apply the + centre-closest tie-break of `util::compute_minimizer`: take the tied + locus closest to the window centre, and between two equally close ones + -- mirror images i and k-m-i -- the smaller index if kmer <= rc(kmer), + the larger otherwise. Forward and mirror-equivariant. It runs on ~1e-4 + of the windows, so the rescan below costs nothing overall. + + Tied loci cannot precede the leftmost minimum, so the scan starts there. + */ + void break_tie(Kmer kmer, Kmer const& kmer_rc, minimizer_info& mini_info) const { + assert(m_num_mins > 1); + + const uint64_t window_begin = m_min_position - m_min_pos_in_kmer; + const uint64_t two_c = m_k - m_m; + uint64_t best_dist = constants::invalid_uint64; + uint64_t lo = m_min_pos_in_kmer; + uint64_t hi = m_min_pos_in_kmer; + + Kmer window = kmer; + window.drop_chars(m_min_pos_in_kmer); + for (uint64_t i = m_min_pos_in_kmer; i != m_k - m_m + 1; ++i) { + const uint64_t v = util::canonical_mmer_at(window, kmer_rc, m_k, m_m, i); + if (m_hasher.hash(v) == m_min_hash) { + assert(v == m_min_value); + const uint64_t dist = 2 * i > two_c ? 2 * i - two_c : two_c - 2 * i; + if (dist < best_dist) { + best_dist = dist; + lo = i; + hi = i; + } else if (dist == best_dist) { + hi = i; + } } + window.drop_char(); } - assert(minimizer_info(m_min_value, m_min_pos_in_kmer) == - util::compute_minimizer(kmer, m_k, m_m, m_hasher)); - - return {m_min_value, m_min_position, m_min_pos_in_kmer}; - } - -private: - uint64_t m_k, m_m; - uint64_t m_position, m_min_pos_in_kmer; - uint64_t m_min_value, m_min_position, m_min_hash; - hasher_type m_hasher; - - void rescan(Kmer kmer) { - m_min_hash = constants::invalid_uint64; - m_min_value = constants::invalid_uint64; - m_min_pos_in_kmer = 0; - uint64_t begin = m_position; - for (int64_t i = m_k - m_m; i >= 0; --i, ++m_position) { - Kmer mmer = kmer; - mmer.drop_chars(i); - mmer.take_chars(m_m); - uint64_t hash = m_hasher.hash(uint64_t(mmer)); - if (hash <= m_min_hash) { // rightmost - m_min_hash = hash; - m_min_value = uint64_t(mmer); - m_min_pos_in_kmer = i; - } - } - m_position -= 1; - m_min_position = begin + (m_k - m_min_pos_in_kmer - m_m); + const uint64_t chosen = (lo == hi or !(kmer_rc < kmer)) ? lo : hi; + mini_info = minimizer_info(m_min_value, window_begin + chosen, chosen); } }; diff --git a/include/sparse_and_skew_index.hpp b/include/sparse_and_skew_index.hpp index 4509558..96f73fc 100644 --- a/include/sparse_and_skew_index.hpp +++ b/include/sparse_and_skew_index.hpp @@ -109,7 +109,9 @@ struct sparse_and_skew_index // bits::compact_vector::iterator m_it; }; - bucket_iterator lookup(const Kmer uint_kmer, const minimizer_info mini_info) const // + bucket_iterator lookup(const Kmer uint_kmer, // + const Kmer uint_kmer_rc, // + const minimizer_info mini_info) const // { uint64_t code = codewords.lookup(mini_info.minimizer); @@ -132,7 +134,8 @@ struct sparse_and_skew_index // } assert(status == bucket_t::HEAVYLOAD); // minimizer is part of the skew index - uint64_t offset = ski.lookup(uint_kmer, code); + const Kmer uint_kmer_canon = std::min(uint_kmer, uint_kmer_rc); + uint64_t offset = ski.lookup(uint_kmer_canon, code); return {this, offset, 1, bucket_t::HEAVYLOAD}; } diff --git a/include/spectrum_preserving_string_set.hpp b/include/spectrum_preserving_string_set.hpp index feca86b..b5fd5eb 100644 --- a/include/spectrum_preserving_string_set.hpp +++ b/include/spectrum_preserving_string_set.hpp @@ -26,30 +26,58 @@ struct spectrum_preserving_string_set // return util::read_kmer_at(strings, k - 1, Kmer::bits_per_char * (string_end - k + 1)); } + /* + The decoded locate set of one minimizer, filled by `lookup` (which + decodes straight into it, at no extra cost) and consumed by + `lookup_from_positions`: it lets streaming queries verify kmers + sharing the minimizer of their last seed without a codeword lookup + nor offset decoding. Only buckets of known size up to + `max_cached_size` are cached: HEAVYLOAD buckets have unknown size, + and larger buckets are too rare to be worth a bigger cache. + Buckets whose minimizer is not found are never cached (size = 0). + */ + struct bucket_cache { + static constexpr uint64_t max_cached_size = 8; + uint64_t size = 0; // 0 = nothing cached + std::array positions; + }; + template - lookup_result lookup_regular(Iterator it, // - const Kmer kmer, // - const minimizer_info mini_info) const // + lookup_result lookup(Iterator it, // + const Kmer kmer, const Kmer kmer_rc, // + const minimizer_info mini_info, // + bucket_cache* cache = nullptr) const // { const uint64_t size = it.size(); assert(size > 0); static thread_local // std::array - v; + tl; + + auto* v = tl.data(); + if (cache != nullptr) { + cache->size = 0; + if (it.bucket_type() != bucket_t::HEAVYLOAD and size <= bucket_cache::max_cached_size) { + v = cache->positions.data(); + } + } for (uint64_t i = 0; i != size; ++i, ++it) { uint64_t minimizer_offset = *it; v[i] = strings_offsets.decode(minimizer_offset); } - /* check minimizer first */ + /* Check minimizer first. The minimizer's value is the canonical m-mer at + the anchored locus, so the m-mer stored at that offset is either the + minimizer itself or its reverse complement. */ if (uint64_t read_mmer = uint64_t( util::read_kmer_at(strings, m, Kmer::bits_per_char * v[0].absolute_offset)); - read_mmer != mini_info.minimizer) // + read_mmer != mini_info.minimizer and + util::canonical_mmer(read_mmer, m) != mini_info.minimizer) // { /* - The function `lookup_regular` determines if the minimizer is found at the + This function determines if the minimizer is found at the offset `Kmer::bits_per_char * p.absolute_offset`, not whether the minimizer does not appear at all. In fact, it can happen that the minimizer appear but not at the specified offset, so it would be wrong to set `res.minimizer_found` @@ -64,50 +92,26 @@ struct spectrum_preserving_string_set // return lookup_result(it.bucket_type() != bucket_t::HEAVYLOAD ? false : true); } - lookup_result res; - for (uint64_t i = 0; i != size; ++i) { - if (_lookup_regular(res, v[i], kmer, mini_info)) return res; - } + if (v != tl.data()) cache->size = size; - return lookup_result(); + return lookup_from_positions(v, size, kmer, kmer_rc, mini_info); } - template - lookup_result lookup_canonical(Iterator it, // - const Kmer kmer, const Kmer kmer_rc, // - const minimizer_info mini_info) const // + /* + The verification half of `lookup`: try the candidate loci of each + position of the minimizer's locate set. The positions must be the + decoded locate set of `mini_info.minimizer`, whose presence at + `positions[0]` has already been checked. + */ + lookup_result lookup_from_positions(typename Offsets::decoded_offset const* positions, // + const uint64_t size, // + const Kmer kmer, const Kmer kmer_rc, // + const minimizer_info mini_info) const // { - const uint64_t size = it.size(); - assert(size > 0); - - static thread_local // - std::array - v; - - for (uint64_t i = 0; i != size; ++i, ++it) { - uint64_t minimizer_offset = *it; - v[i] = strings_offsets.decode(minimizer_offset); - } - - /* check minimizer first */ - if (uint64_t read_mmer = uint64_t( - util::read_kmer_at(strings, m, Kmer::bits_per_char * v[0].absolute_offset)); - read_mmer != mini_info.minimizer) // - { - Kmer tmp = mini_info.minimizer; - tmp.reverse_complement_inplace(m); - uint64_t minimizer_rc = uint64_t(tmp); - if (read_mmer != minimizer_rc) { - /* Same note as for the function `lookup_regular`. */ - return lookup_result(it.bucket_type() != bucket_t::HEAVYLOAD ? false : true); - } - } - lookup_result res; for (uint64_t i = 0; i != size; ++i) { - if (_lookup_canonical(res, v[i], kmer, kmer_rc, mini_info)) return res; + if (_lookup(res, positions[i], kmer, kmer_rc, mini_info)) return res; } - return lookup_result(); } @@ -210,47 +214,34 @@ struct spectrum_preserving_string_set // visitor.visit(t.strings); } - bool _lookup_regular(lookup_result& res, // - typename Offsets::decoded_offset p, // - const Kmer kmer, // - const minimizer_info mini_info) const // + /* + The minimizer is anchored at locus `pos_in_kmer` of the kmer and, by + mirror-equivariance, at locus `k - m - pos_in_kmer` of its reverse + complement. So a minimizer occurrence at offset j is the anchor of a kmer + starting either at j - pos_in_kmer or at j - (k - m - pos_in_kmer): + always exactly two candidates, with no case analysis. They lie within + 2k-m characters of each other, hence usually on the same cache line. + */ + bool _lookup(lookup_result& res, // + typename Offsets::decoded_offset p, // + const Kmer kmer, // + const Kmer kmer_rc, // + const minimizer_info mini_info) const // { - if (p.absolute_offset < mini_info.pos_in_kmer) return false; - - res.kmer_offset = p.absolute_offset - mini_info.pos_in_kmer; - - if (kmer != util::read_kmer_at(strings, k, Kmer::bits_per_char * res.kmer_offset)) { - return false; - } - - if (res.kmer_offset >= res.string_begin and res.kmer_offset < res.string_end - k + 1) { - res.kmer_id = res.kmer_offset - res.string_id * (k - 1); // absolute kmer id - res.kmer_id_in_string = res.kmer_offset - res.string_begin; // relative kmer id + if constexpr (Kmer::has_reverse_complement) { + if (_lookup_at(res, p, kmer, kmer_rc, mini_info.pos_in_kmer)) return true; + return _lookup_at(res, p, kmer, kmer_rc, k - m - mini_info.pos_in_kmer); } else { - strings_offsets.offset_to_id(res, p, k); + /* no reverse complement: the two candidates would coincide */ + return _lookup_at(res, p, kmer, kmer_rc, mini_info.pos_in_kmer); } - - if (res.kmer_offset < res.string_end - k + 1) return true; - return false; - } - - bool _lookup_canonical(lookup_result& res, // - typename Offsets::decoded_offset p, // - const Kmer kmer, // - const Kmer kmer_rc, // - const minimizer_info mini_info) const // - { - uint64_t pos_in_kmer = mini_info.pos_in_kmer; - if (__lookup_canonical(res, p, kmer, kmer_rc, pos_in_kmer)) return true; - pos_in_kmer = k - m - mini_info.pos_in_kmer; - return __lookup_canonical(res, p, kmer, kmer_rc, pos_in_kmer); } - bool __lookup_canonical(lookup_result& res, // - typename Offsets::decoded_offset p, // - const Kmer kmer, // - const Kmer kmer_rc, // - const uint64_t pos_in_kmer) const // + bool _lookup_at(lookup_result& res, // + typename Offsets::decoded_offset p, // + const Kmer kmer, // + const Kmer kmer_rc, // + const uint64_t pos_in_kmer) const // { if (p.absolute_offset < pos_in_kmer) return false; diff --git a/include/streaming_query.hpp b/include/streaming_query.hpp index e5b8bd7..d9c4a8d 100644 --- a/include/streaming_query.hpp +++ b/include/streaming_query.hpp @@ -6,7 +6,7 @@ namespace sshash { -template +template struct streaming_query // { using kmer_t = typename Dict::kmer_type; @@ -22,35 +22,31 @@ struct streaming_query // , m_m(dict->m_m) , m_minimizer_it(dict->m_k, dict->m_m, dict->m_hasher) - , m_minimizer_it_rc(dict->m_k, dict->m_m, dict->m_hasher) , m_curr_mini_info() - , m_prev_mini_info() - , m_curr_mini_info_rc() - , m_prev_mini_info_rc() , m_it(dict->m_spss.strings, m_k) , m_remaining_string_bases(0) + , m_last_seed_mini_info() + , m_last_seed_positive(false) + , m_num_searches(0) , m_num_extensions(0) , m_num_invalid(0) , m_num_negative(0) + , m_num_skipped_singleton_lookups(0) + , m_num_bucket_cache_hits(0) - { - if (canonical != m_dict->m_canonical) { - std::stringstream ss; - ss << "dict.canonical() = " << (m_dict->canonical() ? "true" : "false") - << " but required " << (canonical ? "true" : "false"); - throw std::runtime_error(ss.str()); - } - } + {} void reset() { m_start = true; m_remaining_string_bases = 0; m_res = lookup_result(); m_minimizer_it.reset(); - m_minimizer_it_rc.reset(); + m_last_seed_mini_info = minimizer_info(); + m_last_seed_positive = false; + m_bucket_cache.size = 0; } lookup_result lookup(char const* kmer) // @@ -80,7 +76,6 @@ struct streaming_query // } m_curr_mini_info = m_minimizer_it.next(m_kmer); - m_curr_mini_info_rc = m_minimizer_it_rc.next(m_kmer_rc); /* 3. compute result */ if (m_remaining_string_bases == 0) { @@ -100,8 +95,6 @@ struct streaming_query // } /* 4. update state */ - m_prev_mini_info = m_curr_mini_info; - m_prev_mini_info_rc = m_curr_mini_info_rc; m_start = false; assert(equal_lookup_result(m_dict->lookup(kmer), m_res)); @@ -113,6 +106,8 @@ struct streaming_query // uint64_t num_positive_lookups() const { return num_searches() + num_extensions(); } uint64_t num_negative_lookups() const { return m_num_negative; } uint64_t num_invalid_lookups() const { return m_num_invalid; } + uint64_t num_skipped_singleton_lookups() const { return m_num_skipped_singleton_lookups; } + uint64_t num_bucket_cache_hits() const { return m_num_bucket_cache_hits; } private: Dict const* m_dict; @@ -127,58 +122,118 @@ struct streaming_query // /* minimizer state */ minimizer_iterator m_minimizer_it; - minimizer_iterator_rc m_minimizer_it_rc; - minimizer_info m_curr_mini_info, m_prev_mini_info; - minimizer_info m_curr_mini_info_rc, m_prev_mini_info_rc; + minimizer_info m_curr_mini_info; /* string state */ kmer_iterator m_it; uint64_t m_remaining_string_bases; + /* + Last real seed state, for the shortcuts of `seed`: the minimizer of + the last `m_dict->lookup` call (with, for the singleton shortcut, the stream + position of its occurrence and whether a positive match anchors it), + and the cached decoded locate set of that minimizer. + */ + minimizer_info m_last_seed_mini_info; + bool m_last_seed_positive; + /* performance counts */ uint64_t m_num_searches; uint64_t m_num_extensions; uint64_t m_num_invalid; uint64_t m_num_negative; + uint64_t m_num_skipped_singleton_lookups; + uint64_t m_num_bucket_cache_hits; + + /* large and cold: keep it last, away from the per-kmer members above */ + typename Dict::spss_type::bucket_cache m_bucket_cache; - void seed() // + /* Whether the minimizer of the last seed is its own reverse complement + (possible for even m only). Rarely needed, so computed on demand. */ + bool self_rc_minimizer() const { + if constexpr (kmer_t::has_reverse_complement) { + return kmer_t::reverse_complement_mmer(m_last_seed_mini_info.minimizer, m_m) == + m_last_seed_mini_info.minimizer; + } + return false; + } + + /* Outlined: `seed` inlines into `lookup` at two call sites, and letting + its body (grown by the shortcuts) bloat the extension fast path measurably + slows down streams that rarely seed. */ + __attribute__((noinline)) void seed() // { m_remaining_string_bases = 0; - /* if minimizer does not change and previous minimizer was not found, - surely any kmer having the same minimizer cannot be found as well */ - if (m_curr_mini_info.minimizer == m_prev_mini_info.minimizer and // - m_curr_mini_info_rc.minimizer == m_prev_mini_info_rc.minimizer and // - m_res.minimizer_found == false) // + if (m_curr_mini_info.minimizer == m_last_seed_mini_info.minimizer) // { - assert(m_res.kmer_id == constants::invalid_uint64); - m_num_negative += 1; - return; - } + /* The minimizer was absent at the last seed: any kmer having the + same minimizer is surely absent as well. */ + if (m_res.minimizer_found == false) { + assert(m_res.kmer_id == constants::invalid_uint64); + m_num_negative += 1; + return; + } - if constexpr (canonical) { - if (m_curr_mini_info.minimizer < m_curr_mini_info_rc.minimizer) { - m_res = m_dict->lookup_canonical(m_kmer, m_kmer_rc, m_curr_mini_info); - } else if (m_curr_mini_info_rc.minimizer < m_curr_mini_info.minimizer) { - m_res = m_dict->lookup_canonical(m_kmer, m_kmer_rc, m_curr_mini_info_rc); - } else { - m_res = m_dict->lookup_canonical(m_kmer, m_kmer_rc, m_curr_mini_info); + /* + The minimizer is present but the kmer may not be -- the most + common negative in a positive stream, e.g., a sequencing error + inside the kmer. Two shortcuts avoid the full lookup. + + 1. The last seed matched via a singleton bucket, at position L + say, and the current kmer still carries the very same + minimizer occurrence (same `pos_in_seq`: the scheme is + forward, so a sampled position, once abandoned, is never + re-selected). Then the only locus the kmer could occupy is + L - pos_in_kmer, which the just-failed extension step + already compared and rejected (or it lies beyond the string + boundary): negative, with no text access at all. The mirror + locus L - (k - m - pos_in_kmer) hosts the reverse + complement of the minimizer occurrence, not the occurrence + itself, so it is excluded too -- unless the minimizer is + its own reverse complement (possible for even m only), in + which case we fall through to the verification below. + Only a singleton bucket admits this free skip: a larger + bucket has candidate loci the failed extension never + examined, so those must be verified against the text. + */ + if (m_last_seed_positive and m_bucket_cache.size == 1 and + m_curr_mini_info.pos_in_seq == m_last_seed_mini_info.pos_in_seq and + !self_rc_minimizer()) // + { + m_res = lookup_result(); + m_num_negative += 1; + m_num_skipped_singleton_lookups += 1; + return; + } + + /* 2. The decoded locate set of the bucket is cached: verify the + candidate loci directly against the text, with no codeword + lookup nor offset decoding. */ + if (m_bucket_cache.size != 0) { + m_res = m_dict->m_spss.lookup_from_positions( // + m_bucket_cache.positions.data(), m_bucket_cache.size, m_kmer, m_kmer_rc, + m_curr_mini_info); + m_num_bucket_cache_hits += 1; if (m_res.kmer_id == constants::invalid_uint64) { - m_res = m_dict->lookup_canonical(m_kmer, m_kmer_rc, m_curr_mini_info_rc); + m_num_negative += 1; + return; } + /* keep the singleton shortcut anchored to the latest match */ + m_last_seed_mini_info = m_curr_mini_info; + m_last_seed_positive = true; + m_num_searches += 1; + begin_extension(); + return; } - } else { - m_res = m_dict->lookup_regular(m_kmer, m_curr_mini_info); - bool minimizer_found = m_res.minimizer_found; - if (m_res.kmer_id == constants::invalid_uint64) { - assert(m_res.kmer_orientation == constants::forward_orientation); - m_res = m_dict->lookup_regular(m_kmer_rc, m_curr_mini_info_rc); - m_res.kmer_orientation = constants::backward_orientation; - bool minimizer_rc_found = m_res.minimizer_found; - m_res.minimizer_found = minimizer_rc_found or minimizer_found; - } + + /* heavy bucket: fall through to the full lookup */ } + m_res = m_dict->lookup(m_kmer, m_kmer_rc, m_curr_mini_info, &m_bucket_cache); + m_last_seed_mini_info = m_curr_mini_info; + m_last_seed_positive = m_res.kmer_id != constants::invalid_uint64; + if (m_res.kmer_id == constants::invalid_uint64) { m_num_negative += 1; return; @@ -186,6 +241,10 @@ struct streaming_query // assert(m_res.minimizer_found == true); m_num_searches += 1; + begin_extension(); + } + + void begin_extension() { uint64_t kmer_offset = 2 * (m_res.kmer_id + m_res.string_id * (m_k - 1)); m_remaining_string_bases = (m_res.string_end - m_res.string_begin - m_k) - m_res.kmer_id_in_string; diff --git a/include/util.hpp b/include/util.hpp index bf9bebd..3d7ee87 100644 --- a/include/util.hpp +++ b/include/util.hpp @@ -25,7 +25,9 @@ struct streaming_query_report { , num_negative_kmers(0) , num_invalid_kmers(0) , num_searches(0) - , num_extensions(0) {} + , num_extensions(0) + , num_skipped_singleton_lookups(0) + , num_bucket_cache_hits(0) {} uint64_t num_kmers; uint64_t num_positive_kmers; @@ -33,6 +35,12 @@ struct streaming_query_report { uint64_t num_invalid_kmers; uint64_t num_searches; uint64_t num_extensions; + + /* seeds answered from the state of the last seed, saving a full lookup: + negatives implied by a singleton bucket (skipped entirely), and seeds + verified from the cached locate set (see `streaming_query::seed`) */ + uint64_t num_skipped_singleton_lookups; + uint64_t num_bucket_cache_hits; }; struct lookup_result { @@ -74,10 +82,10 @@ inline std::ostream& operator<<(std::ostream& os, lookup_result const& res) { return os; } -template +template struct neighbourhood { - std::array forward; - std::array backward; + std::array forward; + std::array backward; }; struct minimizer_info { @@ -150,7 +158,6 @@ struct build_configuration { , lambda(constants::lambda) - , canonical(false) , weighted(false) , verbose(true) @@ -166,7 +173,6 @@ struct build_configuration { double lambda; // drive PTHash trade-off - bool canonical; bool weighted; bool verbose; @@ -179,7 +185,6 @@ struct build_configuration { << ", num_threads = " << num_threads // << ", ram_limit_in_GiB = " << ram_limit_in_GiB // << ", lambda = " << lambda // - << ", canonical = " << (canonical ? "true" : "false") // << ", weighted = " << (weighted ? "true" : "false") // << ", verbose = " << (verbose ? "true" : "false") // << ", tmp_dirname = '" << tmp_dirname << "'" << std::endl; // @@ -204,82 +209,266 @@ static inline uint64_t get_seed_for_hash_function(build_configuration const& bui return std::equal(pattern.begin(), pattern.end(), str.end() - pattern.size()); } -template -[[maybe_unused]] static kmer_t string_to_uint_kmer(char const* str, uint64_t k) { - assert(k <= kmer_t::max_k); - kmer_t x = 0; - for (uint64_t i = 0; i != k; ++i) x.set(i, kmer_t::char_to_uint(str[i])); +template +[[maybe_unused]] static Kmer string_to_uint_kmer(char const* str, uint64_t k) { + assert(k <= Kmer::max_k); + Kmer x = 0; + for (uint64_t i = 0; i != k; ++i) x.set(i, Kmer::char_to_uint(str[i])); return x; } -template -static void uint_kmer_to_string(kmer_t x, char* str, uint64_t k) { - assert(k <= kmer_t::max_k); - for (uint64_t i = 0; i != k; ++i) str[i] = kmer_t::uint64_to_char(x.pop_char()); +template +static void uint_kmer_to_string(Kmer x, char* str, uint64_t k) { + assert(k <= Kmer::max_k); + for (uint64_t i = 0; i != k; ++i) str[i] = Kmer::uint64_to_char(x.pop_char()); } -template -[[maybe_unused]] static std::string uint_kmer_to_string(kmer_t x, uint64_t k) { - assert(k <= kmer_t::max_k); +template +[[maybe_unused]] static std::string uint_kmer_to_string(Kmer x, uint64_t k) { + assert(k <= Kmer::max_k); std::string str; str.resize(k); uint_kmer_to_string(x, str.data(), k); return str; } -template +template [[maybe_unused]] static std::string uint_minimizer_to_string(uint64_t minimizer, uint64_t m) { - assert(m <= kmer_t::max_m); + assert(m <= Kmer::max_m); std::string str; str.resize(m); - kmer_t x = minimizer; + Kmer x = minimizer; uint_kmer_to_string(x, str.data(), m); return str; } -template +template [[maybe_unused]] static bool is_valid(char const* str, uint64_t size) { for (uint64_t i = 0; i != size; ++i) { - if (!kmer_t::is_valid(str[i])) return false; + if (!Kmer::is_valid(str[i])) return false; } return true; } -template -static kmer_t read_kmer_at(bits::bit_vector const& bv, const uint64_t k, const uint64_t pos) { - static_assert(kmer_t::uint_kmer_bits % 64 == 0); - kmer_t kmer = 0; - for (int i = kmer_t::uint_kmer_bits - 64; i >= 0; i -= 64) { +template +static Kmer read_kmer_at(bits::bit_vector const& bv, const uint64_t k, const uint64_t pos) { + static_assert(Kmer::uint_kmer_bits % 64 == 0); + Kmer kmer = 0; + for (int i = Kmer::uint_kmer_bits - 64; i >= 0; i -= 64) { if (pos + i < bv.num_bits()) kmer.append64(bv.get_word64(pos + i)); } - kmer.take(kmer_t::bits_per_char * k); + kmer.take(Kmer::bits_per_char * k); return kmer; } /* - This implements the random minimizer. + The canonical m-mer: the smaller of the m-mer and its reverse + complement, under the numeric order on the packed encoding. For alphabets + that have no reverse complement (e.g. amino acids) this is the identity. +*/ +template +inline uint64_t canonical_mmer(const uint64_t mmer, const uint64_t m) { + if constexpr (!Kmer::has_reverse_complement) { + (void)m; + return mmer; + } else { + return std::min(mmer, Kmer::reverse_complement_mmer(mmer, m)); + } +} + +/* + The canonical m-mer at locus i of a kmer x, given rc(x). + + Uses rc(x_i) = rc(x)_{k-m-i}: the reverse complements of all k-m+1 loci are + windows of the single reverse complement of the whole kmer, so the kmer is + reverse-complemented once instead of once per locus. `window` is x shifted so + that locus i sits at its low end -- the forward sliding window the callers + already carry along. + + When the whole kmer fits in a word -- always so for the 64-bit kmer type, and + for k <= 32 with the wider one -- the reverse window is extracted with a plain + 64-bit shift; otherwise with a shift of the wide type, which costs a few + instructions more but still beats reverse-complementing every locus. +*/ +template +inline uint64_t canonical_mmer_at(Kmer window, [[maybe_unused]] Kmer const& kmer_rc, + [[maybe_unused]] const uint64_t k, const uint64_t m, + [[maybe_unused]] const uint64_t i) // +{ + window.take_chars(m); + if constexpr (!Kmer::has_reverse_complement) { + return uint64_t(window); + } else { + constexpr uint64_t b = Kmer::bits_per_char; + assert(i + m <= k); + uint64_t rc; + if (k * b <= 64) { + const uint64_t mask = m * b == 64 ? ~uint64_t(0) : (uint64_t(1) << (m * b)) - 1; + rc = (uint64_t(kmer_rc) >> (b * (k - m - i))) & mask; + } else { + Kmer mmer_rc = kmer_rc; + mmer_rc.drop_chars(k - m - i); + mmer_rc.take_chars(m); + rc = uint64_t(mmer_rc); + } + assert(rc == Kmer::reverse_complement_mmer(uint64_t(window), m)); + return std::min(uint64_t(window), rc); + } +} + +/* + True if `kmer` is the canonical one of the pair (kmer, rc(kmer)). A kmer that + equals its own reverse complement (possible only for even k) is deemed + canonical, so that the answer is always well defined. +*/ +template +inline bool is_canonical(Kmer kmer, const uint64_t k) { + if constexpr (!Kmer::has_reverse_complement) { + (void)k; + return true; + } else { + Kmer kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(k); + return !(kmer_rc < kmer); + } +} + +/* + This implements the random minimizer of a kmer x: the locus i* in [0, k-m] + minimizing h(kappa(i)), where kappa(i) := min(x_i, rc(x_i)) is the canonical + m-mer at locus i. The minimizer's value is kappa(i*). + + Because kappa(i) is invariant under reverse-complementation and the sequence + (kappa(0), ..., kappa(k-m)) of rc(x) is the reversal of that of x, the locus + selected for rc(x) is the mirror image of the one selected for x, and the + value is literally the same m-mer. So x and rc(x) always land in the same + bucket -- a lookup costs a single probe -- while the minimizer remains a + minimizer of the window under a single random order, which is what keeps the + density (and hence the number of super-kmers) at that of the plain forward + minimizer, 2/(k-m+2). + + When the alphabet has no reverse complement, kappa is the identity and this + reduces exactly to the plain forward minimizer. + + Ties -- h(kappa(i)) == h(kappa(j)) for i != j, which happens when x_i = x_j + or x_i = rc(x_j) -- must be broken in a mirror-equivariant way, or x and + rc(x) would be sent to different buckets, which is a correctness failure and + not merely a density one. No leftmost rule can be equivariant: the mirror + reverses the order of the tied loci, so "leftmost" names different loci on + the two strands. We use the tie-break of [Cologni and Pibiri, "Canonical + Schemes and Minimizers", Proposition 24]: among the tied loci, take the one + closest to the centre (k-m)/2 of the window -- the one reference point both + strands agree on -- and when two are equally close (they are then mirror + images i and k-m-i), take the smaller index if x <= rc(x) and the larger + otherwise. This rule is mirror-equivariant (for odd k, where x != rc(x) + always) and, unlike the previous frame-of-the-canonical-kmer rule, it is + also *forward*: the distance of a fixed locus to the centre is monotone as + the window slides right, so a locus at least as close as another stays so, + and the anchor never moves backwards. Forwardness makes the number of + super-kmers equal the number of sampled positions -- nothing is ever + abandoned and re-opened. + + Since rc(x_i) = rc(x)_{k-m-i}, the reverse complements of all k-m+1 loci are + just windows of the single reverse complement of the whole kmer: `kmer_rc` is + reverse-complemented once and each rc(x_i) is read out of it with a shift, + instead of reverse-complementing every locus separately. The caller usually + has rc(x) already -- a lookup needs it anyway to recognise which orientation + it found -- in which case the reverse complementation is free. */ -template -minimizer_info compute_minimizer(kmer_t kmer, const uint64_t k, const uint64_t m, +/* + A tie at the minimum (rare, rate ~ sigma^{-m/2}): all tied loci carry the + same class, so only the sampled position is at stake, not the value. Take + the tied locus closest to the centre of the window; the doubled distance + |2i - (k-m)| avoids the half-integer centre of an odd window. + + Deliberately not inlined: this runs on ~1e-4 of the windows, and letting its + loop inline into the lookup path costs measurably more than the call. +*/ +template +__attribute__((noinline)) minimizer_info resolve_tie(Kmer const& kmer, Kmer const& kmer_rc, + const uint64_t k, const uint64_t m, + hasher_type const& hasher, + const uint64_t min_hash, + const uint64_t minimizer, + const uint64_t leftmost, + const uint64_t rightmost) // +{ + const uint64_t two_c = k - m; + uint64_t best_dist = constants::invalid_uint64; + uint64_t lo = leftmost; + uint64_t hi = leftmost; + Kmer window = kmer; + window.drop_chars(leftmost); + for (uint64_t i = leftmost; i <= rightmost; ++i) { + if (hasher.hash(canonical_mmer_at(window, kmer_rc, k, m, i)) == min_hash) { + assert(canonical_mmer_at(window, kmer_rc, k, m, i) == minimizer); + const uint64_t dist = 2 * i > two_c ? 2 * i - two_c : two_c - 2 * i; + if (dist < best_dist) { + best_dist = dist; + lo = i; + hi = i; + } else if (dist == best_dist) { + hi = i; + } + } + window.drop_char(); + } + /* lo and hi are the two equally-closest loci (mirror images), or one locus */ + const uint64_t chosen = (lo == hi or !(kmer_rc < kmer)) ? lo : hi; + return {minimizer, chosen}; +} + +template +minimizer_info compute_minimizer(Kmer kmer, Kmer const& kmer_rc, const uint64_t k, const uint64_t m, hasher_type const& hasher) // { - assert(m <= kmer_t::max_m); + assert(m <= Kmer::max_m); assert(m <= k); - uint64_t min_hash = constants::invalid_uint64; - kmer_t minimizer = kmer_t(-1); - uint64_t pos = 0; - for (uint64_t i = 0; i != k - m + 1; ++i) { - kmer_t mmer = kmer; - mmer.take_chars(m); - uint64_t hash = hasher.hash(uint64_t(mmer)); - if (hash < min_hash) { + + auto kappa = [&](Kmer const& window, const uint64_t i) { + return canonical_mmer_at(window, kmer_rc, k, m, i); + }; + + /* The first locus is peeled off the loop so that `min_hash` starts out at a + real hash value: initializing it to invalid_uint64 would make an actual + hash of invalid_uint64 register as a tie rather than as the minimum. */ + Kmer window = kmer; + uint64_t minimizer = kappa(window, 0); + uint64_t min_hash = hasher.hash(minimizer); + uint64_t leftmost = 0; + uint64_t rightmost = 0; + bool tie = false; + window.drop_char(); + + for (uint64_t i = 1; i != k - m + 1; ++i) { + uint64_t value = kappa(window, i); + uint64_t hash = hasher.hash(value); + if (hash < min_hash) { // leftmost min_hash = hash; - minimizer = mmer; - pos = i; + minimizer = value; + leftmost = i; + rightmost = i; + tie = false; + } else if (hash == min_hash) { + rightmost = i; + tie = true; } - kmer.drop_char(); + window.drop_char(); + } + + if (__builtin_expect(tie, false)) { + return resolve_tie(kmer, kmer_rc, k, m, hasher, min_hash, minimizer, leftmost, rightmost); } - return {uint64_t(minimizer), pos}; + return {minimizer, leftmost}; +} + +template +minimizer_info compute_minimizer(Kmer kmer, const uint64_t k, const uint64_t m, + hasher_type const& hasher) // +{ + Kmer kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(k); + return compute_minimizer(kmer, kmer_rc, k, m, hasher); } } // namespace util diff --git a/script/bench.py b/script/bench.py index 2b334e3..f082b41 100644 --- a/script/bench.py +++ b/script/bench.py @@ -47,19 +47,17 @@ def build_project(max_k63: bool): run_cmd(["make", "-j"]) -def run_bench(k, canonical, runs = 3): +def run_bench(k, runs = 3): """Run SSHASH benchmark for all datasets.""" - mode = "canon" if canonical else "regular" out_dir = results_dir / f"k{k}" out_dir.mkdir(parents=True, exist_ok=True) - log_file = out_dir / f"{mode}-bench.log" - json_file = out_dir / f"{mode}-bench.json" + log_file = out_dir / "bench.log" + json_file = out_dir / "bench.json" for dataset in datasets: - suffix = f".k{k}.canon.sshash" if canonical else f".k{k}.sshash" - index_path = index_dir / f"{dataset}{suffix}" + index_path = index_dir / f"{dataset}.k{k}.sshash" - print(f"\n>>> Benchmarking {dataset} (k={k}, mode={mode})\n") + print(f"\n>>> Benchmarking {dataset} (k={k})\n") for i in range(runs): print(f" ==> run {i+1}/{runs}") cmd = ["./sshash", "bench", "-i", str(index_path)] @@ -82,13 +80,11 @@ def run_bench(k, canonical, runs = 3): # --- Build for k=31 --- build_project(max_k63=False) -run_bench(31, False) -run_bench(31, True) +run_bench(31) # --- Build for k=63 --- build_project(max_k63=True) -run_bench(63, False) -run_bench(63, True) +run_bench(63) # --- Restore to default --- build_project(max_k63=False) diff --git a/script/build.py b/script/build.py index c80973b..8964a79 100644 --- a/script/build.py +++ b/script/build.py @@ -53,23 +53,20 @@ def build_project(max_k63: bool): run_cmd(["make", "-j"]) -def build_sshash(k, canonical, m_values): - mode_dir = results_dir / f"k{k}" - mode_dir.mkdir(parents=True, exist_ok=True) +def build_sshash(k, m_values): + out_dir = results_dir / f"k{k}" + out_dir.mkdir(parents=True, exist_ok=True) - mode = "canon" if canonical else "regular" - log_file = mode_dir / f"{mode}-build.log" - json_file = mode_dir / f"{mode}-build.json" - time_file = mode_dir / f"{mode}-build.time.log" + log_file = out_dir / "build.log" + json_file = out_dir / "build.json" + time_file = out_dir / "build.time.log" for dataset in datasets: m_val = m_values[dataset] input_file = datasets_dir / f"{dataset}.k{k}.eulertigs.fa.gz" output_file = index_dir / f"{dataset}.k{k}" - if canonical: - output_file = str(output_file) + ".canon" - print(f"\n>>> Building {dataset} (k={k}, m={m_val}, mode={mode})\n") + print(f"\n>>> Building {dataset} (k={k}, m={m_val})\n") # Clean tmp directory (should be empty after each build anyway) subprocess.run(f"rm -rf {tmp_dir}/*", shell=True, check=True) @@ -86,8 +83,6 @@ def build_sshash(k, canonical, m_values): "-d", str(tmp_dir), "-o", f"{output_file}.sshash" ] - if canonical: - cmd.append("--canonical") # Append stdout to .log, stderr to .json with open(log_file, "a") as log, open(json_file, "a") as js: @@ -100,13 +95,11 @@ def build_sshash(k, canonical, m_values): # k = 31 build_project(max_k63=False) -build_sshash(31, False, m_values_k31) -build_sshash(31, True, m_values_k31) +build_sshash(31, m_values_k31) # k = 63 build_project(max_k63=True) -build_sshash(63, False, m_values_k63) -build_sshash(63, True, m_values_k63) +build_sshash(63, m_values_k63) # rebuild back to default build_project(max_k63=False) diff --git a/script/plot-trade-off-l.py b/script/plot-trade-off-l.py index 8b56ea2..424f0c6 100644 --- a/script/plot-trade-off-l.py +++ b/script/plot-trade-off-l.py @@ -32,102 +32,96 @@ def parse_results(results_dir): if not k_match: continue k_val = int(k_match.group(1)) - for mode in ['regular', 'canon']: - build_json = k_dir / f"{mode}-build.json" - bench_json = k_dir / f"{mode}-bench.json" - - if not build_json.exists() or not bench_json.exists(): - continue - - # 1. Parse Build JSON for Space (bits/k-mer) - space_dict = {} - with open(build_json, 'r') as f: - for line in f: - try: - j = json.loads(line) - - # Extract dataset name from "/mnt/.../human.k31.eulertigs.fa.gz" - filename = os.path.basename(j.get("input_filename", "")) - ds = filename.split('.')[0] if filename else "unknown" - - if "index_size_in_bytes" in j and "num_kmers" in j: - bytes_size = float(j["index_size_in_bytes"]) - num_kmers = float(j["num_kmers"]) - # Calculate bits per k-mer - bits_per_kmer = (bytes_size * 8.0) / num_kmers - space_dict[ds] = bits_per_kmer - except json.JSONDecodeError: - continue - - # 2. Parse Bench JSON for Query Time (ns/kmer) - time_dict = {} - count_dict = {} - with open(bench_json, 'r') as f: - for line in f: - try: - j = json.loads(line) - - # Extract dataset name from ".../human.k31.l4.sshash" - filename = os.path.basename(j.get("index_filename", "")) - ds = filename.split('.')[0] if filename else "unknown" - - # Use positive lookup time - t_str = j.get("positive lookup (avg_nanosec_per_kmer)") - if t_str is not None: - t = float(t_str) - time_dict[ds] = time_dict.get(ds, 0.0) + t - count_dict[ds] = count_dict.get(ds, 0) + 1 - except json.JSONDecodeError: - continue - - # 3. Combine and store - for ds in space_dict.keys(): - if ds in time_dict and count_dict[ds] > 0: - # Average the benchmark runs - avg_time = time_dict[ds] / count_dict[ds] - data.append({ - 'Dataset': ds, - 'k': k_val, - 'l': l_val, - 'Mode': mode, - 'Space (bits/k-mer)': space_dict[ds], - 'Query Time (ns/k-mer)': avg_time - }) + build_json = k_dir / "build.json" + bench_json = k_dir / "bench.json" + + if not build_json.exists() or not bench_json.exists(): + continue + + # 1. Parse Build JSON for Space (bits/k-mer) + space_dict = {} + with open(build_json, 'r') as f: + for line in f: + try: + j = json.loads(line) + + # Extract dataset name from "/mnt/.../human.k31.eulertigs.fa.gz" + filename = os.path.basename(j.get("input_filename", "")) + ds = filename.split('.')[0] if filename else "unknown" + + if "index_size_in_bytes" in j and "num_kmers" in j: + bytes_size = float(j["index_size_in_bytes"]) + num_kmers = float(j["num_kmers"]) + # Calculate bits per k-mer + bits_per_kmer = (bytes_size * 8.0) / num_kmers + space_dict[ds] = bits_per_kmer + except json.JSONDecodeError: + continue + + # 2. Parse Bench JSON for Query Time (ns/kmer) + time_dict = {} + count_dict = {} + with open(bench_json, 'r') as f: + for line in f: + try: + j = json.loads(line) + + # Extract dataset name from ".../human.k31.l4.sshash" + filename = os.path.basename(j.get("index_filename", "")) + ds = filename.split('.')[0] if filename else "unknown" + + # Use positive lookup time + t_str = j.get("positive lookup (avg_nanosec_per_kmer)") + if t_str is not None: + t = float(t_str) + time_dict[ds] = time_dict.get(ds, 0.0) + t + count_dict[ds] = count_dict.get(ds, 0) + 1 + except json.JSONDecodeError: + continue + + # 3. Combine and store + for ds in space_dict.keys(): + if ds in time_dict and count_dict[ds] > 0: + # Average the benchmark runs + avg_time = time_dict[ds] / count_dict[ds] + data.append({ + 'Dataset': ds, + 'k': k_val, + 'l': l_val, + 'Space (bits/k-mer)': space_dict[ds], + 'Query Time (ns/k-mer)': avg_time + }) return pd.DataFrame(data) def plot_tradeoff(df, output_img="tradeoff_plot_l.pdf"): """ Generates a space-time trade-off plot. - Different lines for datasets/k/modes, points vary by 'l'. + Different lines for datasets/k, points vary by 'l'. """ if df.empty: print("No data parsed! Please check the JSON keys in the script.") return # Enforce categorical order so the legend is populated exactly how we want: - # Human before SE, and regular before canon + # Human before SE df['Dataset'] = pd.Categorical(df['Dataset'], categories=['human', 'se'], ordered=True) - df['Mode'] = pd.Categorical(df['Mode'], categories=['regular', 'canon'], ordered=True) - df = df.sort_values(by=['Dataset', 'k', 'Mode']) + df = df.sort_values(by=['Dataset', 'k']) plt.figure(figsize=(5, 10)) plt.style.use('seaborn-v0_8-whitegrid') - # Group by Dataset, k, and Mode with sort=False to preserve our categorical ordering - groups = df.groupby(['Dataset', 'k', 'Mode'], sort=False) + # Group by Dataset and k with sort=False to preserve our categorical ordering + groups = df.groupby(['Dataset', 'k'], sort=False) - for (dataset, k, mode), group in groups: + for (dataset, k), group in groups: # Sort by l to make the line connect logically group = group.sort_values(by='l') - label = f"{'Human' if dataset == 'human' else 'SE'} (k={k}, {mode})" - - # Color logic: Red for Human, Blue for SE. Darker if canonical. - if dataset == 'human': - color = 'firebrick' if mode == 'canon' else 'lightcoral' - else: # se - color = 'royalblue' if mode == 'canon' else 'lightskyblue' + label = f"{'Human' if dataset == 'human' else 'SE'} (k={k})" + + # Color logic: Red for Human, Blue for SE. + color = 'firebrick' if dataset == 'human' else 'royalblue' # Marker logic: Circle for k=31, Square for k=63 marker = 'o' if k == 31 else 's' diff --git a/script/plot-trade-off-m.py b/script/plot-trade-off-m.py index 0a1d079..cfffbee 100644 --- a/script/plot-trade-off-m.py +++ b/script/plot-trade-off-m.py @@ -32,102 +32,96 @@ def parse_results(results_dir): if not k_match: continue k_val = int(k_match.group(1)) - for mode in ['regular', 'canon']: - build_json = k_dir / f"{mode}-build.json" - bench_json = k_dir / f"{mode}-bench.json" - - if not build_json.exists() or not bench_json.exists(): - continue - - # 1. Parse Build JSON for Space (bits/k-mer) - space_dict = {} - with open(build_json, 'r') as f: - for line in f: - try: - j = json.loads(line) - - # Extract dataset name from "/mnt/.../human.k31.eulertigs.fa.gz" - filename = os.path.basename(j.get("input_filename", "")) - ds = filename.split('.')[0] if filename else "unknown" - - if "index_size_in_bytes" in j and "num_kmers" in j: - bytes_size = float(j["index_size_in_bytes"]) - num_kmers = float(j["num_kmers"]) - # Calculate bits per k-mer - bits_per_kmer = (bytes_size * 8.0) / num_kmers - space_dict[ds] = bits_per_kmer - except json.JSONDecodeError: - continue - - # 2. Parse Bench JSON for Query Time (ns/kmer) - time_dict = {} - count_dict = {} - with open(bench_json, 'r') as f: - for line in f: - try: - j = json.loads(line) - - # Extract dataset name from ".../human.k31.m17.sshash" - filename = os.path.basename(j.get("index_filename", "")) - ds = filename.split('.')[0] if filename else "unknown" - - # Use positive lookup time - t_str = j.get("positive lookup (avg_nanosec_per_kmer)") - if t_str is not None: - t = float(t_str) - time_dict[ds] = time_dict.get(ds, 0.0) + t - count_dict[ds] = count_dict.get(ds, 0) + 1 - except json.JSONDecodeError: - continue - - # 3. Combine and store - for ds in space_dict.keys(): - if ds in time_dict and count_dict[ds] > 0: - # Average the 3 benchmark runs - avg_time = time_dict[ds] / count_dict[ds] - data.append({ - 'Dataset': ds, - 'k': k_val, - 'm': m_val, - 'Mode': mode, - 'Space (bits/k-mer)': space_dict[ds], - 'Query Time (ns/k-mer)': avg_time - }) + build_json = k_dir / "build.json" + bench_json = k_dir / "bench.json" + + if not build_json.exists() or not bench_json.exists(): + continue + + # 1. Parse Build JSON for Space (bits/k-mer) + space_dict = {} + with open(build_json, 'r') as f: + for line in f: + try: + j = json.loads(line) + + # Extract dataset name from "/mnt/.../human.k31.eulertigs.fa.gz" + filename = os.path.basename(j.get("input_filename", "")) + ds = filename.split('.')[0] if filename else "unknown" + + if "index_size_in_bytes" in j and "num_kmers" in j: + bytes_size = float(j["index_size_in_bytes"]) + num_kmers = float(j["num_kmers"]) + # Calculate bits per k-mer + bits_per_kmer = (bytes_size * 8.0) / num_kmers + space_dict[ds] = bits_per_kmer + except json.JSONDecodeError: + continue + + # 2. Parse Bench JSON for Query Time (ns/kmer) + time_dict = {} + count_dict = {} + with open(bench_json, 'r') as f: + for line in f: + try: + j = json.loads(line) + + # Extract dataset name from ".../human.k31.m17.sshash" + filename = os.path.basename(j.get("index_filename", "")) + ds = filename.split('.')[0] if filename else "unknown" + + # Use positive lookup time + t_str = j.get("positive lookup (avg_nanosec_per_kmer)") + if t_str is not None: + t = float(t_str) + time_dict[ds] = time_dict.get(ds, 0.0) + t + count_dict[ds] = count_dict.get(ds, 0) + 1 + except json.JSONDecodeError: + continue + + # 3. Combine and store + for ds in space_dict.keys(): + if ds in time_dict and count_dict[ds] > 0: + # Average the 3 benchmark runs + avg_time = time_dict[ds] / count_dict[ds] + data.append({ + 'Dataset': ds, + 'k': k_val, + 'm': m_val, + 'Space (bits/k-mer)': space_dict[ds], + 'Query Time (ns/k-mer)': avg_time + }) return pd.DataFrame(data) def plot_tradeoff(df, output_img="tradeoff_plot.pdf"): """ Generates a space-time trade-off plot. - Different lines for datasets/k/modes, points vary by 'm'. + Different lines for datasets/k, points vary by 'm'. """ if df.empty: print("No data parsed! Please check the JSON keys in the script.") return # Enforce categorical order so the legend is populated exactly how we want: - # Human before SE, and regular before canon + # Human before SE df['Dataset'] = pd.Categorical(df['Dataset'], categories=['human', 'se'], ordered=True) - df['Mode'] = pd.Categorical(df['Mode'], categories=['regular', 'canon'], ordered=True) - df = df.sort_values(by=['Dataset', 'k', 'Mode']) + df = df.sort_values(by=['Dataset', 'k']) plt.figure(figsize=(5, 10)) plt.style.use('seaborn-v0_8-whitegrid') - # Group by Dataset, k, and Mode with sort=False to preserve our categorical ordering - groups = df.groupby(['Dataset', 'k', 'Mode'], sort=False) + # Group by Dataset and k with sort=False to preserve our categorical ordering + groups = df.groupby(['Dataset', 'k'], sort=False) - for (dataset, k, mode), group in groups: + for (dataset, k), group in groups: # Sort by m to make the line connect logically group = group.sort_values(by='m') - label = f"{'Human' if dataset == 'human' else 'SE'} (k={k}, {mode})" + label = f"{'Human' if dataset == 'human' else 'SE'} (k={k})" - # Color logic: Red for Human, Blue for SE. Darker if canonical. - if dataset == 'human': - color = 'firebrick' if mode == 'canon' else 'lightcoral' - else: # se - color = 'royalblue' if mode == 'canon' else 'lightskyblue' + # Color logic: Red for Human, Blue for SE. + color = 'firebrick' if dataset == 'human' else 'royalblue' # Marker logic: Circle for k=31, Square for k=63 marker = 'o' if k == 31 else 's' diff --git a/script/streaming-query.py b/script/streaming-query.py index e62dc08..5776e3b 100644 --- a/script/streaming-query.py +++ b/script/streaming-query.py @@ -57,19 +57,17 @@ def build_project(max_k63: bool): run_cmd(["make", "-j"]) -def run_bench(k, canonical, runs = 1): +def run_bench(k, runs = 1): """Run SSHASH benchmark for all datasets.""" - mode = "canon" if canonical else "regular" out_dir = results_dir / f"k{k}" out_dir.mkdir(parents=True, exist_ok=True) - log_file = out_dir / f"{mode}-streaming-queries.log" - json_file = out_dir / f"{mode}-streaming-queries.json" + log_file = out_dir / "streaming-queries.log" + json_file = out_dir / "streaming-queries.json" for dataset in datasets: - suffix = f".k{k}.canon.sshash" if canonical else f".k{k}.sshash" - index_path = index_dir / f"{dataset}{suffix}" + index_path = index_dir / f"{dataset}.k{k}.sshash" - print(f"\n>>> Benchmarking {dataset} (k={k}, mode={mode})\n") + print(f"\n>>> Benchmarking {dataset} (k={k})\n") for i in range(runs): print(f" ==> run {i+1}/{runs}") cmd = ["./sshash", "query", "-i", str(index_path), "-q", str(query_dir) + "/" + queries[dataset] + ".fastq.gz"] @@ -92,13 +90,11 @@ def run_bench(k, canonical, runs = 1): # --- Build for k=31 --- build_project(max_k63=False) -run_bench(31, False) -run_bench(31, True) +run_bench(31) # --- Build for k=63 --- build_project(max_k63=True) -run_bench(63, False) -run_bench(63, True) +run_bench(63) # --- Restore to default --- build_project(max_k63=False) diff --git a/script/sweep-m.py b/script/sweep-m.py index 2110f92..8c411be 100644 --- a/script/sweep-m.py +++ b/script/sweep-m.py @@ -53,24 +53,21 @@ def build_project(max_k63: bool): ]) run_cmd(["make", "-j"]) -def build_sshash(k, canonical, dataset, m_val): +def build_sshash(k, dataset, m_val): # Differentiate results dir by m_val - mode_dir = results_dir / f"m{m_val}" / f"k{k}" - mode_dir.mkdir(parents=True, exist_ok=True) + out_dir = results_dir / f"m{m_val}" / f"k{k}" + out_dir.mkdir(parents=True, exist_ok=True) - mode = "canon" if canonical else "regular" - log_file = mode_dir / f"{mode}-build.log" - json_file = mode_dir / f"{mode}-build.json" - time_file = mode_dir / f"{mode}-build.time.log" + log_file = out_dir / "build.log" + json_file = out_dir / "build.json" + time_file = out_dir / "build.time.log" input_file = datasets_dir / f"{dataset}.k{k}.eulertigs.fa.gz" - + # Append m_val to the output filename output_file = index_dir / f"{dataset}.k{k}.m{m_val}" - if canonical: - output_file = str(output_file) + ".canon" - print(f"\n>>> Building {dataset} (k={k}, m={m_val}, mode={mode})\n") + print(f"\n>>> Building {dataset} (k={k}, m={m_val})\n") # Clean tmp directory (should be empty after each build anyway) subprocess.run(f"rm -rf {tmp_dir}/*", shell=True, check=True) @@ -87,28 +84,23 @@ def build_sshash(k, canonical, dataset, m_val): "-d", str(tmp_dir), "-o", f"{output_file}.sshash" ] - if canonical: - cmd.append("--canonical") # Append stdout to .log, stderr to .json with open(log_file, "a") as log, open(json_file, "a") as js: subprocess.run(cmd, stdout=log, stderr=js, check=True) -def run_bench(k, canonical, dataset, m_val, runs=3): +def run_bench(k, dataset, m_val, runs=3): """Run SSHASH benchmark for a specific dataset and m.""" - mode = "canon" if canonical else "regular" - # Store results in the specific m_val / k folder out_dir = results_dir / f"m{m_val}" / f"k{k}" out_dir.mkdir(parents=True, exist_ok=True) - log_file = out_dir / f"{mode}-bench.log" - json_file = out_dir / f"{mode}-bench.json" + log_file = out_dir / "bench.log" + json_file = out_dir / "bench.json" - # Match the new index naming scheme that includes m_val - suffix = f".k{k}.m{m_val}.canon.sshash" if canonical else f".k{k}.m{m_val}.sshash" - index_path = index_dir / f"{dataset}{suffix}" + # Match the index naming scheme that includes m_val + index_path = index_dir / f"{dataset}.k{k}.m{m_val}.sshash" - print(f"\n>>> Benchmarking {dataset} (k={k}, m={m_val}, mode={mode})\n") + print(f"\n>>> Benchmarking {dataset} (k={k}, m={m_val})\n") for i in range(runs): print(f" ==> run {i+1}/{runs}") cmd = ["./sshash", "bench", "-i", str(index_path)] @@ -131,26 +123,16 @@ def run_bench(k, canonical, dataset, m_val, runs=3): for dataset in datasets: for current_m in m_sweeps_k31[dataset]: - # Regular - build_sshash(31, False, dataset, current_m) - run_bench(31, False, dataset, current_m) - - # Canonical - build_sshash(31, True, dataset, current_m) - run_bench(31, True, dataset, current_m) + build_sshash(31, dataset, current_m) + run_bench(31, dataset, current_m) # --- k = 63 Sweep --- build_project(max_k63=True) for dataset in datasets: for current_m in m_sweeps_k63[dataset]: - # Regular - build_sshash(63, False, dataset, current_m) - run_bench(63, False, dataset, current_m) - - # Canonical - build_sshash(63, True, dataset, current_m) - run_bench(63, True, dataset, current_m) + build_sshash(63, dataset, current_m) + run_bench(63, dataset, current_m) # Restore default compilation at the end print("\nRestoring default compilation (max_k63=False)...") diff --git a/script/sweep-min-l.py b/script/sweep-min-l.py index f16729a..c2bf862 100644 --- a/script/sweep-min-l.py +++ b/script/sweep-min-l.py @@ -87,26 +87,23 @@ def build_project(max_k63: bool): ]) run_cmd(["make", "-j"]) -def build_sshash(k, canonical, m_values, l_val): +def build_sshash(k, m_values, l_val): # Differentiate results dir by l_val - mode_dir = results_dir / f"l{l_val}" / f"k{k}" - mode_dir.mkdir(parents=True, exist_ok=True) + out_dir = results_dir / f"l{l_val}" / f"k{k}" + out_dir.mkdir(parents=True, exist_ok=True) - mode = "canon" if canonical else "regular" - log_file = mode_dir / f"{mode}-build.log" - json_file = mode_dir / f"{mode}-build.json" - time_file = mode_dir / f"{mode}-build.time.log" + log_file = out_dir / "build.log" + json_file = out_dir / "build.json" + time_file = out_dir / "build.time.log" for dataset in datasets: m_val = m_values[dataset] input_file = datasets_dir / f"{dataset}.k{k}.eulertigs.fa.gz" - + # 2. Append l_val to the output filename output_file = index_dir / f"{dataset}.k{k}.l{l_val}" - if canonical: - output_file = str(output_file) + ".canon" - print(f"\n>>> Building {dataset} (k={k}, m={m_val}, l={l_val}, mode={mode})\n") + print(f"\n>>> Building {dataset} (k={k}, m={m_val}, l={l_val})\n") # Clean tmp directory (should be empty after each build anyway) subprocess.run(f"rm -rf {tmp_dir}/*", shell=True, check=True) @@ -123,29 +120,24 @@ def build_sshash(k, canonical, m_values, l_val): "-d", str(tmp_dir), "-o", f"{output_file}.sshash" ] - if canonical: - cmd.append("--canonical") # Append stdout to .log, stderr to .json with open(log_file, "a") as log, open(json_file, "a") as js: subprocess.run(cmd, stdout=log, stderr=js, check=True) -def run_bench(k, canonical, l_val, runs=3): +def run_bench(k, l_val, runs=3): """Run SSHASH benchmark for all datasets.""" - mode = "canon" if canonical else "regular" - # Store results in the specific l_val / k folder out_dir = results_dir / f"l{l_val}" / f"k{k}" out_dir.mkdir(parents=True, exist_ok=True) - log_file = out_dir / f"{mode}-bench.log" - json_file = out_dir / f"{mode}-bench.json" + log_file = out_dir / "bench.log" + json_file = out_dir / "bench.json" for dataset in datasets: - # Match the new index naming scheme that includes l_val - suffix = f".k{k}.l{l_val}.canon.sshash" if canonical else f".k{k}.l{l_val}.sshash" - index_path = index_dir / f"{dataset}{suffix}" + # Match the index naming scheme that includes l_val + index_path = index_dir / f"{dataset}.k{k}.l{l_val}.sshash" - print(f"\n>>> Benchmarking {dataset} (k={k}, l={l_val}, mode={mode})\n") + print(f"\n>>> Benchmarking {dataset} (k={k}, l={l_val})\n") for i in range(runs): print(f" ==> run {i+1}/{runs}") cmd = ["./sshash", "bench", "-i", str(index_path)] @@ -166,23 +158,15 @@ def run_bench(k, canonical, l_val, runs=3): # Update the header file update_constants_hpp(current_l) - # Build and benchmark for k = 31 (Regular) + # Build and benchmark for k = 31 build_project(max_k63=False) - build_sshash(31, False, m_values_k31, current_l) - run_bench(31, False, current_l) - - # Build and benchmark for k = 31 (Canonical) - build_sshash(31, True, m_values_k31, current_l) - run_bench(31, True, current_l) + build_sshash(31, m_values_k31, current_l) + run_bench(31, current_l) - # Build and benchmark for k = 63 (Regular) + # Build and benchmark for k = 63 build_project(max_k63=True) - build_sshash(63, False, m_values_k63, current_l) - run_bench(63, False, current_l) - - # Build and benchmark for k = 63 (Canonical) - build_sshash(63, True, m_values_k63, current_l) - run_bench(63, True, current_l) + build_sshash(63, m_values_k63, current_l) + run_bench(63, current_l) # Restore default constants file at the end print("\nRestoring default constants.hpp (min_l = 6)...") diff --git a/src/builder/build_sparse_and_skew_index.cpp b/src/builder/build_sparse_and_skew_index.cpp index 7ed9886..2ef51dc 100644 --- a/src/builder/build_sparse_and_skew_index.cpp +++ b/src/builder/build_sparse_and_skew_index.cpp @@ -21,7 +21,6 @@ void dictionary_builder::build_sparse_and_skew_index( uint64_t num_buckets_larger_than_1_not_in_skew_index = 0; uint64_t num_buckets_in_skew_index = 0; - uint64_t num_super_kmers_in_buckets_larger_than_1 = 0; uint64_t num_minimizer_positions_of_buckets_larger_than_1 = 0; uint64_t num_minimizer_positions_of_buckets_in_skew_index = 0; @@ -41,7 +40,6 @@ void dictionary_builder::build_sparse_and_skew_index( ++num_buckets_in_skew_index; num_minimizer_positions_of_buckets_in_skew_index += bucket_size; } - num_super_kmers_in_buckets_larger_than_1 += bucket.num_super_kmers(); } for (auto mt : bucket) { @@ -97,7 +95,8 @@ void dictionary_builder::build_sparse_and_skew_index( std::vector buckets; buckets.reserve(num_buckets_larger_than_1_not_in_skew_index + num_buckets_in_skew_index); std::vector tuples; // backed memory - tuples.reserve(num_super_kmers_in_buckets_larger_than_1); + tuples.reserve(num_minimizer_positions_of_buckets_larger_than_1 + + num_minimizer_positions_of_buckets_in_skew_index); // Second pass: collect buckets > 1 for sorting AND handle size-1 buckets for (minimizers_tuples_iterator it(input.data(), input.data() + input.size()); // @@ -108,20 +107,13 @@ void dictionary_builder::build_sparse_and_skew_index( const uint64_t bucket_size = bucket.size(); if (bucket_size == 1) { - // Handle size-1 buckets: encode directly into control codewords - uint64_t prev_pos_in_seq = constants::invalid_uint64; - for (auto mt : bucket) { - if (mt.pos_in_seq != prev_pos_in_seq) { - /* - For minimizers occurring once, store a (log(N)+1)-bit - code, as follows: |offset|0|, i.e., the LSB is 0. - */ - uint64_t code = mt.pos_in_seq << 1; // first LS bit encodes status code: 0 - assert(code < (uint64_t(1) << num_bits_for_control)); - control_codewords_builder.set(bucket_id, code); - prev_pos_in_seq = mt.pos_in_seq; - } - } + /* + For minimizers occurring once, store a (log(N)+1)-bit + code, as follows: |offset|0|, i.e., the LSB is 0. + */ + const uint64_t code = (*bucket.begin()).pos_in_seq << 1; // LSB encodes status: 0 + assert(code < (uint64_t(1) << num_bits_for_control)); + control_codewords_builder.set(bucket_id, code); } else { // Collect buckets > 1 for later processing minimizer_tuple const* begin = tuples.data() + tuples.size(); @@ -202,37 +194,21 @@ void dictionary_builder::build_sparse_and_skew_index( } if (curr_bucket_size <= min_size) { - uint64_t prev_pos_in_seq = constants::invalid_uint64; - for (auto mt : bucket) { - if (prev_pos_in_seq == constants::invalid_uint64) { // only once - uint64_t p = (list_id << constants::min_l) | (curr_bucket_size - 2); - uint64_t code = (p << 2) | 1; // first two LS bits encode status code: 01 - assert(code < (uint64_t(1) << num_bits_for_control)); - control_codewords_builder.set(mt.minimizer, code); - } - if (mt.pos_in_seq != prev_pos_in_seq) { - mid_load_buckets_builder.push_back(mt.pos_in_seq); - prev_pos_in_seq = mt.pos_in_seq; - mid_load_buckets_size += 1; - } - } + uint64_t p = (list_id << constants::min_l) | (curr_bucket_size - 2); + uint64_t code = (p << 2) | 1; // first two LS bits encode status code: 01 + assert(code < (uint64_t(1) << num_bits_for_control)); + control_codewords_builder.set(bucket.begin_ptr()->minimizer, code); + for (auto mt : bucket) mid_load_buckets_builder.push_back(mt.pos_in_seq); + mid_load_buckets_size += bucket_size; ++list_id; } else { - uint64_t prev_pos_in_seq = constants::invalid_uint64; - for (auto mt : bucket) { - if (prev_pos_in_seq == constants::invalid_uint64) { // only once - assert(partition_id < 8); - uint64_t p = (heavy_load_buckets_size << 3) | partition_id; - uint64_t code = (p << 2) | 3; // first two LS bits encode status code: 11 - assert(code < (uint64_t(1) << num_bits_for_control)); - control_codewords_builder.set(mt.minimizer, code); - } - if (mt.pos_in_seq != prev_pos_in_seq) { - heavy_load_buckets_builder.push_back(mt.pos_in_seq); - prev_pos_in_seq = mt.pos_in_seq; - heavy_load_buckets_size += 1; - } - } + assert(partition_id < 8); + uint64_t p = (heavy_load_buckets_size << 3) | partition_id; + uint64_t code = (p << 2) | 3; // first two LS bits encode status code: 11 + assert(code < (uint64_t(1) << num_bits_for_control)); + control_codewords_builder.set(bucket.begin_ptr()->minimizer, code); + for (auto mt : bucket) heavy_load_buckets_builder.push_back(mt.pos_in_seq); + heavy_load_buckets_size += bucket_size; } } @@ -443,13 +419,9 @@ void dictionary_builder::build_sparse_and_skew_index( assert(bucket.size() > lower and bucket.size() <= upper); uint64_t pos_in_bucket = -1; - uint64_t prev_pos_in_seq = constants::invalid_uint64; for (auto mt : bucket) // { - if (mt.pos_in_seq != prev_pos_in_seq) { - prev_pos_in_seq = mt.pos_in_seq; - ++pos_in_bucket; - } + ++pos_in_bucket; // one position per tuple: the scheme is forward assert(mt.pos_in_seq >= mt.pos_in_kmer); mt.pos_in_seq = d.m_spss.strings_offsets.decode(mt.pos_in_seq).absolute_offset; @@ -459,12 +431,11 @@ void dictionary_builder::build_sparse_and_skew_index( d.m_spss.strings, k, Kmer::bits_per_char * starting_pos_of_super_kmer); for (uint64_t i = 0; i != mt.num_kmers_in_super_kmer; ++i) { auto kmer = it.get(); - if (build_config.canonical) { /* take the canonical kmer */ - auto kmer_rc = kmer; - kmer_rc.reverse_complement_inplace(k); - kmer = std::min(kmer, kmer_rc); - } - kmers.push_back(kmer); + /* take the canonical kmer: a kmer and its reverse complement + share a bucket, hence a skew-index partition */ + auto kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(k); + kmers.push_back(std::min(kmer, kmer_rc)); positions_in_bucket.push_back(pos_in_bucket); it.next(); } diff --git a/src/builder/compute_minimizer_tuples.cpp b/src/builder/compute_minimizer_tuples.cpp index 94916d6..4070c57 100644 --- a/src/builder/compute_minimizer_tuples.cpp +++ b/src/builder/compute_minimizer_tuples.cpp @@ -28,14 +28,9 @@ void dictionary_builder::compute_minimizer_tuples() // uint64_t num_kmers_in_super_kmer) // { assert(num_kmers_in_super_kmer <= k - m + 1 /* max num kmers in super-kmer */); - if (!buffer.empty() and // - buffer.back().minimizer == mini_info.minimizer and // - buffer.back().pos_in_seq == mini_info.pos_in_seq and // - buffer.back().pos_in_kmer == mini_info.pos_in_kmer) // - { - buffer.back().num_kmers_in_super_kmer += num_kmers_in_super_kmer; - return; - } + /* the scheme is forward: a minimizer position is never saved twice */ + assert(buffer.empty() or !(buffer.back().minimizer == mini_info.minimizer and + buffer.back().pos_in_seq == mini_info.pos_in_seq)); if (buffer.size() == buffer_size) { minimizers.sort_and_flush(buffer); buffer.clear(); @@ -50,7 +45,6 @@ void dictionary_builder::compute_minimizer_tuples() // kmer_iterator kmer_it(strings_builder, k); hasher_type hasher(build_config.seed); minimizer_iterator minimizer_it(k, m, hasher); - minimizer_iterator_rc minimizer_it_rc(k, m, hasher); for (uint64_t i = index_begin; i < index_end; ++i) // { @@ -65,7 +59,6 @@ void dictionary_builder::compute_minimizer_tuples() // kmer_it.at(Kmer::bits_per_char * begin); minimizer_it.set_position(begin); - minimizer_it_rc.set_position(begin); for (uint64_t j = 0; j != sequence_len - k + 1; ++j) { auto uint_kmer = kmer_it.get(); @@ -73,18 +66,6 @@ void dictionary_builder::compute_minimizer_tuples() // assert(mini_info.pos_in_seq < end - m + 1); assert(mini_info.pos_in_kmer < k - m + 1); - if (build_config.canonical) { - auto uint_kmer_rc = uint_kmer; - uint_kmer_rc.reverse_complement_inplace(k); - auto mini_info_rc = minimizer_it_rc.next(uint_kmer_rc); - assert(mini_info_rc.pos_in_seq < end - m + 1); - assert(mini_info_rc.pos_in_kmer < k - m + 1); - if (mini_info_rc.minimizer < mini_info.minimizer) { - mini_info = mini_info_rc; - mini_info.pos_in_kmer = k - m - mini_info.pos_in_kmer; - } - } - mini_info.pos_in_seq = strings_offsets_builder.encode(mini_info.pos_in_seq, begin, i); @@ -92,8 +73,13 @@ void dictionary_builder::compute_minimizer_tuples() // prev_mini_info = mini_info; } - if (mini_info.minimizer != prev_mini_info.minimizer or - mini_info.pos_in_seq != prev_mini_info.pos_in_seq) // + /* + The minimizer's value is the canonical m-mer at the + anchored locus, hence a function of `pos_in_seq` alone: + comparing positions is enough to detect a super-kmer + break, no need to also compare values. + */ + if (mini_info.pos_in_seq != prev_mini_info.pos_in_seq) // { save(prev_mini_info, num_kmers_in_super_kmer); prev_mini_info = mini_info; diff --git a/src/dictionary.cpp b/src/dictionary.cpp index c1b43f7..711d2ca 100644 --- a/src/dictionary.cpp +++ b/src/dictionary.cpp @@ -4,57 +4,6 @@ namespace sshash { -template -lookup_result dictionary::lookup_regular(const Kmer uint_kmer) const { - auto mini_info = util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher); - return lookup_regular(uint_kmer, mini_info); -} - -template -lookup_result dictionary::lookup_regular(const Kmer uint_kmer, // - const minimizer_info mini_info) const // -{ - assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == - util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher)); - - auto it = m_ssi.lookup(uint_kmer, mini_info); - return m_spss.lookup_regular(it, uint_kmer, mini_info); -} - -template -lookup_result dictionary::lookup_canonical(Kmer uint_kmer) const // -{ - Kmer uint_kmer_rc = uint_kmer; - uint_kmer_rc.reverse_complement_inplace(m_k); - auto mini_info = util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher); - auto mini_info_rc = util::compute_minimizer(uint_kmer_rc, m_k, m_m, m_hasher); - if (mini_info.minimizer < mini_info_rc.minimizer) { - return lookup_canonical(uint_kmer, uint_kmer_rc, mini_info); - } else if (mini_info_rc.minimizer < mini_info.minimizer) { - return lookup_canonical(uint_kmer, uint_kmer_rc, mini_info_rc); - } else { - auto res = lookup_canonical(uint_kmer, uint_kmer_rc, mini_info); - if (res.kmer_id == constants::invalid_uint64) { - res = lookup_canonical(uint_kmer, uint_kmer_rc, mini_info_rc); - } - return res; - } -} - -template -lookup_result dictionary::lookup_canonical(const Kmer uint_kmer, // - const Kmer uint_kmer_rc, // - const minimizer_info mini_info) const // -{ - assert(mini_info.minimizer == - std::min(util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher).minimizer, - util::compute_minimizer(uint_kmer_rc, m_k, m_m, m_hasher).minimizer)); - - const Kmer uint_kmer_canon = std::min(uint_kmer, uint_kmer_rc); - auto it = m_ssi.lookup(uint_kmer_canon, mini_info); - return m_spss.lookup_canonical(it, uint_kmer, uint_kmer_rc, mini_info); -} - template lookup_result dictionary::lookup(char const* string_kmer, bool check_reverse_complement) const { @@ -65,17 +14,28 @@ template lookup_result dictionary::lookup(Kmer uint_kmer, bool check_reverse_complement) const // { - if (m_canonical) return lookup_canonical(uint_kmer); - auto res = lookup_regular(uint_kmer); - assert(res.kmer_orientation == constants::forward_orientation); - if (check_reverse_complement and res.kmer_id == constants::invalid_uint64) { - Kmer uint_kmer_rc = uint_kmer; - uint_kmer_rc.reverse_complement_inplace(m_k); - res = lookup_regular(uint_kmer_rc); - res.kmer_orientation = constants::backward_orientation; + Kmer uint_kmer_rc = uint_kmer; + uint_kmer_rc.reverse_complement_inplace(m_k); + /* the reverse complement is needed anyway to tell the two orientations + apart, so computing the minimizer from it costs nothing extra */ + auto mini_info = util::compute_minimizer(uint_kmer, uint_kmer_rc, m_k, m_m, m_hasher); + auto res = lookup(uint_kmer, uint_kmer_rc, mini_info); + if (!check_reverse_complement and res.kmer_orientation == constants::backward_orientation) { // + return lookup_result(); } return res; } +template +lookup_result dictionary::lookup(const Kmer uint_kmer, // + const Kmer uint_kmer_rc, // + const minimizer_info mini_info, // + typename spss_type::bucket_cache* cache) const // +{ + assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == + util::compute_minimizer(uint_kmer, uint_kmer_rc, m_k, m_m, m_hasher)); + auto it = m_ssi.lookup(uint_kmer, uint_kmer_rc, mini_info); + return m_spss.lookup(it, uint_kmer, uint_kmer_rc, mini_info, cache); +} template bool dictionary::is_member(char const* string_kmer, @@ -203,7 +163,7 @@ neighbourhood dictionary::string_neighbours( template uint64_t dictionary::num_bits() const { return 8 * (sizeof(m_vnum) + sizeof(m_num_kmers) + sizeof(m_num_strings) + sizeof(m_k) + - sizeof(m_m) + sizeof(m_canonical) + sizeof(m_hasher)) + + sizeof(m_m) + sizeof(m_hasher)) + m_spss.num_bits() + m_ssi.num_bits() + m_weights.num_bits(); } diff --git a/src/info.cpp b/src/info.cpp index cf214fc..a1c5334 100644 --- a/src/info.cpp +++ b/src/info.cpp @@ -55,7 +55,6 @@ void dictionary::print_info() const { std::cout << "k = " << k() << '\n'; std::cout << "num_minimizers = " << m_ssi.codewords.size() << std::endl; std::cout << "m = " << m() << '\n'; - std::cout << "canonical = " << (canonical() ? "true" : "false") << '\n'; std::cout << "weighted = " << (weighted() ? "true" : "false") << '\n'; print_space_breakdown(); } diff --git a/src/query.cpp b/src/query.cpp index e367eee..3d96517 100644 --- a/src/query.cpp +++ b/src/query.cpp @@ -38,6 +38,8 @@ streaming_query_report streaming_query_from_fasta_file_multiline(Dict const* dic } report.num_searches = query.num_searches(); report.num_extensions = query.num_extensions(); + report.num_skipped_singleton_lookups = query.num_skipped_singleton_lookups(); + report.num_bucket_cache_hits = query.num_bucket_cache_hits(); report.num_positive_kmers = query.num_positive_lookups(); report.num_negative_kmers = query.num_negative_lookups(); report.num_invalid_kmers = query.num_invalid_lookups(); @@ -67,6 +69,8 @@ streaming_query_report streaming_query_from_fasta_file(Dict const* dict, std::is } report.num_searches = query.num_searches(); report.num_extensions = query.num_extensions(); + report.num_skipped_singleton_lookups = query.num_skipped_singleton_lookups(); + report.num_bucket_cache_hits = query.num_bucket_cache_hits(); report.num_positive_kmers = query.num_positive_lookups(); report.num_negative_kmers = query.num_negative_lookups(); report.num_invalid_kmers = query.num_invalid_lookups(); @@ -99,6 +103,8 @@ streaming_query_report streaming_query_from_fastq_file(Dict const* dict, std::is } report.num_searches = query.num_searches(); report.num_extensions = query.num_extensions(); + report.num_skipped_singleton_lookups = query.num_skipped_singleton_lookups(); + report.num_bucket_cache_hits = query.num_bucket_cache_hits(); report.num_positive_kmers = query.num_positive_lookups(); report.num_negative_kmers = query.num_negative_lookups(); report.num_invalid_kmers = query.num_invalid_lookups(); @@ -121,8 +127,7 @@ dictionary::streaming_query_from_file(std::string const& filename bool multiline) const // { using dictionary_type = dictionary; - using regular_query = streaming_query; - using canonical_query = streaming_query; + using query_type = streaming_query; std::ifstream is(filename.c_str()); if (!is.good()) throw std::runtime_error("error in opening the file '" + filename + "'"); @@ -130,42 +135,22 @@ dictionary::streaming_query_from_file(std::string const& filename if (util::ends_with(filename, ".fa.gz") or util::ends_with(filename, ".fasta.gz")) { zip_istream zis(is); - if (canonical()) { - report = streaming_query_from_fasta_file(this, zis, - multiline); - } else { - report = streaming_query_from_fasta_file(this, zis, - multiline); - } + report = streaming_query_from_fasta_file(this, zis, multiline); } else if (util::ends_with(filename, ".fq.gz") or util::ends_with(filename, ".fastq.gz")) { if (multiline) { std::cout << "==> Warning: option 'multiline' is only valid for FASTA files, not FASTQ." << std::endl; } zip_istream zis(is); - if (canonical()) { - report = streaming_query_from_fastq_file(this, zis); - } else { - report = streaming_query_from_fastq_file(this, zis); - } + report = streaming_query_from_fastq_file(this, zis); } else if (util::ends_with(filename, ".fa") or util::ends_with(filename, ".fasta")) { - if (canonical()) { - report = streaming_query_from_fasta_file(this, is, - multiline); - } else { - report = streaming_query_from_fasta_file(this, is, - multiline); - } + report = streaming_query_from_fasta_file(this, is, multiline); } else if (util::ends_with(filename, ".fq") or util::ends_with(filename, ".fastq")) { if (multiline) { std::cout << "==> Warning: option 'multiline' is only valid for FASTA files, not FASTQ." << std::endl; } - if (canonical()) { - report = streaming_query_from_fastq_file(this, is); - } else { - report = streaming_query_from_fastq_file(this, is); - } + report = streaming_query_from_fastq_file(this, is); } else { std::cerr << "unsupported query file format" << std::endl; } diff --git a/test/test_minimizer.cpp b/test/test_minimizer.cpp new file mode 100644 index 0000000..484bc0d --- /dev/null +++ b/test/test_minimizer.cpp @@ -0,0 +1,187 @@ +#include +#include +#include +#include + +#include "include/util.hpp" +#include "include/kmer.hpp" +#include "include/minimizer_iterator.hpp" + +using namespace sshash; + +static uint64_t num_failed = 0; +static uint64_t num_checked = 0; +static uint64_t num_ties = 0; + +static void fail(std::string const& what) { + if (num_failed < 10) std::cerr << "FAILED: " << what << std::endl; + ++num_failed; +} + +/* + The property the whole design rests on: the minimizer must be + mirror-equivariant. The m-mer selected for a kmer and for its reverse + complement must be literally the same string -- otherwise the two would land + in different buckets and a lookup by one orientation could not find the other + -- and the selected locus must be the mirror image. +*/ +template +static void check_equivariance(kmer_t kmer, uint64_t k, uint64_t m, hasher_type const& hasher) { + kmer_t kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(k); + + auto a = util::compute_minimizer(kmer, k, m, hasher); + auto b = util::compute_minimizer(kmer_rc, k, m, hasher); + + ++num_checked; + + if (a.minimizer != b.minimizer) { + fail("minimizer of kmer '" + util::uint_kmer_to_string(kmer, k) + "' (m=" + + std::to_string(m) + ") is " + util::uint_minimizer_to_string(a.minimizer, m) + + " but that of its reverse complement is " + + util::uint_minimizer_to_string(b.minimizer, m)); + return; + } + if (a.pos_in_kmer + b.pos_in_kmer != k - m) { + fail("locus of kmer '" + util::uint_kmer_to_string(kmer, k) + "' is " + + std::to_string(a.pos_in_kmer) + " but that of its reverse complement is " + + std::to_string(b.pos_in_kmer) + " (should mirror to " + + std::to_string(k - m - a.pos_in_kmer) + ")"); + return; + } + + /* the reported value must really be the canonical m-mer at the reported locus */ + kmer_t anchored = kmer; + anchored.drop_chars(a.pos_in_kmer); + anchored.take_chars(m); + if (a.minimizer != util::canonical_mmer(uint64_t(anchored), m)) { + fail("minimizer value does not match the m-mer at the selected locus"); + } +} + +/* Count the loci attaining the minimum hash, to report how often ties fire. */ +template +static bool has_tie(kmer_t kmer, uint64_t k, uint64_t m, hasher_type const& hasher) { + uint64_t min_hash = constants::invalid_uint64; + uint64_t count = 0; + for (uint64_t i = 0; i != k - m + 1; ++i) { + kmer_t mmer = kmer; + mmer.take_chars(m); + uint64_t hash = hasher.hash(util::canonical_mmer(uint64_t(mmer), m)); + if (i == 0 or hash < min_hash) { + min_hash = hash; + count = 1; + } else if (hash == min_hash) { + ++count; + } + kmer.drop_char(); + } + return count > 1; +} + +template +static void test_kmers(uint64_t k, uint64_t m, uint64_t num_kmers, std::mt19937_64& gen) { + hasher_type hasher(1234); + std::string s(k, 0); + for (uint64_t n = 0; n != num_kmers; ++n) { + for (uint64_t i = 0; i != k; ++i) s[i] = "ACGT"[gen() & 3]; + kmer_t kmer = util::string_to_uint_kmer(s.data(), k); + if (has_tie(kmer, k, m, hasher)) ++num_ties; + check_equivariance(kmer, k, m, hasher); + } +} + +/* + The incremental "re-scan" iterator must agree with the brute-force reference + on every window of a sequence, including the windows where a tie fires. + + The scheme must also be *forward*: the sampled position must never decrease + as the window slides, so that the number of super-kmers equals the number of + sampled positions. The centre-closest tie-break is what guarantees this + (leftmost-if-canonical / rightmost-otherwise, the previous rule, is + mirror-equivariant too but moves the anchor backwards on ~1e-5 of the + windows). +*/ +template +static void test_iterator(uint64_t k, uint64_t m, uint64_t length, std::mt19937_64& gen) { + hasher_type hasher(1234); + std::string s(length, 0); + for (uint64_t i = 0; i != length; ++i) s[i] = "ACGT"[gen() & 3]; + + minimizer_iterator it(k, m, hasher); + it.set_position(0); + + uint64_t prev_pos_in_seq = 0; + for (uint64_t i = 0; i + k <= length; ++i) { + kmer_t kmer = util::string_to_uint_kmer(s.data() + i, k); + auto got = it.next(kmer); + auto expected = util::compute_minimizer(kmer, k, m, hasher); + ++num_checked; + if (got.minimizer != expected.minimizer or got.pos_in_kmer != expected.pos_in_kmer) { + fail("iterator disagrees with the reference at position " + std::to_string(i) + + " (k=" + std::to_string(k) + ", m=" + std::to_string(m) + ")"); + return; + } + /* pos_in_seq must be the absolute position of the selected locus */ + if (got.pos_in_seq != i + got.pos_in_kmer) { + fail("iterator reports pos_in_seq " + std::to_string(got.pos_in_seq) + + " but expected " + std::to_string(i + got.pos_in_kmer)); + return; + } + /* forwardness: the sampled position never decreases */ + if (got.pos_in_seq < prev_pos_in_seq) { + fail("scheme is not forward: sampled position " + std::to_string(got.pos_in_seq) + + " after " + std::to_string(prev_pos_in_seq) + " (k=" + std::to_string(k) + + ", m=" + std::to_string(m) + ")"); + return; + } + prev_pos_in_seq = got.pos_in_seq; + } +} + +int main() { + using kmer_t = dna_uint_kmer_t; + using wide_kmer_t = dna_uint_kmer_t<__uint128_t>; + + std::mt19937_64 gen(42); + + std::cout << "checking mirror-equivariance of the minimizer..." << std::endl; + for (uint64_t k : {5, 15, 21, 31}) { + for (uint64_t m = 1; m <= std::min(k, kmer_t::max_m); ++m) { + test_kmers(k, m, 20000, gen); + } + } + for (uint64_t k : {31, 47, 63}) { + for (uint64_t m : {1, 2, 3, 7, 13, 21, 31}) { + if (m > k) continue; + test_kmers(k, m, 20000, gen); + } + } + + /* + Small m makes the m-mer universe tiny, so a window very often contains + two loci carrying the same canonical m-mer: this is what stresses the + tie-breaking rule, which on real data with m=13 fires on ~1e-5 of the + windows only. + */ + std::cout << " " << num_ties << "/" << num_checked << " of the kmers checked had a tie" + << std::endl; + + std::cout << "checking the incremental iterator against the reference..." << std::endl; + for (uint64_t k : {5, 15, 31}) { + for (uint64_t m = 1; m <= std::min(k, kmer_t::max_m); ++m) { + test_iterator(k, m, 20000, gen); + } + } + for (uint64_t k : {47, 63}) { + for (uint64_t m : {2, 3, 13, 31}) { test_iterator(k, m, 20000, gen); } + } + + std::cout << "checked " << num_checked << " kmers" << std::endl; + if (num_failed != 0) { + std::cerr << num_failed << " CHECKS FAILED" << std::endl; + return 1; + } + std::cout << "EVERYTHING OK!" << std::endl; + return 0; +} diff --git a/tools/build.cpp b/tools/build.cpp index 6630386..977c4dc 100644 --- a/tools/build.cpp +++ b/tools/build.cpp @@ -39,9 +39,6 @@ int build(int argc, char** argv) { "RAM limit in GiB. Default value is " + std::to_string(constants::default_ram_limit_in_GiB) + " GiB.", "-g", false); - parser.add("canonical", - "This option results in a trade-off between index space and lookup time.", - "--canonical", false, true); parser.add("weighted", "Also store the weights in compressed format.", "--weighted", false, true); parser.add("check", "Check correctness after construction.", "--check", false, true); @@ -59,7 +56,6 @@ int build(int argc, char** argv) { if (parser.parsed("seed")) build_config.seed = parser.get("seed"); if (parser.parsed("lambda")) build_config.lambda = parser.get("lambda"); - build_config.canonical = parser.get("canonical"); build_config.weighted = parser.get("weighted"); build_config.verbose = parser.get("verbose"); if (parser.parsed("tmp_dirname")) { diff --git a/tools/query.cpp b/tools/query.cpp index 75bdcdd..6793dcf 100644 --- a/tools/query.cpp +++ b/tools/query.cpp @@ -44,6 +44,8 @@ int query(int argc, char** argv) { query_stats.add("num_invalid_kmers", report.num_invalid_kmers); query_stats.add("num_searches", report.num_searches); query_stats.add("num_extensions", report.num_extensions); + query_stats.add("num_skipped_singleton_lookups", report.num_skipped_singleton_lookups); + query_stats.add("num_bucket_cache_hits", report.num_bucket_cache_hits); query_stats.add("elapsed_millisec", uint64_t(t.elapsed())); std::cout << "==== query report:\n"; @@ -60,6 +62,9 @@ int query(int argc, char** argv) { std::cout << "num_extensions = " << report.num_extensions << "/" << report.num_positive_kmers << " (" << (report.num_extensions * 100.0) / report.num_positive_kmers << "%)" << std::endl; + std::cout << "num_skipped_singleton_lookups = " << report.num_skipped_singleton_lookups + << std::endl; + std::cout << "num_bucket_cache_hits = " << report.num_bucket_cache_hits << std::endl; std::cout << "elapsed = " << t.elapsed() / 1000 << " sec / "; std::cout << t.elapsed() / 1000 / 60 << " min / "; std::cout << (t.elapsed() * 1e6) / report.num_kmers << " ns/kmer" << std::endl; diff --git a/tools/sshash.cpp b/tools/sshash.cpp index 8138a93..5d14509 100644 --- a/tools/sshash.cpp +++ b/tools/sshash.cpp @@ -53,7 +53,6 @@ int bench(int argc, char** argv) { perf_stats.add("index_filename", index_filename.c_str()); perf_stats.add("k", dict.k()); perf_stats.add("m", dict.m()); - perf_stats.add("canonical", dict.canonical() ? "true" : "false"); perf_test_lookup_access(dict, perf_stats); if (dict.weighted()) perf_test_lookup_weight(dict, perf_stats);