Merge pull request #402 from PowerGridModel/feature/refactor-replace-…

…indptrs

Replace indptr with sparse and dense grouped index vector

CMakePresets.json

@@ -89,7 +89,7 @@
       "environment": {
         "CC": "clang-cl.exe",
         "CXX": "clang-cl.exe",
-        "CXXFLAGS": "$penv{CXXFLAGS} -Wno-uninitialized-const-reference"
+        "CXXFLAGS": "$penv{CXXFLAGS} -Wno-uninitialized-const-reference -Wno-unknown-attributes"
       },
       "inherits": "windows-base"
     },
@@ -115,7 +115,7 @@
         "CMAKE_C_CLANG_TIDY": "clang-tidy.exe;--extra-arg=/EHsc",
         "CMAKE_CXX_CLANG_TIDY": "clang-tidy.exe;--extra-arg=/EHsc"
       },
-      "inherits": "windows-base",
+      "inherits": "clang-cl-base",
       "hidden": true
     },
     {

README.md

@@ -3,8 +3,8 @@ SPDX-FileCopyrightText: 2022 Contributors to the Power Grid Model project <dynam
 
 SPDX-License-Identifier: MPL-2.0
 -->
-[![PyPI version](https://badge.fury.io/py/power-grid-model.svg)](https://badge.fury.io/py/power-grid-model)
-[![Anaconda-Server Badge](https://anaconda.org/conda-forge/power-grid-model/badges/version.svg)](https://anaconda.org/conda-forge/power-grid-model)
+[![PyPI version](https://badge.fury.io/py/power-grid-model.svg?no-cache)](https://badge.fury.io/py/power-grid-model)
+[![Anaconda-Server Badge](https://anaconda.org/conda-forge/power-grid-model/badges/version.svg?no-cache)](https://anaconda.org/conda-forge/power-grid-model)
 [![License: MIT](https://img.shields.io/badge/License-MPL2.0-informational.svg)](https://github.com/PowerGridModel/power-grid-model/blob/main/LICENSE)
 [![Build and Test C++ and Python](https://github.com/PowerGridModel/power-grid-model/actions/workflows/main.yml/badge.svg)](https://github.com/PowerGridModel/power-grid-model/actions/workflows/main.yml)
 [![Check Code Quality](https://github.com/PowerGridModel/power-grid-model/actions/workflows/check-code-quality.yml/badge.svg)](https://github.com/PowerGridModel/power-grid-model/actions/workflows/check-code-quality.yml)

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -7,6 +7,7 @@
 #define POWER_GRID_MODEL_CALCULATION_PARAMETERS_HPP
 
 #include "enum.hpp"
+#include "grouped_index_vector.hpp"
 #include "power_grid_model.hpp"
 #include "three_phase_tensor.hpp"
 
@@ -92,41 +93,41 @@ struct MathModelTopology {
     std::vector<double> phase_shift;
     std::vector<BranchIdx> branch_bus_idx;
     std::vector<BranchIdx> fill_in;
-    IdxVector source_bus_indptr;
-    IdxVector shunt_bus_indptr;
-    IdxVector load_gen_bus_indptr;
+    DenseGroupedIdxVector sources_per_bus;
+    SparseGroupedIdxVector shunts_per_bus;
+    SparseGroupedIdxVector load_gens_per_bus;
     std::vector<LoadGenType> load_gen_type;
-    IdxVector voltage_sensor_indptr;
-    IdxVector source_power_sensor_indptr;      // indptr of the source
-    IdxVector load_gen_power_sensor_indptr;    // indptr of the load_gen
-    IdxVector shunt_power_sensor_indptr;       // indptr of the shunt
-    IdxVector branch_from_power_sensor_indptr; // indptr of the branch
-    IdxVector branch_to_power_sensor_indptr;   // indptr of the branch
-    IdxVector bus_power_sensor_indptr;         // indptr of the bus
+    SparseGroupedIdxVector voltage_sensors_per_bus;
+    SparseGroupedIdxVector power_sensors_per_source;
+    SparseGroupedIdxVector power_sensors_per_load_gen;
+    SparseGroupedIdxVector power_sensors_per_shunt;
+    SparseGroupedIdxVector power_sensors_per_branch_from;
+    SparseGroupedIdxVector power_sensors_per_branch_to;
+    SparseGroupedIdxVector power_sensors_per_bus;
 
     Idx n_bus() const { return static_cast<Idx>(phase_shift.size()); }
 
     Idx n_branch() const { return static_cast<Idx>(branch_bus_idx.size()); }
 
-    Idx n_source() const { return source_bus_indptr.back(); }
+    Idx n_source() const { return sources_per_bus.element_size(); }
 
-    Idx n_shunt() const { return shunt_bus_indptr.back(); }
+    Idx n_shunt() const { return shunts_per_bus.element_size(); }
 
-    Idx n_load_gen() const { return load_gen_bus_indptr.back(); }
+    Idx n_load_gen() const { return load_gens_per_bus.element_size(); }
 
-    Idx n_voltage_sensor() const { return voltage_sensor_indptr.back(); }
+    Idx n_voltage_sensor() const { return voltage_sensors_per_bus.element_size(); }
 
-    Idx n_source_power_sensor() const { return source_power_sensor_indptr.back(); }
+    Idx n_source_power_sensor() const { return power_sensors_per_source.element_size(); }
 
-    Idx n_load_gen_power_sensor() const { return load_gen_power_sensor_indptr.back(); }
+    Idx n_load_gen_power_sensor() const { return power_sensors_per_load_gen.element_size(); }
 
-    Idx n_shunt_power_power_sensor() const { return shunt_power_sensor_indptr.back(); }
+    Idx n_shunt_power_power_sensor() const { return power_sensors_per_shunt.element_size(); }
 
-    Idx n_branch_from_power_sensor() const { return branch_from_power_sensor_indptr.back(); }
+    Idx n_branch_from_power_sensor() const { return power_sensors_per_branch_from.element_size(); }
 
-    Idx n_branch_to_power_sensor() const { return branch_to_power_sensor_indptr.back(); }
+    Idx n_branch_to_power_sensor() const { return power_sensors_per_branch_to.element_size(); }
 
-    Idx n_bus_power_sensor() const { return bus_power_sensor_indptr.back(); }
+    Idx n_bus_power_sensor() const { return power_sensors_per_bus.element_size(); }
 };
 
 template <bool sym> struct MathModelParam {

power_grid_model_c/power_grid_model/include/power_grid_model/grouped_index_vector.hpp

@@ -7,6 +7,7 @@
 #define POWER_GRID_MODEL_SPARSE_IDX_VECTOR_HPP
 
 #include "power_grid_model.hpp"
+#include "typing.hpp"
 
 #include <boost/iterator/iterator_facade.hpp>
 #include <boost/iterator/zip_iterator.hpp>
@@ -51,14 +52,6 @@ inline auto sparse_decode(IdxVector const& indptr) {
     return result;
 }
 
-template <typename T, typename U> constexpr auto narrow_cast(U value) {
-    if constexpr (std::same_as<T, U>) {
-        return value;
-    }
-    assert(std::in_range<T>(value));
-    return static_cast<T>(value);
-}
-
 // TODO(mgovers): replace the below relevant subset here ourselves with the STD equivalent when we have std::ranges.
 // boost::counting_iterator does not satisfy all requirements std::*_iterator concepts:
 static_assert(!std::random_access_iterator<IdxCount>);
@@ -97,6 +90,8 @@ concept index_range_iterator =
 
 template <typename T>
 concept grouped_index_vector_type = std::default_initializable<T> && requires(T const t, Idx const idx) {
+                                                                         typename T::iterator;
+
                                                                          { t.size() } -> std::same_as<Idx>;
 
                                                                          { t.begin() } -> index_range_iterator;
@@ -152,6 +147,8 @@ class SparseGroupedIdxVector {
     auto group_iterator(Idx group) const { return GroupIterator{indptr_, group}; }
 
   public:
+    using iterator = GroupIterator;
+
     auto size() const { return static_cast<Idx>(indptr_.size()) - 1; }
     auto begin() const { return group_iterator(0); }
     auto end() const { return group_iterator(size()); }
@@ -205,8 +202,8 @@ class DenseGroupedIdxVector {
         auto dereference() const -> iterator {
             assert(dense_vector_ != nullptr);
             return boost::counting_range(
-                detail::narrow_cast<Idx>(std::distance(std::cbegin(*dense_vector_), group_range_.first)),
-                detail::narrow_cast<Idx>(std::distance(std::cbegin(*dense_vector_), group_range_.second)));
+                narrow_cast<Idx>(std::distance(std::cbegin(*dense_vector_), group_range_.first)),
+                narrow_cast<Idx>(std::distance(std::cbegin(*dense_vector_), group_range_.second)));
         }
         constexpr auto equal(GroupIterator const& other) const { return group_ == other.group_; }
         constexpr auto distance_to(GroupIterator const& other) const { return other.group_ - group_; }
@@ -225,6 +222,8 @@ class DenseGroupedIdxVector {
     auto group_iterator(Idx group) const { return GroupIterator{dense_vector_, group}; }
 
   public:
+    using iterator = GroupIterator;
+
     auto size() const { return num_groups_; }
     auto begin() const { return group_iterator(Idx{}); }
     auto end() const { return group_iterator(size()); }

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/common_solver_functions.hpp

@@ -13,10 +13,9 @@
 namespace power_grid_model::common_solver_functions {
 
 template <bool sym>
-void add_sources(IdxVector const& source_bus_indptr, Idx const& bus_number, YBus<sym> const& y_bus,
+void add_sources(grouped_idx_vector_type auto const& sources_per_bus, Idx const& bus_number, YBus<sym> const& y_bus,
                  ComplexVector const& u_source_vector, ComplexTensor<sym>& diagonal_element, ComplexValue<sym>& u_bus) {
-    for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
-         ++source_number) {
+    for (Idx const source_number : sources_per_bus.get_element_range(bus_number)) {
         ComplexTensor<sym> const y_source = y_bus.math_model_param().source_param[source_number];
         diagonal_element += y_source; // add y_source to the diagonal of Ybus
         u_bus += dot(y_source, ComplexValue<sym>{u_source_vector[source_number]}); // rhs += Y_source * U_source
@@ -35,8 +34,8 @@ template <bool sym> void copy_y_bus(YBus<sym> const& y_bus, ComplexTensorVector<
 
 template <bool sym>
 void calculate_source_result(Idx const& bus_number, YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
-                             MathOutput<sym>& output, IdxVector const& source_bus_indptr) {
-    for (Idx source = source_bus_indptr[bus_number]; source != source_bus_indptr[bus_number + 1]; ++source) {
+                             MathOutput<sym>& output, grouped_idx_vector_type auto const& sources_per_bus) {
+    for (Idx const source : sources_per_bus.get_element_range(bus_number)) {
         ComplexValue<sym> const u_ref{input.source[source]};
         ComplexTensor<sym> const y_ref = y_bus.math_model_param().source_param[source];
         output.source[source].i = dot(y_ref, u_ref - output.u[bus_number]);
@@ -46,8 +45,9 @@ void calculate_source_result(Idx const& bus_number, YBus<sym> const& y_bus, Powe
 
 template <bool sym, class LoadGenFunc>
 void calculate_load_gen_result(Idx const& bus_number, PowerFlowInput<sym> const& input, MathOutput<sym>& output,
-                               IdxVector const& load_gen_bus_indptr, LoadGenFunc const& load_gen_func) {
-    for (Idx load_gen = load_gen_bus_indptr[bus_number]; load_gen != load_gen_bus_indptr[bus_number + 1]; ++load_gen) {
+                               grouped_idx_vector_type auto const& load_gens_per_bus,
+                               LoadGenFunc const& load_gen_func) {
+    for (auto load_gen : load_gens_per_bus.get_element_range(bus_number)) {
         switch (LoadGenType const type = load_gen_func(load_gen); type) {
             using enum LoadGenType;
 

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp

@@ -85,7 +85,7 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
 
     // Add source admittance to Y bus and set variable for prepared y bus to true
     void initialize_derived_solver(YBus<sym> const& y_bus, MathOutput<sym> const& /* output */) {
-        IdxVector const& source_bus_indptr = *this->source_bus_indptr_;
+        auto const& sources_per_bus = *this->sources_per_bus_;
         IdxVector const& bus_entry = y_bus.lu_diag();
         // if Y bus is not up to date
         // re-build matrix and prefactorize Build y bus data with source admittance
@@ -97,8 +97,7 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
             for (Idx bus_number = 0; bus_number != this->n_bus_; ++bus_number) {
                 Idx const data_sequence = bus_entry[bus_number];
                 // loop sources
-                for (Idx source_number = source_bus_indptr[bus_number];
-                     source_number != source_bus_indptr[bus_number + 1]; ++source_number) {
+                for (auto source_number : sources_per_bus.get_element_range(bus_number)) {
                     // YBus_diag += Y_source
                     mat_data[data_sequence] += y_bus.math_model_param().source_param[source_number];
                 }
@@ -117,17 +116,17 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
     // Prepare matrix calculates injected current ie. RHS of solver for each iteration.
     void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
                                 ComplexValueVector<sym> const& u) {
-        IdxVector const& load_gen_bus_indptr = *this->load_gen_bus_indptr_;
-        IdxVector const& source_bus_indptr = *this->source_bus_indptr_;
+        auto const& load_gens_per_bus = *this->load_gens_per_bus_;
+        auto const& sources_per_bus = *this->sources_per_bus_;
         std::vector<LoadGenType> const& load_gen_type = *this->load_gen_type_;
 
         // set rhs to zero for iteration start
         std::fill(rhs_u_.begin(), rhs_u_.end(), ComplexValue<sym>{0.0});
 
         // loop buses: i
         for (Idx bus_number = 0; bus_number != this->n_bus_; ++bus_number) {
-            add_loads(bus_number, input, load_gen_bus_indptr, load_gen_type, u);
-            add_sources(bus_number, y_bus, input, source_bus_indptr);
+            add_loads(bus_number, input, load_gens_per_bus, load_gen_type, u);
+            add_sources(bus_number, y_bus, input, sources_per_bus);
         }
     }
 
@@ -158,10 +157,10 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
     SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>> sparse_solver_;
     std::shared_ptr<BlockPermArray const> perm_;
 
-    void add_loads(Idx const& bus_number, PowerFlowInput<sym> const& input, IdxVector const& load_gen_bus_indptr,
-                   std::vector<LoadGenType> const& load_gen_type, ComplexValueVector<sym> const& u) {
-        for (Idx load_number = load_gen_bus_indptr[bus_number]; load_number != load_gen_bus_indptr[bus_number + 1];
-             ++load_number) {
+    void add_loads(Idx const& bus_number, PowerFlowInput<sym> const& input,
+                   grouped_idx_vector_type auto const& load_gens_per_bus, std::vector<LoadGenType> const& load_gen_type,
+                   ComplexValueVector<sym> const& u) {
+        for (auto load_number : load_gens_per_bus.get_element_range(bus_number)) {
             // load type
             LoadGenType const type = load_gen_type[load_number];
             switch (type) {
@@ -186,9 +185,8 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
     }
 
     void add_sources(Idx const& bus_number, YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
-                     IdxVector const& source_bus_indptr) {
-        for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
-             ++source_number) {
+                     grouped_idx_vector_type auto const& sources_per_bus) {
+        for (Idx const source_number : sources_per_bus.get_element_range(bus_number)) {
             // I_inj_i += Y_source_j * U_ref_j
             rhs_u_[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
                                       ComplexValue<sym>{input.source[source_number]});