Feature: support multi-dimensional batch calculation[experimental]

power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/dataset.hpp

@@ -287,6 +287,9 @@ template <dataset_type_tag dataset_type_> class Dataset {
     }
     constexpr bool is_dense(Idx const i) const { return is_dense(buffers_[i]); }
     constexpr bool is_dense(Buffer const& buffer) const { return buffer.indptr.empty(); }
+    constexpr bool is_dense() const {
+        return std::ranges::all_of(buffers_, [this](Buffer const& buffer) { return is_dense(buffer); });
+    }
     constexpr bool is_sparse(std::string_view component, bool with_attribute_buffers = false) const {
         Idx const idx = find_component(component, false);
         if (idx == invalid_index) {
@@ -510,6 +513,35 @@ template <dataset_type_tag dataset_type_> class Dataset {
         return result;
     }
 
+    // get slice dataset from batch
+    Dataset get_slice_scenario(Idx begin, Idx end) const
+        requires(!is_indptr_mutable_v<dataset_type>)
+    {
+        assert(begin <= end);
+        assert(0 <= begin);
+        assert(end <= batch_size());
+        assert(is_batch());
+        assert(is_dense());
+
+        // empty slice
+        if (begin == end) {
+            Dataset result{true, 0, dataset_info_.dataset->name, *meta_data_};
+            result.add_buffer("node", 0, 0, nullptr, nullptr);
+            return result;
+        }
+
+        // start with begin
+        Dataset result = get_individual_scenario(begin);
+        Idx const batch_size = end - begin;
+        result.dataset_info_.is_batch = true;
+        result.dataset_info_.batch_size = batch_size;
+        for (auto&& [buffer, component_info] : std::views::zip(result.buffers_, result.dataset_info_.component_info)) {
+            Idx const size = component_info.elements_per_scenario * batch_size;
+            component_info.total_elements = size;
+        }
+        return result;
+    }
+
   private:
     MetaData const* meta_data_;
     DatasetInfo dataset_info_;

power_grid_model_c/power_grid_model_c/include/power_grid_model_c/basics.h

@@ -137,6 +137,12 @@ typedef struct PGM_WritableDataset PGM_WritableDataset;
  * @brief Opaque struct for the information of the dataset.
  */
 typedef struct PGM_DatasetInfo PGM_DatasetInfo;
+
+/**
+ * @brief Opaque struct for the multi dimensional dataset class.
+ *
+ */
+typedef struct PGM_MultiDimensionalDataset PGM_MultiDimensionalDataset;
 #endif
 
 // NOLINTEND(modernize-use-using)

power_grid_model_c/power_grid_model_c/include/power_grid_model_c/dataset.h

@@ -309,6 +309,36 @@ PGM_API void PGM_dataset_mutable_add_attribute_buffer(PGM_Handle* handle, PGM_Mu
  */
 PGM_API PGM_DatasetInfo const* PGM_dataset_mutable_get_info(PGM_Handle* handle, PGM_MutableDataset const* dataset);
 
+/**
+ * @brief Create a PGM_MultiDimensionalDataset from multiple PGM_ConstDataset instances
+ *
+ * @param handle
+ * @param const_datasets
+ * @param n_datasets
+ * @return
+ */
+PGM_API PGM_MultiDimensionalDataset*
+PGM_dataset_create_multidimensional_from_const(PGM_Handle* handle, PGM_ConstDataset const** const_datasets,
+                                               PGM_Idx n_datasets);
+
+/**
+ * @brief Get the array pointer from a PGM_MultiDimensionalDataset
+ *
+ * @param handle
+ * @param multidimensional_dataset
+ * @return
+ */
+PGM_API PGM_ConstDataset const*
+PGM_get_array_pointer_from_multidimensional(PGM_Handle* handle,
+                                            PGM_MultiDimensionalDataset const* multidimensional_dataset);
+
+/**
+ * @brief destroy the multidimensional dataset object
+ *
+ * @param multidimensional_dataset
+ */
+PGM_API void PGM_destroy_multidimensional_dataset(PGM_MultiDimensionalDataset* multidimensional_dataset);
+
 #ifdef __cplusplus
 }
 #endif

power_grid_model_c/power_grid_model_c/include/power_grid_model_c/options.h

@@ -121,6 +121,15 @@ PGM_API void PGM_set_short_circuit_voltage_scaling(PGM_Handle* handle, PGM_Optio
  */
 PGM_API void PGM_set_tap_changing_strategy(PGM_Handle* handle, PGM_Options* opt, PGM_Idx tap_changing_strategy);
 
+/**
+ * @brief Specify the batch dimension for batch calculations
+ *
+ * @param handle
+ * @param opt pointer to option instance
+ * @param batch_dimension dimension of the batch calculation
+ */
+PGM_API void PGM_set_batch_dimension(PGM_Handle* handle, PGM_Options* opt, PGM_Idx batch_dimension);
+
 /**
  * @brief Enable/disable experimental features.
  *

power_grid_model_c/power_grid_model_c/src/dataset.cpp

@@ -160,3 +160,28 @@
 PGM_DatasetInfo const* PGM_dataset_mutable_get_info(PGM_Handle* /*unused*/, PGM_MutableDataset const* dataset) {
     return &dataset->get_description();
 }
+
+PGM_MultiDimensionalDataset* PGM_dataset_create_multidimensional_from_const(PGM_Handle* handle,
+                                                                            PGM_ConstDataset const** const_datasets,
+                                                                            PGM_Idx n_datasets) {
+    return call_with_catch(
+        handle,
+        [const_datasets, n_datasets]() {
+            auto* multidimensional_dataset = new PGM_MultiDimensionalDataset();
+            for (PGM_Idx i = 0; i < n_datasets; ++i) {
+                multidimensional_dataset->emplace_back(*const_datasets[i]);
+            }
+            return multidimensional_dataset;
+        },
+        PGM_regular_error);
+}
+
+PGM_ConstDataset const*
+PGM_get_array_pointer_from_multidimensional(PGM_Handle* /*unused*/,
+                                            PGM_MultiDimensionalDataset const* multidimensional_dataset) {
+    return multidimensional_dataset->data();
+}
+
+void PGM_destroy_multidimensional_dataset(PGM_MultiDimensionalDataset* multidimensional_dataset) {
+    delete multidimensional_dataset;
+}

power_grid_model_c/power_grid_model_c/src/forward_declarations.hpp

@@ -10,6 +10,8 @@
 
 #include <power_grid_model/auxiliary/dataset_fwd.hpp>
 
+#include <vector>
+
 // forward declare all referenced struct/class in C++ core
 // alias them in the root namespace
 
@@ -36,3 +38,4 @@ using PGM_ConstDataset = power_grid_model::meta_data::Dataset<power_grid_model::
 using PGM_MutableDataset = power_grid_model::meta_data::Dataset<power_grid_model::mutable_dataset_t>;
 using PGM_WritableDataset = power_grid_model::meta_data::Dataset<power_grid_model::writable_dataset_t>;
 using PGM_DatasetInfo = power_grid_model::meta_data::DatasetInfo;
+using PGM_MultiDimensionalDataset = std::vector<PGM_ConstDataset>;

power_grid_model_c/power_grid_model_c/src/model.cpp

@@ -15,6 +15,9 @@
 #include <power_grid_model/common/common.hpp>
 #include <power_grid_model/main_model.hpp>
 
+#include <memory>
+#include <numeric>
+
 namespace {
 using namespace power_grid_model;
 } // namespace
@@ -55,6 +58,7 @@
         PGM_regular_error);
 }
 
