apply indptr replace source bus (still segfaults)

Signed-off-by: Martijn Govers <Martijn.Govers@Alliander.com>

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,7 +93,7 @@ struct MathModelTopology {
     std::vector<double> phase_shift;
     std::vector<BranchIdx> branch_bus_idx;
     std::vector<BranchIdx> fill_in;
-    IdxVector source_bus_indptr;
+    SparseIdxVector source_buses; // TODO(mgovers) replace with DenseIdxVector
     IdxVector shunt_bus_indptr;
     IdxVector load_gen_bus_indptr;
     std::vector<LoadGenType> load_gen_type;
@@ -108,7 +109,7 @@ struct MathModelTopology {
 
     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 source_buses.size(); }
 
     Idx n_shunt() const { return shunt_bus_indptr.back(); }
 

power_grid_model_c/power_grid_model/include/power_grid_model/grouped_index_vector.hpp

@@ -58,13 +58,20 @@ class SparseIdxVector {
     explicit SparseIdxVector(IdxVector indptr) : indptr_(indptr.empty() ? IdxVector{0} : indptr) {
         assert(!indptr.empty());
     }
+    SparseIdxVector(SparseIdxVector const&) = default;
+    SparseIdxVector(SparseIdxVector&&) = default;
+    SparseIdxVector& operator=(SparseIdxVector const&) = default;
+    SparseIdxVector& operator=(SparseIdxVector&&) = default;
+    ~SparseIdxVector() = default;
 
     constexpr auto size() const { return static_cast<Idx>(indptr_.size()) - 1; }
     constexpr auto begin() const { return GroupIterator<Idx>(indptr_, 0); }
     constexpr auto end() const { return GroupIterator<Idx>(indptr_, size()); }
 
     constexpr auto element_size() const { return indptr_.back(); }
-    auto get_element_range(Idx group) const { return boost::iterator_range<IdxCount>(indptr_[group], indptr_[group + 1]); }
+    auto get_element_range(Idx group) const {
+        return boost::iterator_range<IdxCount>(indptr_[group], indptr_[group + 1]);
+    }
     auto get_group(Idx element) const {
         assert(element < element_size());
         return static_cast<Idx>(std::upper_bound(indptr_.begin(), indptr_.end(), element) - indptr_.begin() - 1);
@@ -133,7 +140,7 @@ class DenseIdxVector {
 template <typename T>
 concept sparse_type = std::is_same<T, SparseIdxVector>::value;
 
-template <sparse_type First, sparse_type... Rest> auto zip_sequence(First& first, Rest&... rest) {
+template <sparse_type First, sparse_type... Rest> auto zip_sequence(First const& first, Rest const&... rest) {
 
     assert((first.size() == rest.size()) && ...);
     // TODO Add common index as count at the first postiion

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(SparseIdxVector const& source_buses, 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 source_number : source_buses.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, SparseIdxVector const& source_buses) {
+    for (Idx source : source_buses.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]);

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& source_buses = *this->source_buses_;
         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 : source_buses.get_element_range(bus_number)) {
                     // YBus_diag += Y_source
                     mat_data[data_sequence] += y_bus.math_model_param().source_param[source_number];
                 }
@@ -118,7 +117,7 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
     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_;
+        SparseIdxVector const& source_buses = *this->source_buses_;
         std::vector<LoadGenType> const& load_gen_type = *this->load_gen_type_;
 
         // set rhs to zero for iteration start
@@ -127,7 +126,7 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
         // 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_sources(bus_number, y_bus, input, source_buses);
         }
     }
 
@@ -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) {
+                     SparseIdxVector const& source_buses) {
+        for (Idx source_number : source_buses.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]});

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp

@@ -175,25 +175,24 @@ template <bool sym> class MeasuredValues {
         for (Idx bus = 0; bus != math_topology_->n_bus(); ++bus) {
             Idx const load_gen_begin = math_topology_->load_gen_bus_indptr[bus];
             Idx const load_gen_end = math_topology_->load_gen_bus_indptr[bus + 1];
-            Idx const source_begin = math_topology_->source_bus_indptr[bus];
-            Idx const source_end = math_topology_->source_bus_indptr[bus + 1];
+            auto const sources = math_topology_->source_buses.get_element_range(bus);
 
             // under-determined or exactly determined
             if (bus_injection_[bus].n_unmeasured_appliances > 0) {
-                calculate_non_over_determined_injection(
-                    bus_injection_[bus].n_unmeasured_appliances, load_gen_begin, load_gen_end, source_begin, source_end,
-                    bus_appliance_injection_[bus], s[bus], load_gen_flow, source_flow);
+                calculate_non_over_determined_injection(bus_injection_[bus].n_unmeasured_appliances, load_gen_begin,
+                                                        load_gen_end, sources, bus_appliance_injection_[bus], s[bus],
+                                                        load_gen_flow, source_flow);
             }
             // over-determined
             else {
-                calculate_over_determined_injection(load_gen_begin, load_gen_end, source_begin, source_end,
+                calculate_over_determined_injection(load_gen_begin, load_gen_end, sources,
                                                     bus_appliance_injection_[bus], s[bus], load_gen_flow, source_flow);
             }
             // current injection
             for (Idx load_gen = load_gen_begin; load_gen != load_gen_end; ++load_gen) {
                 load_gen_flow[load_gen].i = conj(load_gen_flow[load_gen].s / u[bus]);
             }
-            for (Idx source = source_begin; source != source_end; ++source) {
+            for (Idx source : sources) {
                 source_flow[source].i = conj(source_flow[source].s / u[bus]);
             }
         }
@@ -305,7 +304,7 @@ template <bool sym> class MeasuredValues {
             process_bus_objects(bus, topo.load_gen_bus_indptr, topo.load_gen_power_sensor_indptr, input.load_gen_status,
                                 input.measured_load_gen_power, extra_value_, idx_load_gen_power_);
             // source
-            process_bus_objects(bus, topo.source_bus_indptr, topo.source_power_sensor_indptr, input.source_status,
+            process_bus_objects(bus, topo.source_buses, topo.source_power_sensor_indptr, input.source_status,
                                 input.measured_source_power, extra_value_, idx_source_power_);
 
             combine_appliances_to_injection_measurements(input, topo, bus);
@@ -325,7 +324,7 @@ template <bool sym> class MeasuredValues {
             add_appliance_measurements(idx_load_gen_power_[load_gen], appliance_injection_measurement, n_unmeasured);
         }
 
-        for (Idx source = topo.source_bus_indptr[bus]; source != topo.source_bus_indptr[bus + 1]; ++source) {
+        for (Idx source : topo.source_buses.get_element_range(bus)) {
             add_appliance_measurements(idx_source_power_[source], appliance_injection_measurement, n_unmeasured);
         }
 
@@ -460,6 +459,16 @@ template <bool sym> class MeasuredValues {
         }
     }
 
+    // process objects in batch for shunt, load_gen, source
+    // return the status of the object type, if all the connected objects are measured
+    static void process_bus_objects(Idx const bus, SparseIdxVector const& objects, IdxVector const& sensor_indptr,
+                                    IntSVector const& obj_status, std::vector<SensorCalcParam<sym>> const& input_data,
+                                    std::vector<SensorCalcParam<sym>>& result_data, IdxVector& result_idx) {
+        for (Idx obj : objects.get_element_range(bus)) {
+            result_idx[obj] = process_one_object(obj, sensor_indptr, obj_status, input_data, result_data);
+        }
+    }
+
     // process one object
     static constexpr auto default_status_checker = [](auto x) -> bool { return x; };
     template <class TS, class StatusChecker = decltype(default_status_checker)>
@@ -497,20 +506,21 @@ template <bool sym> class MeasuredValues {
     }
 
     void calculate_non_over_determined_injection(Idx n_unmeasured, Idx load_gen_begin, Idx load_gen_end,
-                                                 Idx source_begin, Idx source_end,
+                                                 boost::iterator_range<IdxCount> sources,
                                                  SensorCalcParam<sym> const& bus_appliance_injection,
                                                  ComplexValue<sym> const& s, FlowVector& load_gen_flow,
                                                  FlowVector& source_flow) const {
         // calculate residual, divide, and assign to unmeasured (but connected) appliances
-        ComplexValue<sym> const s_residual_per_appliance = (s - bus_appliance_injection.value) / (double)n_unmeasured;
+        ComplexValue<sym> const s_residual_per_appliance =
+            (s - bus_appliance_injection.value) / static_cast<double>(n_unmeasured);
         for (Idx load_gen = load_gen_begin; load_gen != load_gen_end; ++load_gen) {
             if (has_load_gen(load_gen)) {
                 load_gen_flow[load_gen].s = load_gen_power(load_gen).value;
             } else if (idx_load_gen_power_[load_gen] == unmeasured) {
                 load_gen_flow[load_gen].s = s_residual_per_appliance;
             }
         }
-        for (Idx source = source_begin; source != source_end; ++source) {
+        for (Idx source : sources) {
             if (has_source(source)) {
                 source_flow[source].s = source_power(source).value;
             } else if (idx_source_power_[source] == unmeasured) {
@@ -519,7 +529,8 @@ template <bool sym> class MeasuredValues {
         }
     }
 
-    void calculate_over_determined_injection(Idx load_gen_begin, Idx load_gen_end, Idx source_begin, Idx source_end,
+    void calculate_over_determined_injection(Idx load_gen_begin, Idx load_gen_end,
+                                             boost::iterator_range<IdxCount> sources,
                                              SensorCalcParam<sym> const& bus_appliance_injection,
                                              ComplexValue<sym> const& s, FlowVector& load_gen_flow,
                                              FlowVector& source_flow) const {
@@ -532,7 +543,7 @@ template <bool sym> class MeasuredValues {
                 load_gen_flow[load_gen].s = load_gen_power(load_gen).value - (load_gen_power(load_gen).variance) * mu;
             }
         }
-        for (Idx source = source_begin; source != source_end; ++source) {
+        for (Idx source : sources) {
             if (has_source(source)) {
                 source_flow[source].s = source_power(source).value - (source_power(source).variance) * mu;
             }
@@ -608,7 +619,7 @@ template <bool sym> class IterativeLinearSESolver {
         main_timer.stop();
 
         const auto key = Timer::make_key(2228, "Max number of iterations");
-        calculation_info[key] = std::max(calculation_info[key], (double)num_iter);
+        calculation_info[key] = std::max(calculation_info[key], static_cast<double>(num_iter));
 
         return output;
     }

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_pf_solver.hpp

@@ -30,7 +30,7 @@ template <bool sym, typename DerivedSolver> class IterativePFSolver {
                                    Idx max_iter, CalculationInfo& calculation_info) {
         // get derived reference for derived solver class
         auto derived_solver = static_cast<DerivedSolver&>(*this);
-        IdxVector const& source_bus_indptr = *source_bus_indptr_;
+        auto const& source_buses = *source_buses_;
         std::vector<double> const& phase_shift = *phase_shift_;
 
         // prepare
@@ -47,7 +47,7 @@ template <bool sym, typename DerivedSolver> class IterativePFSolver {
             DoubleComplex const u_ref = [&]() {
                 DoubleComplex sum_u_ref = 0.0;
                 for (Idx bus = 0; bus != n_bus_; ++bus) {
-                    for (Idx source = source_bus_indptr[bus]; source != source_bus_indptr[bus + 1]; ++source) {
+                    for (Idx source : source_buses.get_element_range(bus)) {
                         sum_u_ref += input.source[source] * std::exp(1.0i * -phase_shift[bus]); // offset phase shift
                     }
                 }
@@ -109,13 +109,12 @@ template <bool sym, typename DerivedSolver> class IterativePFSolver {
         output.shunt = y_bus.template calculate_shunt_flow<ApplianceMathOutput<sym>>(output.u);
 
         // prepare source, load gen and bus_injection
-        output.source.resize(source_bus_indptr_->back());
+        output.source.resize(source_buses_->size());
         output.load_gen.resize(load_gen_bus_indptr_->back());
         output.bus_injection.resize(n_bus_);
 
         for (Idx bus_number = 0; bus_number != n_bus_; ++bus_number) {
-            common_solver_functions::calculate_source_result<sym>(bus_number, y_bus, input, output,
-                                                                  *source_bus_indptr_);
+            common_solver_functions::calculate_source_result<sym>(bus_number, y_bus, input, output, *source_buses_);
             common_solver_functions::calculate_load_gen_result<sym>(bus_number, input, output, *load_gen_bus_indptr_,
                                                                     [this](Idx i) { return (*load_gen_type_)[i]; });
         }
@@ -126,13 +125,13 @@ template <bool sym, typename DerivedSolver> class IterativePFSolver {
     Idx n_bus_;
     std::shared_ptr<DoubleVector const> phase_shift_;
     std::shared_ptr<IdxVector const> load_gen_bus_indptr_;
-    std::shared_ptr<IdxVector const> source_bus_indptr_;
+    std::shared_ptr<SparseIdxVector const> source_buses_;
     std::shared_ptr<std::vector<LoadGenType> const> load_gen_type_;
     IterativePFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
         : n_bus_{y_bus.size()},
           phase_shift_{topo_ptr, &topo_ptr->phase_shift},
           load_gen_bus_indptr_{topo_ptr, &topo_ptr->load_gen_bus_indptr},
-          source_bus_indptr_{topo_ptr, &topo_ptr->source_bus_indptr},
+          source_buses_{topo_ptr, &topo_ptr->source_buses},
           load_gen_type_{topo_ptr, &topo_ptr->load_gen_type} {}
 };