Reduce usage of shared_ptr in the core

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

@@ -80,17 +80,17 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym_type, IterativeCur
 
     static constexpr auto is_iterative = true;
 
-    IterativeCurrentPFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
+    IterativeCurrentPFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
         : IterativePFSolver<sym, IterativeCurrentPFSolver>{y_bus, topo_ptr},
           rhs_u_(y_bus.size()),
-          sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()} {}
+          sparse_solver_{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()} {}
 
     // Add source admittance to Y bus and set variable for prepared y bus to true
     void initialize_derived_solver(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
                                    SolverOutput<sym>& output) {
         make_flat_start(input, output.u);
 
-        auto const& sources_per_bus = *this->sources_per_bus_;
+        auto const& sources_per_bus = this->sources_per_bus_.get();
         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
@@ -119,14 +119,14 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym_type, IterativeCur
     // Prepare matrix calculates injected current, i.e., RHS of solver for each iteration.
     void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
                                 ComplexValueVector<sym> const& u) {
-        std::vector<LoadGenType> const& load_gen_type = *this->load_gen_type_;
+        std::vector<LoadGenType> const& load_gen_type = this->load_gen_type_.get();
 
         // set rhs to zero for iteration start
         std::fill(rhs_u_.begin(), rhs_u_.end(), ComplexValue<sym>{0.0});
 
         // loop buses: i
         for (auto const& [bus_number, load_gens, sources] :
-             enumerated_zip_sequence(*this->load_gens_per_bus_, *this->sources_per_bus_)) {
+             enumerated_zip_sequence(this->load_gens_per_bus_.get(), this->sources_per_bus_.get())) {
             add_loads(load_gens, bus_number, input, load_gen_type, u);
             add_sources(sources, bus_number, y_bus, input);
         }
@@ -197,11 +197,11 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym_type, IterativeCur
     }
 
     void make_flat_start(PowerFlowInput<sym> const& input, ComplexValueVector<sym>& output_u) {
-        std::vector<double> const& phase_shift = *this->phase_shift_;
+        std::vector<double> const& phase_shift = this->phase_shift_.get();
         // average u_ref of all sources
         DoubleComplex const u_ref = [&]() {
             DoubleComplex sum_u_ref = 0.0;
-            for (auto const& [bus, sources] : enumerated_zip_sequence(*this->sources_per_bus_)) {
+            for (auto const& [bus, sources] : enumerated_zip_sequence(this->sources_per_bus_.get())) {
                 for (Idx const source : sources) {
                     sum_u_ref += input.source[source] * std::exp(1.0i * -phase_shift[bus]); // offset phase shift
                 }

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

@@ -75,12 +75,12 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
     static constexpr Idx bsr_block_size_ = is_symmetric_v<sym> ? 2 : 6;
 
   public:
-    IterativeLinearSESolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> topo_ptr)
+    IterativeLinearSESolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
         : n_bus_{y_bus.size()},
-          math_topo_{std::move(topo_ptr)},
+          math_topo_{topo_ptr},
           data_gain_(y_bus.nnz_lu()),
           x_rhs_(y_bus.size()),
-          sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()},
+          sparse_solver_{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()},
           perm_(y_bus.size()) {}
 
     SolverOutput<sym> run_state_estimation(YBus<sym> const& y_bus, StateEstimationInput<sym> const& input,
@@ -111,7 +111,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
         sub_timer = Timer{log, LogEvent::initialize_voltages}; // TODO(mgovers): make scoped subtimers
         RealValue<sym> const mean_angle_shift = measured_values.mean_angle_shift();
         for (Idx bus = 0; bus != n_bus_; ++bus) {
-            output.u[bus] = exp(1.0i * (mean_angle_shift + math_topo_->phase_shift[bus]));
+            output.u[bus] = exp(1.0i * (mean_angle_shift + math_topo_.get().phase_shift[bus]));
         }
 
         // loop to iterate
@@ -155,7 +155,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
 
     Idx n_bus_;
     // shared topo data
-    std::shared_ptr<MathModelTopology const> math_topo_;
+    std::reference_wrapper<MathModelTopology const> math_topo_;
 
     // data for gain matrix
     std::vector<ILSEGainBlock<sym>> data_gain_;
@@ -359,7 +359,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
             if (has_angle) {
                 return 1.0;
             }
-            auto const& voltage = x_rhs_[math_topo_->slack_bus].u();
+            auto const& voltage = x_rhs_[math_topo_.get().slack_bus].u();
             auto const& voltage_a = [&voltage]() -> auto const& {
                 if constexpr (is_symmetric_v<sym>) {
                     return voltage;

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

@@ -81,22 +81,22 @@ template <symmetry_tag sym, typename DerivedSolver> class IterativePFSolver {
     }
 
     void calculate_result(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, SolverOutput<sym>& output) {
-        detail::calculate_pf_result(y_bus, input, *sources_per_bus_, *load_gens_per_bus_, output,
-                                    [this](Idx i) { return (*load_gen_type_)[i]; });
+        detail::calculate_pf_result(y_bus, input, sources_per_bus_.get(), load_gens_per_bus_.get(), output,
+                                    [this](Idx i) { return (load_gen_type_.get())[i]; });
     }
 
   private:
     Idx n_bus_;
-    std::shared_ptr<DoubleVector const> phase_shift_;
-    std::shared_ptr<SparseGroupedIdxVector const> load_gens_per_bus_;
-    std::shared_ptr<DenseGroupedIdxVector const> sources_per_bus_;
-    std::shared_ptr<std::vector<LoadGenType> const> load_gen_type_;
-    IterativePFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
+    std::reference_wrapper<DoubleVector const> phase_shift_;
+    std::reference_wrapper<SparseGroupedIdxVector const> load_gens_per_bus_;
+    std::reference_wrapper<DenseGroupedIdxVector const> sources_per_bus_;
+    std::reference_wrapper<std::vector<LoadGenType> const> load_gen_type_;
+    IterativePFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
         : n_bus_{y_bus.size()},
-          phase_shift_{topo_ptr, &topo_ptr->phase_shift},
-          load_gens_per_bus_{topo_ptr, &topo_ptr->load_gens_per_bus},
-          sources_per_bus_{topo_ptr, &topo_ptr->sources_per_bus},
-          load_gen_type_{topo_ptr, &topo_ptr->load_gen_type} {}
+          phase_shift_{std::cref(topo_ptr.phase_shift)},
+          load_gens_per_bus_{std::cref(topo_ptr.load_gens_per_bus)},
+          sources_per_bus_{std::cref(topo_ptr.sources_per_bus)},
+          load_gen_type_{std::cref(topo_ptr.load_gen_type)} {}
 };
 
 } // namespace power_grid_model::math_solver

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/linear_pf_solver.hpp

@@ -54,12 +54,12 @@ template <symmetry_tag sym_type> class LinearPFSolver {
 
     static constexpr auto is_iterative = false;
 
-    LinearPFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
+    LinearPFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
         : n_bus_{y_bus.size()},
-          load_gens_per_bus_{topo_ptr, &topo_ptr->load_gens_per_bus},
-          sources_per_bus_{topo_ptr, &topo_ptr->sources_per_bus},
+          load_gens_per_bus_{std::cref(topo_ptr.load_gens_per_bus)},
+          sources_per_bus_{std::cref(topo_ptr.sources_per_bus)},
           mat_data_(y_bus.nnz_lu()),
-          sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()},
+          sparse_solver_{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()},
           perm_(n_bus_) {}
 
     SolverOutput<sym> run_power_flow(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, Logger& log) {
@@ -92,20 +92,21 @@ template <symmetry_tag sym_type> class LinearPFSolver {
   private:
     Idx n_bus_;
     // shared topo data
-    std::shared_ptr<SparseGroupedIdxVector const> load_gens_per_bus_;
-    std::shared_ptr<DenseGroupedIdxVector const> sources_per_bus_;
+    std::reference_wrapper<SparseGroupedIdxVector const> load_gens_per_bus_;
+    std::reference_wrapper<DenseGroupedIdxVector const> sources_per_bus_;
     // sparse linear equation
     ComplexTensorVector<sym> mat_data_;
     // sparse solver
     SparseSolverType sparse_solver_;
     BlockPermArray perm_;
 
     void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, SolverOutput<sym>& output) {
-        detail::prepare_linear_matrix_and_rhs(y_bus, input, *load_gens_per_bus_, *sources_per_bus_, output, mat_data_);
+        detail::prepare_linear_matrix_and_rhs(y_bus, input, load_gens_per_bus_.get(), sources_per_bus_.get(), output,
+                                              mat_data_);
     }
 
     void calculate_result(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, SolverOutput<sym>& output) {
-        detail::calculate_pf_result(y_bus, input, *sources_per_bus_, *load_gens_per_bus_, output,
+        detail::calculate_pf_result(y_bus, input, sources_per_bus_.get(), load_gens_per_bus_.get(), output,
                                     [](Idx /*i*/) { return LoadGenType::const_y; });
     }
 };

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

@@ -87,7 +87,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
         // construct model if needed
         if (!iec60909_sc_solver_.has_value()) {
             Timer const timer{log, LogEvent::create_math_solver};
-            iec60909_sc_solver_.emplace(y_bus, topo_ptr_);
+            iec60909_sc_solver_.emplace(y_bus, *topo_ptr_);
         }
 
         // call calculation
@@ -121,7 +121,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
                                                     Logger& log, YBus<sym> const& y_bus) {
         if (!newton_raphson_pf_solver_.has_value()) {
             Timer const timer{log, LogEvent::create_math_solver};
-            newton_raphson_pf_solver_.emplace(y_bus, topo_ptr_);
+            newton_raphson_pf_solver_.emplace(y_bus, *topo_ptr_);
         }
         return newton_raphson_pf_solver_.value().run_power_flow(y_bus, input, err_tol, max_iter, log);
     }
@@ -130,7 +130,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
                                             Logger& log, YBus<sym> const& y_bus) {
         if (!linear_pf_solver_.has_value()) {
             Timer const timer{log, LogEvent::create_math_solver};
-            linear_pf_solver_.emplace(y_bus, topo_ptr_);
+            linear_pf_solver_.emplace(y_bus, *topo_ptr_);
         }
         return linear_pf_solver_.value().run_power_flow(y_bus, input, log);
     }
@@ -139,7 +139,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
                                                        Logger& log, YBus<sym> const& y_bus) {
         if (!iterative_current_pf_solver_.has_value()) {
             Timer const timer{log, LogEvent::create_math_solver};
-            iterative_current_pf_solver_.emplace(y_bus, topo_ptr_);
+            iterative_current_pf_solver_.emplace(y_bus, *topo_ptr_);
         }
         return iterative_current_pf_solver_.value().run_power_flow(y_bus, input, err_tol, max_iter, log);
     }
@@ -154,7 +154,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
         // construct model if needed
         if (!iterative_linear_se_solver_.has_value()) {
             Timer const timer{log, LogEvent::create_math_solver};
-            iterative_linear_se_solver_.emplace(y_bus, topo_ptr_);
+            iterative_linear_se_solver_.emplace(y_bus, *topo_ptr_);
         }
 
         // call calculation
@@ -166,7 +166,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
         // construct model if needed
         if (!newton_raphson_se_solver_.has_value()) {
             Timer const timer{log, LogEvent::create_math_solver};
-            newton_raphson_se_solver_.emplace(y_bus, topo_ptr_);
+            newton_raphson_se_solver_.emplace(y_bus, *topo_ptr_);
         }
 
         // call calculation

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp

@@ -36,7 +36,7 @@ template <symmetry_tag sym> class MeasuredValues {
   public:
     // construct
     MeasuredValues(std::shared_ptr<MathModelTopology const> topo, StateEstimationInput<sym> const& input)
-        : math_topology_{std::move(topo)},
+        : math_topology_{*topo},
           bus_appliance_injection_(math_topology().n_bus()),
           idx_voltage_(math_topology().n_bus()),
           bus_injection_(math_topology().n_bus()),
@@ -112,12 +112,12 @@ template <symmetry_tag sym> class MeasuredValues {
 
     LoadGenSourceFlow calculate_load_gen_source(ComplexValueVector<sym> const& u,
                                                 ComplexValueVector<sym> const& s) const {
-        std::vector<ApplianceSolverOutput<sym>> load_gen_flow(math_topology_->n_load_gen());
-        std::vector<ApplianceSolverOutput<sym>> source_flow(math_topology_->n_source());
+        std::vector<ApplianceSolverOutput<sym>> load_gen_flow(math_topology_.n_load_gen());
+        std::vector<ApplianceSolverOutput<sym>> source_flow(math_topology_.n_source());
 
         // loop all buses
         for (auto const& [bus, load_gens, sources] :
-             enumerated_zip_sequence(math_topology_->load_gens_per_bus, math_topology_->sources_per_bus)) {
+             enumerated_zip_sequence(math_topology_.load_gens_per_bus, math_topology_.sources_per_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_gens, sources,
@@ -173,7 +173,7 @@ template <symmetry_tag sym> class MeasuredValues {
 
   private:
     // cache topology
-    std::shared_ptr<MathModelTopology const> math_topology_;
+    MathModelTopology const& math_topology_;
 
     // flat arrays of all the relevant measurement for the main calculation
     // branch/shunt flow, bus voltage, injection flow
@@ -214,7 +214,7 @@ template <symmetry_tag sym> class MeasuredValues {
     // the lowest bus index with a voltage measurement
     Idx first_voltage_measurement_{};
 
-    constexpr MathModelTopology const& math_topology() const { return *math_topology_; }
+    constexpr MathModelTopology const& math_topology() const { return math_topology_; }
 
     void process_bus_related_measurements(StateEstimationInput<sym> const& input) {
         /*

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_pf_solver.hpp

@@ -211,12 +211,12 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF
 
     static constexpr auto is_iterative = true;
 
-    NewtonRaphsonPFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
+    NewtonRaphsonPFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
         : IterativePFSolver<sym, NewtonRaphsonPFSolver>{y_bus, topo_ptr},
           data_jac_(y_bus.nnz_lu()),
           x_(y_bus.size()),
           del_x_pq_(y_bus.size()),
-          sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()},
+          sparse_solver_{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()},
           perm_(y_bus.size()) {}
 
     // Initilize the unknown variable in polar form
@@ -225,13 +225,12 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF
         using LinearSparseSolverType = SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>>;
 
         ComplexTensorVector<sym> linear_mat_data(y_bus.nnz_lu());
-        LinearSparseSolverType linear_sparse_solver{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(),
-                                                    y_bus.shared_diag_lu()};
+        LinearSparseSolverType linear_sparse_solver{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()};
         typename LinearSparseSolverType::BlockPermArray linear_perm(y_bus.size());
 
         detail::copy_y_bus<sym>(y_bus, linear_mat_data);
-        detail::prepare_linear_matrix_and_rhs(y_bus, input, *this->load_gens_per_bus_, *this->sources_per_bus_, output,
-                                              linear_mat_data);
+        detail::prepare_linear_matrix_and_rhs(y_bus, input, this->load_gens_per_bus_.get(),
+                                              this->sources_per_bus_.get(), output, linear_mat_data);
         linear_sparse_solver.prefactorize_and_solve(linear_mat_data, linear_perm, output.u, output.u);
 
         // get magnitude and angle of start voltage
@@ -244,13 +243,13 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF
     // Calculate the Jacobian and deviation
     void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
                                 ComplexValueVector<sym> const& u) {
-        std::vector<LoadGenType> const& load_gen_type = *this->load_gen_type_;
+        std::vector<LoadGenType> const& load_gen_type = this->load_gen_type_.get();
         IdxVector const& bus_entry = y_bus.lu_diag();
 
         prepare_matrix_and_rhs_from_network_perspective(y_bus, u, bus_entry);
 
         for (auto const& [bus_number, load_gens, sources] :
-             enumerated_zip_sequence(*this->load_gens_per_bus_, *this->sources_per_bus_)) {
+             enumerated_zip_sequence(this->load_gens_per_bus_.get(), this->sources_per_bus_.get())) {
             Idx const diagonal_position = bus_entry[bus_number];
             add_loads(load_gens, bus_number, diagonal_position, input, load_gen_type);
             add_sources(sources, bus_number, diagonal_position, y_bus, input, u);