diff --git a/BenchmarkCode/benchmark.cpp b/BenchmarkCode/benchmark.cpp index 107ec2d..90fe249 100644 --- a/BenchmarkCode/benchmark.cpp +++ b/BenchmarkCode/benchmark.cpp @@ -148,6 +148,168 @@ void BM_EditDistance(benchmark::State& state, TSearchScheme scheme) } } +template +inline bool runStrataBenchmarkStep(TIter & it, TText & read, TOptimalScheme const & scheme, TDistanceTag const & /**/) +{ + unsigned hitsNbr = 0; + auto delegate = [&hitsNbr](TIter const & it, TText const & /*read*/, unsigned const /*errors*/) { + ++hitsNbr; + unsigned x = 0; + for (unsigned i = 0; i < length(getOccurrences(it)); ++i) + x += getOccurrences(it)[i].i2; + }; + + _optimalSearchScheme(delegate, it, read, scheme, TDistanceTag()); + reverseComplement(read); + _optimalSearchScheme(delegate, it, read, scheme, TDistanceTag()); + + return hitsNbr != 0; +} + +void BM_EditDistanceStrata0(benchmark::State& state, uint8_t const maxErrors) +{ + typedef EditDistance TDistanceTag; + + TIter it(fm_index); + + auto s0 = OptimalSearchSchemes<0, 0>::VALUE; + auto s1 = OptimalSearchSchemes<1, 1>::VALUE; + auto s2 = OptimalSearchSchemes<2, 2>::VALUE; + auto s3 = OptimalSearchSchemes<3, 3>::VALUE; + auto s4 = OptimalSearchSchemes<4, 4>::VALUE; + _optimalSearchSchemeComputeFixedBlocklength(s0, length(reads[0])); + _optimalSearchSchemeComputeFixedBlocklength(s1, length(reads[0])); + _optimalSearchSchemeComputeFixedBlocklength(s2, length(reads[0])); + _optimalSearchSchemeComputeFixedBlocklength(s3, length(reads[0])); + _optimalSearchSchemeComputeFixedBlocklength(s4, length(reads[0])); + + unsigned hits[5] = {}; + for (auto _ : state) + { + for (unsigned i = 0; i < length(reads); ++i) + { + if (runStrataBenchmarkStep(it, reads[i], s0, TDistanceTag())) + ++hits[0]; + else if (maxErrors >= 1 && runStrataBenchmarkStep(it, reads[i], s1, TDistanceTag())) + ++hits[1]; + else if (maxErrors >= 2 && runStrataBenchmarkStep(it, reads[i], s2, TDistanceTag())) + ++hits[2]; + else if (maxErrors >= 3 && runStrataBenchmarkStep(it, reads[i], s3, TDistanceTag())) + ++hits[3]; + else if (maxErrors >= 4 && runStrataBenchmarkStep(it, reads[i], s4, TDistanceTag())) + ++hits[4]; + } + } + + for (unsigned e = 0; e <= maxErrors; ++e) + { + if (hits[e] > 0) + std::cout << "hits[" << e << "] = " << hits[e] << std::endl; + } +} + +template +inline unsigned runStrataBenchmarkStep1(TIter & it, TText & read, TOptimalScheme const & scheme, TDistanceTag const & /**/) +{ + unsigned bestErrors = 10; + auto delegate = [&bestErrors](TIter const & it, TText const & /*read*/, unsigned const errors) { + bestErrors = std::min(errors, bestErrors); + unsigned x = 0; + for (unsigned i = 0; i < length(getOccurrences(it)); ++i) + x += getOccurrences(it)[i].i2; + }; + + _optimalSearchScheme(delegate, it, read, scheme, TDistanceTag()); + reverseComplement(read); + _optimalSearchScheme(delegate, it, read, scheme, TDistanceTag()); + + return bestErrors; +} + +void BM_EditDistanceStrata1(benchmark::State& state, uint8_t const maxErrors) +{ + typedef EditDistance TDistanceTag; + + TIter it(fm_index); + + auto s0_0 = OptimalSearchSchemes<0, 0>::VALUE; + auto s0_1 = OptimalSearchSchemes<0, 1>::VALUE; + auto s0_2 = OptimalSearchSchemes<0, 2>::VALUE; + auto s2_2 = OptimalSearchSchemes<2, 2>::VALUE; + auto s2_3 = OptimalSearchSchemes<2, 3>::VALUE; + auto s4_4 = OptimalSearchSchemes<4, 4>::VALUE; + _optimalSearchSchemeComputeFixedBlocklength(s0_0, length(reads[0])); + _optimalSearchSchemeComputeFixedBlocklength(s0_1, length(reads[0])); + _optimalSearchSchemeComputeFixedBlocklength(s0_2, length(reads[0])); + _optimalSearchSchemeComputeFixedBlocklength(s2_2, length(reads[0])); + _optimalSearchSchemeComputeFixedBlocklength(s2_3, length(reads[0])); + _optimalSearchSchemeComputeFixedBlocklength(s4_4, length(reads[0])); + + unsigned tmp = 0; + + unsigned hits[5] = {}; + for (auto _ : state) + { + for (unsigned i = 0; i < length(reads); ++i) + { + unsigned best = runStrataBenchmarkStep1(it, reads[i], s0_1, TDistanceTag()); + if (best == 0) + ++hits[0]; + else if (best == 1) + { + ++hits[1]; + if (maxErrors > 1) + tmp += runStrataBenchmarkStep1(it, reads[i], s2_2, TDistanceTag()); + } + else if (maxErrors > 1) + { + if (maxErrors == 2) + { + best = runStrataBenchmarkStep1(it, reads[i], s2_2, TDistanceTag()); + if (best == 2) + ++hits[2]; + } + else + { + best = runStrataBenchmarkStep1(it, reads[i], s2_3, TDistanceTag()); + if (best == 2) + ++hits[2]; + else if (best == 3) + { + ++hits[3]; + if (maxErrors > 3) + tmp += runStrataBenchmarkStep1(it, reads[i], s4_4, TDistanceTag()); + } + else if (maxErrors == 4) + { + best = runStrataBenchmarkStep1(it, reads[i], s4_4, TDistanceTag()); + if (best == 4) + ++hits[4]; + } + } + } + } + } + std::cout << "tmp: " << tmp << '\n'; + + for (unsigned e = 0; e <= maxErrors; ++e) + { + if (hits[e] > 0) + std::cout << "hits[" << e << "] = " << hits[e] << std::endl; + } +} + +BENCHMARK_CAPTURE(BM_EditDistanceStrata1, errors_1_strata_1, (uint8_t)1)->Unit(benchmark::kMillisecond); +BENCHMARK_CAPTURE(BM_EditDistanceStrata1, errors_2_strata_1, (uint8_t)2)->Unit(benchmark::kMillisecond); +BENCHMARK_CAPTURE(BM_EditDistanceStrata1, errors_3_strata_1, (uint8_t)3)->Unit(benchmark::kMillisecond); +BENCHMARK_CAPTURE(BM_EditDistanceStrata1, errors_4_strata_1, (uint8_t)4)->Unit(benchmark::kMillisecond); + +BENCHMARK_CAPTURE(BM_EditDistanceStrata0, errors_0_strata_0, (uint8_t)0)->Unit(benchmark::kMillisecond); +BENCHMARK_CAPTURE(BM_EditDistanceStrata0, errors_1_strata_0, (uint8_t)1)->Unit(benchmark::kMillisecond); +BENCHMARK_CAPTURE(BM_EditDistanceStrata0, errors_2_strata_0, (uint8_t)2)->Unit(benchmark::kMillisecond); +BENCHMARK_CAPTURE(BM_EditDistanceStrata0, errors_3_strata_0, (uint8_t)3)->Unit(benchmark::kMillisecond); +BENCHMARK_CAPTURE(BM_EditDistanceStrata0, errors_4_strata_0, (uint8_t)4)->Unit(benchmark::kMillisecond); + BENCHMARK_CAPTURE(BM_HammingDistance, errors_1_backtracking , (uint8_t)1)->Unit(benchmark::kMillisecond); BENCHMARK_CAPTURE(BM_HammingDistance, errors_1_parts_k_plus_1, PaperOptimumSearchSchemes<1>::VALUE_plus_one)->Unit(benchmark::kMillisecond); BENCHMARK_CAPTURE(BM_HammingDistance, errors_1_parts_k_plus_2, PaperOptimumSearchSchemes<1>::VALUE_plus_two)->Unit(benchmark::kMillisecond); diff --git a/BenchmarkCode/paper_optimum_schemes.h b/BenchmarkCode/paper_optimum_schemes.h index 6f85395..829d80e 100644 --- a/BenchmarkCode/paper_optimum_schemes.h +++ b/BenchmarkCode/paper_optimum_schemes.h @@ -77,7 +77,7 @@ struct PaperOptimumSearchSchemes<3, TVoidType> { {{2, 3, 4, 1}}, {{0, 0, 0, 0}}, {{1, 2, 2, 3}}, {{0, 0, 0, 0}}, 0}, { {{3, 4, 2, 1}}, {{0, 0, 2, 2}}, {{0, 0, 3, 3}}, {{0, 0, 0, 0}}, 0} }}; - static constexpr std::array, 4> VALUE_plus_two + static constexpr std::array, 3> VALUE_plus_two {{ { {{1, 2, 3, 4, 5}}, {{0, 0, 0, 2, 2}}, {{0, 0, 3, 3, 3}}, {{0, 0, 0, 0, 0}}, 0}, { {{4, 3, 2, 1, 5}}, {{0, 0, 0, 0, 0}}, {{1, 1, 2, 2, 3}}, {{0, 0, 0, 0, 0}}, 0}, @@ -99,7 +99,7 @@ template constexpr std::array, 3> PaperOptimumSearchSchemes<3, TVoidType>::VALUE_plus_one; template -constexpr std::array, 4> PaperOptimumSearchSchemes<3, TVoidType>::VALUE_plus_two; +constexpr std::array, 3> PaperOptimumSearchSchemes<3, TVoidType>::VALUE_plus_two; template constexpr std::array, 3> PaperOptimumSearchSchemes<3, TVoidType>::VALUE_plus_three;