Minor bug fix (infinity loop) in benchmarks for 3 errors and k+2=5 parts

BenchmarkCode/benchmark.cpp

@@ -148,6 +148,168 @@ void BM_EditDistance(benchmark::State& state, TSearchScheme scheme)
     }
 }
 
+template <typename TIter, typename TText, typename TOptimalScheme, typename TDistanceTag>
+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 <typename TIter, typename TText, typename TOptimalScheme, typename TDistanceTag>
+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);

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<OptimalSearch<5>, 4> VALUE_plus_two
+    static constexpr std::array<OptimalSearch<5>, 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 <typename TVoidType>
 constexpr std::array<OptimalSearch<4>, 3> PaperOptimumSearchSchemes<3, TVoidType>::VALUE_plus_one;
 
 template <typename TVoidType>
-constexpr std::array<OptimalSearch<5>, 4> PaperOptimumSearchSchemes<3, TVoidType>::VALUE_plus_two;
+constexpr std::array<OptimalSearch<5>, 3> PaperOptimumSearchSchemes<3, TVoidType>::VALUE_plus_two;
 
 template <typename TVoidType>
 constexpr std::array<OptimalSearch<6>, 3> PaperOptimumSearchSchemes<3, TVoidType>::VALUE_plus_three;