update topology to not necessarily go through sparse mapping

power_grid_model_c/power_grid_model/include/power_grid_model/grouped_index_vector.hpp

@@ -15,6 +15,8 @@
 #include <boost/range/adaptor/indexed.hpp>
 #include <boost/range/counting_range.hpp>
 
+#include <ranges>
+
 /*
 A data-structure for iterating through the indptr, ie. sparse representation of data.
 Indptr can be eg: [0, 3, 6, 7]
@@ -39,7 +41,7 @@ namespace detail {
 inline auto sparse_encode(IdxVector const& element_groups, Idx num_groups) {
     IdxVector result(num_groups + 1);
     auto next_group = std::begin(element_groups);
-    for (Idx group = 0; group < num_groups; group++) {
+    for (auto const group : boost::counting_range(Idx{0}, num_groups)) {
         next_group = std::upper_bound(next_group, std::end(element_groups), group);
         result[group + 1] = std::distance(std::begin(element_groups), next_group);
     }
@@ -48,7 +50,7 @@ inline auto sparse_encode(IdxVector const& element_groups, Idx num_groups) {
 
 inline auto sparse_decode(IdxVector const& indptr) {
     auto result = IdxVector(indptr.back());
-    for (Idx const group : boost::counting_range(Idx{}, static_cast<Idx>(indptr.size()) - 1)) {
+    for (Idx const group : boost::counting_range(Idx{0}, static_cast<Idx>(indptr.size()) - 1)) {
         std::fill(std::begin(result) + indptr[group], std::begin(result) + indptr[group + 1], group);
     }
     return result;

power_grid_model_c/power_grid_model/include/power_grid_model/index_mapping.hpp

@@ -0,0 +1,142 @@
+// SPDX-FileCopyrightText: 2022 Contributors to the Power Grid Model project <dynamic.grid.calculation@alliander.com>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+#pragma once
+#ifndef POWER_GRID_MODEL_INDEX_MAPPING_HPP
+#define POWER_GRID_MODEL_INDEX_MAPPING_HPP
+
+#include "grouped_index_vector.hpp"
+#include "power_grid_model.hpp"
+
+#include <ranges>
+
+namespace power_grid_model {
+
+struct SparseIndexMapping {
+    IdxVector indptr;
+    IdxVector reorder;
+};
+
+/// Sparse mapping technique
+///
+/// using counting sort: https://en.wikipedia.org/wiki/Counting_sort
+/// Given a input idx array: idx_B_in_A[...] to couple an array of A and an array of B,
+/// i.e. idx_B_in_A[i_A] = j_B.
+///
+/// For entry i_A in the array, idx_B_in_A[i_A] is the idx of B which couples the A object #i_A.
+///
+/// This sparse mapping tries to build a CSC sparse matrix so that
+/// only entry (i, idx_B_in_A[i_A] = j_B) is filled with the sequence of A, i.e. #i_A.
+///
+/// In CSC format, the entries are ordered by the idx of B
+/// therefore, in indptr, each range indptr[j_B:j_B+1]
+/// represents the entries of A where B object #j_B is coupled.
+///
+/// The indices array is not interesting here.
+/// The data array is original index of A. We call it re-order array.
+/// This can be used to reorder A objects by the coupling of B idx.
+///
+/// Example.
+/// For original idx_B_in_A == [3, 5, 2, 1, 1, 2]
+/// size of A is 6
+/// size of B is 7
+/// Result matrix
+/// indptr = [0, 0, 2, 4, 5, 5, 6, 6]
+/// data/reorder = [3, 4, 2, 5, 0, 1]
+/// to read:
+///     nothing		coupled to B 0
+///     A 3, 4		coupled to B 1
+///     A 2, 5		coupled to B 2
+///     A 0			coupled to B 3
+///     nothing		coupled to B 4
+///     A 1			coupled	to B 5
+///     nothing	coupled to B 6
+inline SparseIndexMapping build_sparse_mapping(IdxVector const& idx_B_in_A, Idx const n_B) {
+    using SparseEntry = std::pair<Idx, Idx>;
+
+    auto const n_A = static_cast<Idx>(idx_B_in_A.size());
+
+    std::vector<SparseEntry> entries(n_A);
+    std::ranges::transform(idx_B_in_A, boost::counting_range(Idx{0}, static_cast<Idx>(idx_B_in_A.size())),
+                           entries.begin(), [](Idx j_B, Idx i_A) {
+                               return SparseEntry{i_A, j_B};
+                           });
+
+    SparseIndexMapping sparse_mapping;
+    sparse_mapping.indptr.resize(n_B + 1);
+    sparse_mapping.reorder.resize(n_A);
+
+    IdxVector counter(n_B, 0);
+    for (auto const& entry : entries) {
+        ++counter[entry.second];
+    }
+
+    std::inclusive_scan(counter.cbegin(), counter.cend(), sparse_mapping.indptr.begin() + 1);
+
+    std::copy(sparse_mapping.indptr.cbegin() + 1, sparse_mapping.indptr.cend(), counter.begin());
+
+    for (auto it_entry = entries.crbegin(); it_entry != entries.crend(); ++it_entry) {
+        sparse_mapping.reorder[--counter[it_entry->second]] = it_entry->first;
+    }
+
+    assert(sparse_mapping.indptr[0] == 0);
+    assert(sparse_mapping.indptr.back() == n_A);
+    return sparse_mapping;
+}
+
+struct DenseIndexMapping {
+    IdxVector indvector;
+    IdxVector reorder;
+};
+
+namespace detail {
+
+inline auto build_dense_mapping_comparison_sort(IdxVector const& idx_B_in_A, Idx const /* n_B */) {
+    using DenseEntry = std::pair<Idx, Idx>;
+
+    std::vector<DenseEntry> mapping_to_from;
+    mapping_to_from.reserve(idx_B_in_A.size());
+    std::ranges::transform(idx_B_in_A, boost::counting_range(Idx{0}, static_cast<Idx>(idx_B_in_A.size())),
+                           std::back_inserter(mapping_to_from), [](Idx value, Idx orig_idx) {
+                               return std::pair{value, orig_idx};
+                           });
+
+    std::ranges::sort(mapping_to_from);
+
+    DenseIndexMapping result;
+    result.indvector.reserve(mapping_to_from.size());
+    result.reorder.reserve(mapping_to_from.size());
+    std::ranges::transform(mapping_to_from, std::back_inserter(result.indvector),
+                           [](DenseEntry const& to_from) { return to_from.first; });
+
+    std::ranges::transform(mapping_to_from, std::back_inserter(result.reorder),
+                           [](DenseEntry const& to_from) { return to_from.second; });
+
+    return result;
+}
+
+inline auto build_dense_mapping_counting_sort(IdxVector const& idx_B_in_A, Idx const n_B) {
+    auto sparse_result = build_sparse_mapping(idx_B_in_A, n_B);
+
+    return DenseIndexMapping{.indvector = sparse_decode(sparse_result.indptr),
+                             .reorder = std::move(sparse_result.reorder)};
+}
+
+} // namespace detail
+
+inline DenseIndexMapping build_dense_mapping(IdxVector const& idx_B_in_A, Idx const n_B) {
+    constexpr auto relative_complexity_prefactor = 1.0;
+
+    auto const n_A_ = static_cast<double>(idx_B_in_A.size());
+    auto const n_B_ = static_cast<double>(n_B);
+
+    if (n_A_ + n_B_ < relative_complexity_prefactor * n_A_ * log(n_A_)) {
+        return detail::build_dense_mapping_counting_sort(idx_B_in_A, n_B);
+    }
+    return detail::build_dense_mapping_comparison_sort(idx_B_in_A, n_B);
+}
+
+} // namespace power_grid_model
+
+#endif

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/y_bus.hpp

@@ -8,7 +8,6 @@
 
 #include "../calculation_parameters.hpp"
 #include "../power_grid_model.hpp"
-#include "../sparse_mapping.hpp"
 #include "../three_phase_tensor.hpp"
 
 namespace power_grid_model {