+// helper functions
 namespace {
 void check_no_experimental_features_used(MainModel const& model, MainModel::Options const& opt) {
     // optionally add experimental feature checks here
@@ -142,9 +146,11 @@
 }
 } // namespace
 
-// run calculation
-void PGM_calculate(PGM_Handle* handle, PGM_PowerGridModel* model, PGM_Options const* opt,
-                   PGM_MutableDataset const* output_dataset, PGM_ConstDataset const* batch_dataset) {
+// calculation implementation
+namespace {
+
+void calculate_impl(PGM_Handle* handle, PGM_PowerGridModel* model, PGM_Options const* opt,
+                    PGM_MutableDataset const* output_dataset, PGM_ConstDataset const* batch_dataset) {
     PGM_clear_error(handle);
     // check dataset integrity
     if ((batch_dataset != nullptr) && (!batch_dataset->is_batch() || !output_dataset->is_batch())) {
@@ -171,14 +177,92 @@
         handle->err_msg = e.what();
         handle->failed_scenarios = e.failed_scenarios();
         handle->batch_errs = e.err_msgs();
     } catch (std::exception& e) {
         handle->err_code = PGM_regular_error;
         handle->err_msg = e.what();
     } catch (...) {
         handle->err_code = PGM_regular_error;
         handle->err_msg = "Unknown error!\n";
     }
 }
 
+void merge_batch_error_msgs(PGM_Handle* handle, PGM_Handle const& local_handle, Idx scenario_offset, Idx stride_size) {
+    if (local_handle.err_code == PGM_no_error) {
+        return;
+    }
+    handle->err_code = PGM_batch_error;
+    if (local_handle.err_code == PGM_batch_error) {
+        for (auto&& [idx, err_msg] : std::views::zip(local_handle.failed_scenarios, local_handle.batch_errs)) {
+            handle->failed_scenarios.push_back(idx + scenario_offset);
+            handle->batch_errs.push_back(err_msg);
+        }
+    } else {
+        for (Idx i = 0; i < stride_size; ++i) {
+            handle->failed_scenarios.push_back(scenario_offset + i);
+            handle->batch_errs.push_back(local_handle.err_msg);
+        }
+    }
+}
+
+} // namespace
+
+// run calculation
+void PGM_calculate(PGM_Handle* handle, PGM_PowerGridModel* model, PGM_Options const* opt,
+                   PGM_MutableDataset const* output_dataset, PGM_ConstDataset const* batch_dataset) {
+    // if dimension is zero, no batch calculation, force pointer to NULL
+    if (opt->batch_dimension == 0) {
+        batch_dataset = nullptr;
+    }
+
+    // for dimensions which are 1D batch or default (-1), call implementation directly
+    if (opt->batch_dimension < 2) {
+        calculate_impl(handle, model, opt, output_dataset, batch_dataset);
+        return;
+    }
+
+    // get stride size of the rest of dimensions
+    Idx const first_batch_size = batch_dataset->batch_size();
+    Idx const stride_size =
+        std::transform_reduce(batch_dataset + 1, batch_dataset + opt->batch_dimension, Idx{1}, std::multiplies{},
+                              [](PGM_ConstDataset const& ds) { return ds.batch_size(); });
+
+    // loop over the first dimension batche
+    for (Idx i = 0; i < first_batch_size; ++i) {
+        // a new handle
+        PGM_Handle local_handle{};
+        // deep opt is one dimension less
+        PGM_Options deep_opt = *opt;
+        --deep_opt.batch_dimension;
+        // create sliced datasets for the rest of dimensions
+        PGM_ConstDataset const single_update_dataset = batch_dataset->get_individual_scenario(i);
+        PGM_MutableDataset const sliced_output_dataset =
+            output_dataset->get_slice_scenario(i * stride_size, (i + 1) * stride_size);
+
+        // create a model copy
+        std::unique_ptr<PGM_PowerGridModel> const local_model{PGM_copy_model(&local_handle, model)};
+        if (local_handle.err_code != PGM_no_error) {
+            merge_batch_error_msgs(handle, local_handle, i * stride_size, stride_size);
+            continue;
+        }
+
+        // apply the update
+        PGM_update_model(&local_handle, local_model.get(), &single_update_dataset);
+        if (local_handle.err_code != PGM_no_error) {
+            merge_batch_error_msgs(handle, local_handle, i * stride_size, stride_size);
+            continue;
+        }
+
+        // if the deep opt have less than 2 dimensions, call implementation directly
+        if (deep_opt.batch_dimension < 2) {
+            calculate_impl(&local_handle, local_model.get(), &deep_opt, &sliced_output_dataset, batch_dataset + 1);
+        } else {
+            // recursive call
+            PGM_calculate(&local_handle, local_model.get(), &deep_opt, &sliced_output_dataset, batch_dataset + 1);
+        }
+        // merge errors
+        merge_batch_error_msgs(handle, local_handle, i * stride_size, stride_size);
+    }
+}
+
 // destroy model
 void PGM_destroy_model(PGM_PowerGridModel* model) { delete model; }

power_grid_model_c/power_grid_model_c/src/options.cpp

@@ -33,6 +33,9 @@ void PGM_set_short_circuit_voltage_scaling(PGM_Handle* /* handle */, PGM_Options
 void PGM_set_tap_changing_strategy(PGM_Handle* /* handle */, PGM_Options* opt, PGM_Idx tap_changing_strategy) {
     opt->tap_changing_strategy = tap_changing_strategy;
 }
+void PGM_set_batch_dimension(PGM_Handle* /* handle */, PGM_Options* opt, PGM_Idx batch_dimension) {
+    opt->batch_dimension = batch_dimension;
+}
 void PGM_set_experimental_features(PGM_Handle* /* handle */, PGM_Options* opt, PGM_Idx experimental_features) {
     opt->experimental_features = experimental_features;
 }

power_grid_model_c/power_grid_model_c/src/options.hpp

@@ -25,5 +25,6 @@ struct PGM_Options {
     Idx threading{-1};
     Idx short_circuit_voltage_scaling{PGM_short_circuit_voltage_scaling_maximum};
     Idx tap_changing_strategy{PGM_tap_changing_strategy_disabled};
+    Idx batch_dimension{-1};
     Idx experimental_features{PGM_experimental_features_disabled};
 };

power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/model.hpp

@@ -13,6 +13,8 @@
 
 #include "power_grid_model_c/model.h"
 
+#include <vector>
+
 namespace power_grid_model_cpp {
 class Model {
   public:
@@ -54,6 +56,25 @@ class Model {
         handle_.call_with(PGM_calculate, get(), opt.get(), output_dataset.get(), nullptr);
     }
 
+    void calculate(Options const& opt, DatasetMutable const& output_dataset,
+                   std::vector<DatasetConst> const& batch_datasets) {
+        // create multidimensional dataset from the span of datasets
+        std::vector<PGM_ConstDataset const*> dataset_ptrs;
+        dataset_ptrs.reserve(batch_datasets.size());
+        for (auto const& ds : batch_datasets) {
+            dataset_ptrs.push_back(ds.get());
+        }
+        auto multidimensional_dataset =
+            detail::UniquePtr<PGM_MultiDimensionalDataset, &PGM_destroy_multidimensional_dataset>{
+                handle_.call_with(PGM_dataset_create_multidimensional_from_const, dataset_ptrs.data(),
+                                  static_cast<Idx>(dataset_ptrs.size()))};
+        RawConstDataset const* batch_dataset_array_pointer =
+            PGM_get_array_pointer_from_multidimensional(nullptr, multidimensional_dataset.get());
+
+        // call calculate with the multidimensional dataset
+        handle_.call_with(PGM_calculate, get(), opt.get(), output_dataset.get(), batch_dataset_array_pointer);
+    }
+
   private:
     Handle handle_{};
     detail::UniquePtr<PowerGridModel, &PGM_destroy_model> model_;

power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/options.hpp

@@ -39,6 +39,10 @@ class Options {
         handle_.call_with(PGM_set_tap_changing_strategy, get(), tap_changing_strategy);
     }
 
+    void set_batch_dimension(Idx batch_dimension) {
+        handle_.call_with(PGM_set_batch_dimension, get(), batch_dimension);
+    }
+
     void set_experimental_features(Idx experimental_features) {
         handle_.call_with(PGM_set_experimental_features, get(), experimental_features);
     }

src/power_grid_model/_core/options.py

@@ -45,6 +45,7 @@ class Options:
     tap_changing_strategy = OptionSetter(pgc.set_tap_changing_strategy)
     short_circuit_voltage_scaling = OptionSetter(pgc.set_short_circuit_voltage_scaling)
     experimental_features = OptionSetter(pgc.set_experimental_features)
+    batch_dimension = OptionSetter(pgc.set_batch_dimension)
 
     @property
     def opt(self) -> OptionsPtr:

src/power_grid_model/_core/power_grid_core.py

@@ -87,6 +87,12 @@ class WritableDatasetPtr(c_void_p):
     """
 
 
+class MultiDimensionDatasetPtr(c_void_p):
+    """
+    Pointer to multi-dimensional dataset
+    """
+
+
 class DatasetInfoPtr(c_void_p):
     """
     Pointer to dataset info
@@ -338,6 +344,10 @@ def set_max_iter(self, opt: OptionsPtr, max_iter: int) -> None:  # type: ignore[
     def set_threading(self, opt: OptionsPtr, threading: int) -> None:  # type: ignore[empty-body]
         pass  # pragma: no cover
 
+    @make_c_binding
+    def set_batch_dimension(self, opt: OptionsPtr, batch_dimension: int) -> None:  # type: ignore[empty-body]
+        pass  # pragma: no cover
+
     @make_c_binding
     def create_model(  # type: ignore[empty-body]
         self,
@@ -505,6 +515,28 @@ def dataset_writable_set_attribute_buffer(
     ) -> None:  # type: ignore[empty-body]
         pass  # pragma: no cover
 
+    @make_c_binding
+    def dataset_create_multidimensional_from_const(  # type: ignore[empty-body]
+        self,
+        const_datasets: POINTER(ConstDatasetPtr),  # type: ignore[valid-type]
+        n_datasets: int,
+    ) -> MultiDimensionDatasetPtr:
+        pass  # pragma: no cover
+
+    @make_c_binding
+    def get_array_pointer_from_multidimensional(  # type: ignore[empty-body]
+        self,
+        multidimensional_dataset: MultiDimensionDatasetPtr,
+    ) -> ConstDatasetPtr:  # type: ignore[valid-type]
+        pass  # pragma: no cover
+
+    @make_c_binding
+    def destroy_multidimensional_dataset(  # type: ignore[empty-body]
+        self,
+        multidimensional_dataset: MultiDimensionDatasetPtr,
+    ) -> None:
+        pass  # pragma: no cover
+
     @make_c_binding
     def create_deserializer_from_binary_buffer(  # type: ignore[empty-body]
         self, data: bytes, size: int, serialization_format: int