From c5a991c4c98fdc91e20a44f1c9aa0e886cfb1ee9 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Thu, 17 Aug 2023 10:51:54 +0200 Subject: [PATCH] feat: implement SEXP export of Field structure feat: implement NamedVector based on Rcpp::NumericVector feat: remove hard coded checks and substitute by R hook functions, defined in the input script refactor: modify API of DHT_Wrapper/InterpolationModule to expect a work package structure, where input/output values are stored as 2D vectors test: add tests for NamedVector test: extend tests for Field --- app/poet.cpp | 37 +-- bench/dolo/dolo_diffu_inner.R | 106 ++++++--- bench/dolo/dolo_diffu_inner_large.R | 108 ++++++--- bench/dolo/dolo_interp_long.R | 85 +++++-- include/poet/ChemistryModule.hpp | 5 +- include/poet/DHT_Wrapper.hpp | 17 +- include/poet/DataStructures.hpp | 199 ++++++++++++++-- include/poet/DiffusionModule.hpp | 5 +- include/poet/Field.hpp | 159 ------------- include/poet/Interpolation.hpp | 19 +- include/poet/LookupKey.hpp | 4 +- include/poet/RInsidePOET.hpp | 54 +++-- include/poet/SimParams.hpp | 14 +- src/ChemistryModule/ChemistryModule.cpp | 22 +- src/ChemistryModule/MasterFunctions.cpp | 4 +- .../SurrogateModels/DHT_Wrapper.cpp | 52 ++++- .../SurrogateModels/InterpolationModule.cpp | 84 +++---- .../SurrogateModels/ProximityHashTable.cpp | 6 +- src/ChemistryModule/WorkerFunctions.cpp | 15 +- src/DataStructures/Field.cpp | 215 +++++++++++++----- src/DiffusionModule.cpp | 53 ++--- src/SimParams.cpp | 16 +- test/CMakeLists.txt | 4 + test/RInsidePOET_funcs.R | 22 ++ test/testDataStructures.hpp.in | 6 + .../{testDataStructures.cpp => testField.cpp} | 62 ++++- test/testNamedVector.cpp | 63 +++++ 27 files changed, 931 insertions(+), 505 deletions(-) delete mode 100644 include/poet/Field.hpp create mode 100644 test/RInsidePOET_funcs.R create mode 100644 test/testDataStructures.hpp.in rename test/{testDataStructures.cpp => testField.cpp} (71%) create mode 100644 test/testNamedVector.cpp diff --git a/app/poet.cpp b/app/poet.cpp index 51dc51fd9..4095cf858 100644 --- a/app/poet.cpp +++ b/app/poet.cpp @@ -21,12 +21,12 @@ #include #include #include -#include +#include #include #include +#include #include #include -#include #include #include @@ -153,13 +153,15 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, /* displaying iteration number, with C++ and R iterator */ MSG("Going through iteration " + std::to_string(iter)); - MSG("R's $iter: " + std::to_string((uint32_t)(R.parseEval("mysetup$iter"))) + ". Iteration"); + MSG("R's $iter: " + + std::to_string((uint32_t)(R.parseEval("mysetup$iter"))) + + ". Iteration"); /* run transport */ // TODO: transport to diffusion diffusion.simulate(dt); - chem.getField().UpdateFromField(diffusion.getField()); + chem.getField().update(diffusion.getField()); MSG("Chemistry step"); @@ -167,7 +169,7 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, chem.RunCells(); writeFieldsToR(R, diffusion.getField(), chem.GetField()); - diffusion.getField().UpdateFromField(chem.GetField()); + diffusion.getField().update(chem.GetField()); R["req_dt"] = dt; R["simtime"] = (sim_time += dt); @@ -181,7 +183,8 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, // store_result is TRUE) R.parseEvalQ("mysetup <- master_iteration_end(setup=mysetup)"); - MSG("End of *coupling* iteration " + std::to_string(iter) + "/" + std::to_string(maxiter)); + MSG("End of *coupling* iteration " + std::to_string(iter) + "/" + + std::to_string(maxiter)); MSG(); // MPI_Barrier(MPI_COMM_WORLD); @@ -258,15 +261,16 @@ int main(int argc, char *argv[]) { MPI_Comm_rank(MPI_COMM_WORLD, &world_rank); + RInsidePOET &R = RInsidePOET::getInstance(); + if (world_rank == 0) { - MSG("Running POET version " + std::string(poet_version)); + MSG("Running POET version " + std::string(poet_version)); } if (world_rank > 0) { { - SimParams params(world_rank, world_size); - int pret = params.parseFromCmdl(argv, RInsidePOET::getInstance()); + int pret = params.parseFromCmdl(argv, R); if (pret == poet::PARSER_ERROR) { MPI_Finalize(); @@ -286,7 +290,7 @@ int main(int argc, char *argv[]) { } MPI_Barrier(MPI_COMM_WORLD); - + MSG("finished, cleanup of process " + std::to_string(world_rank)); MPI_Finalize(); @@ -294,8 +298,6 @@ int main(int argc, char *argv[]) { return EXIT_SUCCESS; } - RInsidePOET &R = RInsidePOET::getInstance(); - /*Loading Dependencies*/ // TODO: kann raus std::string r_load_dependencies = "source('../R_lib/kin_r_library.R');"; @@ -326,12 +328,12 @@ int main(int argc, char *argv[]) { // MDL: store all parameters if (world_rank == 0) { - MSG("Calling R Function to store calling parameters"); - R.parseEvalQ("StoreSetup(setup=mysetup)"); + MSG("Calling R Function to store calling parameters"); + R.parseEvalQ("StoreSetup(setup=mysetup)"); } if (world_rank == 0) { - MSG("Init done on process with rank " + std::to_string(world_rank)); + MSG("Init done on process with rank " + std::to_string(world_rank)); } // MPI_Barrier(MPI_COMM_WORLD); @@ -351,7 +353,8 @@ int main(int argc, char *argv[]) { r_vis_code = "saveRDS(profiling, file=paste0(fileout,'/timings.rds'));"; R.parseEval(r_vis_code); - MSG("Done! Results are stored as R objects into <" + params.getOutDir() + "/timings.rds>"); + MSG("Done! Results are stored as R objects into <" + params.getOutDir() + + "/timings.rds>"); MPI_Barrier(MPI_COMM_WORLD); @@ -359,7 +362,7 @@ int main(int argc, char *argv[]) { MPI_Finalize(); if (world_rank == 0) { - MSG("done, bye!"); + MSG("done, bye!"); } exit(0); diff --git a/bench/dolo/dolo_diffu_inner.R b/bench/dolo/dolo_diffu_inner.R index 071fe36c1..a16f8a49c 100644 --- a/bench/dolo/dolo_diffu_inner.R +++ b/bench/dolo/dolo_diffu_inner.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2023-08-02 13:59:02 mluebke" +## Time-stamp: "Last modified 2023-08-16 17:04:42 mluebke" database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat") input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi") @@ -62,34 +62,34 @@ alpha_diffu <- c( ## list of boundary conditions/inner nodes vecinj_diffu <- list( - list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0, - "Cl" = 0.002, - "Mg" = 0.001 - ), - list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0.0, - "Cl" = 0.004, - "Mg" = 0.002 - ) + list( + "H" = 110.683, + "O" = 55.3413, + "Charge" = 1.90431e-16, + "C(4)" = 0, + "Ca" = 0, + "Cl" = 0.002, + "Mg" = 0.001 + ), + list( + "H" = 110.683, + "O" = 55.3413, + "Charge" = 1.90431e-16, + "C(4)" = 0, + "Ca" = 0.0, + "Cl" = 0.004, + "Mg" = 0.002 + ) ) vecinj_inner <- list( - l1 = c(1,20,20), - l2 = c(2,80,80), - l3 = c(2,60,80) + l1 = c(1, 20, 20), + l2 = c(2, 80, 80), + l3 = c(2, 60, 80) ) boundary <- list( -# "N" = c(1, rep(0, n-1)), + # "N" = c(1, rep(0, n-1)), "N" = rep(0, n), "E" = rep(0, n), "S" = rep(0, n), @@ -116,21 +116,59 @@ diffusion <- list( ## # Needed when using DHT -dht_species <- c("H" = 10, - "O" = 10, - "Charge" = 3, - "C(4)" = 5, - "Ca" = 5, - "Cl" = 5, - "Mg" = 5, - "Calcite" = 5, - "Dolomite" =5 - ) +dht_species <- c( + "H" = 10, + "O" = 10, + "Charge" = 3, + "C(4)" = 5, + "Ca" = 5, + "Cl" = 5, + "Mg" = 5, + "Calcite" = 5, + "Dolomite" = 5 +) + +check_sign_cal_dol_dht <- function(old, new) { + if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) { + return(TRUE) + } + if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) { + return(TRUE) + } + return(FALSE) +} + +fuzz_input_dht_keys <- function(input) { + return(input[names(dht_species)]) +} + +check_sign_cal_dol_interp <- function(to_interp, data_set) { + data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE) + names(data_set) <- names(dht_species) + cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0) + dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0) + + cal_dol_same_sig <- cal == dol + return(rev(which(!cal_dol_same_sig))) +} + +check_neg_cal_dol <- function(result) { + neg_sign <- (result["Calcite"] <- 0) || (result["Dolomite"] < 0) + return(any(neg_sign)) +} + +hooks <- list( + dht_fill = check_sign_cal_dol_dht, + dht_fuzz = fuzz_input_dht_keys, + interp_pre_func = check_sign_cal_dol_interp, + interp_post_func = check_neg_cal_dol +) chemistry <- list( database = database, input_script = input_script, - dht_species = dht_species + dht_species = dht_species, + hooks = hooks ) ################################################################# diff --git a/bench/dolo/dolo_diffu_inner_large.R b/bench/dolo/dolo_diffu_inner_large.R index 4403ae46f..ce1da8246 100644 --- a/bench/dolo/dolo_diffu_inner_large.R +++ b/bench/dolo/dolo_diffu_inner_large.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2023-08-02 13:59:12 mluebke" +## Time-stamp: "Last modified 2023-08-16 17:05:04 mluebke" database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat") input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi") @@ -62,34 +62,34 @@ alpha_diffu <- c( ## list of boundary conditions/inner nodes vecinj_diffu <- list( - list( - "H" = 0.0001540445, - "O" = 0.002148006, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0, - "Cl" = 0.002, - "Mg" = 0.001 - ), - list( - "H" = 0.0001610193, - "O" = 0.002386934, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0.0, - "Cl" = 0.004, - "Mg" = 0.002 - ) + list( + "H" = 0.0001540445, + "O" = 0.002148006, + "Charge" = 1.90431e-16, + "C(4)" = 0, + "Ca" = 0, + "Cl" = 0.002, + "Mg" = 0.001 + ), + list( + "H" = 0.0001610193, + "O" = 0.002386934, + "Charge" = 1.90431e-16, + "C(4)" = 0, + "Ca" = 0.0, + "Cl" = 0.004, + "Mg" = 0.002 + ) ) vecinj_inner <- list( - l1 = c(1,400,200), - l2 = c(2,1400,800), - l3 = c(2,1600,800) + l1 = c(1, 400, 200), + l2 = c(2, 1400, 800), + l3 = c(2, 1600, 800) ) boundary <- list( -# "N" = c(1, rep(0, n-1)), + # "N" = c(1, rep(0, n-1)), "N" = rep(0, n), "E" = rep(0, m), "S" = rep(0, n), @@ -115,21 +115,59 @@ diffusion <- list( ################################################################# ## # Needed when using DHT -dht_species <- c("H" = 10, - "O" = 10, - "Charge" = 3, - "C(4)" = 5, - "Ca" = 5, - "Cl" = 5, - "Mg" = 5, - "Calcite" = 5, - "Dolomite" =5 - ) +dht_species <- c( + "H" = 10, + "O" = 10, + "Charge" = 3, + "C(4)" = 5, + "Ca" = 5, + "Cl" = 5, + "Mg" = 5, + "Calcite" = 5, + "Dolomite" = 5 +) + +check_sign_cal_dol_dht <- function(old, new) { + if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) { + return(TRUE) + } + if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) { + return(TRUE) + } + return(FALSE) +} + +fuzz_input_dht_keys <- function(input) { + return(input[names(dht_species)]) +} + +check_sign_cal_dol_interp <- function(to_interp, data_set) { + data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE) + names(data_set) <- names(dht_species) + cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0) + dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0) + + cal_dol_same_sig <- cal == dol + return(rev(which(!cal_dol_same_sig))) +} + +check_neg_cal_dol <- function(result) { + neg_sign <- (result["Calcite"] <- 0) || (result["Dolomite"] < 0) + return(any(neg_sign)) +} + +hooks <- list( + dht_fill = check_sign_cal_dol_dht, + dht_fuzz = fuzz_input_dht_keys, + interp_pre_func = check_sign_cal_dol_interp, + interp_post_func = check_neg_cal_dol +) chemistry <- list( database = database, input_script = input_script, - dht_species = dht_species + dht_species = dht_species, + hooks = hooks ) ################################################################# @@ -148,5 +186,5 @@ setup <- list( iterations = iterations, timesteps = rep(dt, iterations), store_result = TRUE, - out_save = seq(5, iterations, by=5) + out_save = seq(5, iterations, by = 5) ) diff --git a/bench/dolo/dolo_interp_long.R b/bench/dolo/dolo_interp_long.R index dedfe1c21..80770f07a 100644 --- a/bench/dolo/dolo_interp_long.R +++ b/bench/dolo/dolo_interp_long.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2023-08-02 13:47:06 mluebke" +## Time-stamp: "Last modified 2023-08-16 14:57:25 mluebke" database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat") input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi") @@ -62,25 +62,25 @@ alpha_diffu <- c( ## list of boundary conditions/inner nodes vecinj_diffu <- list( - list( - "H" = 1.110124E+02, - "O" = 5.550796E+01, - "Charge" = -3.230390327801E-08, - "C(4)" = 0, - "Ca" = 0, - "Cl" = 0.002, - "Mg" = 0.001 - ), - list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0.0, - "Cl" = 0.004, - "Mg" = 0.002 - ), - init_diffu + list( + "H" = 1.110124E+02, + "O" = 5.550796E+01, + "Charge" = -3.230390327801E-08, + "C(4)" = 0, + "Ca" = 0, + "Cl" = 0.002, + "Mg" = 0.001 + ), + list( + "H" = 110.683, + "O" = 55.3413, + "Charge" = 1.90431e-16, + "C(4)" = 0, + "Ca" = 0.0, + "Cl" = 0.004, + "Mg" = 0.002 + ), + init_diffu ) vecinj_inner <- list( @@ -132,19 +132,56 @@ dht_species <- c( ## # Optional when using Interpolation (example with less key species and custom ## # significant digits) -#pht_species <- c( +# pht_species <- c( # "C(4)" = 3, # "Ca" = 3, # "Mg" = 2, # "Calcite" = 2, # "Dolomite" = 2 -#) +# ) + +check_sign_cal_dol_dht <- function(old, new) { + if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) { + return(TRUE) + } + if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) { + return(TRUE) + } + return(FALSE) +} + +fuzz_input_dht_keys <- function(input) { + return(input[names(dht_species)]) +} + +check_sign_cal_dol_interp <- function(to_interp, data_set) { + data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE) + names(data_set) <- names(dht_species) + cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0) + dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0) + + cal_dol_same_sig <- cal == dol + return(rev(which(!cal_dol_same_sig))) +} + +check_neg_cal_dol <- function(result) { + neg_sign <- (result["Calcite"] <- 0) || (result["Dolomite"] < 0) + return(any(neg_sign)) +} + +hooks <- list( + dht_fill = check_sign_cal_dol_dht, + dht_fuzz = fuzz_input_dht_keys, + interp_pre_func = check_sign_cal_dol_interp, + interp_post_func = check_neg_cal_dol +) chemistry <- list( database = database, input_script = input_script, - dht_species = dht_species -# pht_species = pht_species + dht_species = dht_species, + hooks = hooks + # pht_species = pht_species ) ################################################################# diff --git a/include/poet/ChemistryModule.hpp b/include/poet/ChemistryModule.hpp index f5930184d..c0c7f9083 100644 --- a/include/poet/ChemistryModule.hpp +++ b/include/poet/ChemistryModule.hpp @@ -1,11 +1,10 @@ -// Time-stamp: "Last modified 2023-08-08 13:15:49 mluebke" +// Time-stamp: "Last modified 2023-08-15 14:36:28 mluebke" #ifndef CHEMISTRYMODULE_H_ #define CHEMISTRYMODULE_H_ #include "DHT_Wrapper.hpp" #include "DataStructures.hpp" -#include "Field.hpp" #include "Interpolation.hpp" #include "IrmResult.h" #include "PhreeqcRM.h" @@ -382,7 +381,7 @@ protected: uint32_t prop_count = 0; std::vector prop_names; - Field chem_field{0}; + Field chem_field; static constexpr int MODULE_COUNT = 5; diff --git a/include/poet/DHT_Wrapper.hpp b/include/poet/DHT_Wrapper.hpp index a066d3ccd..a74db4f2d 100644 --- a/include/poet/DHT_Wrapper.hpp +++ b/include/poet/DHT_Wrapper.hpp @@ -1,4 +1,4 @@ -// Time-stamp: "Last modified 2023-08-08 18:01:12 mluebke" +// Time-stamp: "Last modified 2023-08-15 14:57:51 mluebke" /* ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of @@ -25,6 +25,8 @@ #include "DataStructures.hpp" #include "LookupKey.hpp" +#include "RInsidePOET.hpp" +#include "SimParams.hpp" #include "enums.hpp" #include "poet/HashFunctions.hpp" #include "poet/LookupKey.hpp" @@ -85,6 +87,8 @@ public: DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size, const NamedVector &key_species, const std::vector &key_indices, + const std::vector &output_names, + const ChemistryParams::Chem_Hook_Functions &hooks, uint32_t data_count); /** * @brief Destroy the dht wrapper object @@ -213,7 +217,11 @@ public: double dt) { dht_results.keys.resize(input_values.size()); for (std::size_t i = 0; i < input_values.size(); i++) { - dht_results.keys[i] = fuzzForDHT(input_values[i], dt); + if (this->hooks.dht_fuzz.isValid()) { + dht_results.keys[i] = fuzzForDHT_R(input_values[i], dt); + } else { + dht_results.keys[i] = fuzzForDHT(input_values[i], dt); + } } } @@ -225,6 +233,7 @@ private: MPI_Comm communicator; LookupKey fuzzForDHT(const std::vector &cell, double dt); + LookupKey fuzzForDHT_R(const std::vector &cell, double dt); std::vector outputToInputAndRates(const std::vector &old_results, @@ -242,6 +251,10 @@ private: std::vector dht_prop_type_vector; std::vector input_key_elements; + const std::vector &output_names; + + const ChemistryParams::Chem_Hook_Functions &hooks; + DHT_ResultObject dht_results; std::array base_totals{0}; diff --git a/include/poet/DataStructures.hpp b/include/poet/DataStructures.hpp index c2ab99d77..834ceca50 100644 --- a/include/poet/DataStructures.hpp +++ b/include/poet/DataStructures.hpp @@ -3,12 +3,16 @@ #include "enums.hpp" +#include + +#include +#include #include -#include +#include +#include #include #include #include - namespace poet { struct WorkPackage { @@ -24,28 +28,193 @@ struct WorkPackage { } }; -template class NamedVector { +template class NamedVector : public Rcpp::NumericVector { public: - void insert(std::pair to_insert) { - this->names.push_back(to_insert.first); - this->values.push_back(to_insert.second); - container_size++; + NamedVector() : Rcpp::NumericVector(){}; + + NamedVector(const std::vector &in_names, + const std::vector &in_values) + : Rcpp::NumericVector(Rcpp::wrap(in_values)) { + this->names() = Rcpp::CharacterVector(Rcpp::wrap(in_names)); } - bool empty() const { return this->container_size == 0; } + NamedVector(const SEXP &s) : Rcpp::NumericVector(s){}; - std::size_t size() const { return this->container_size; } + bool empty() const { return (this->size() == 0); } - value_type &operator[](std::size_t i) { return values[i]; } + std::vector getNames() const { + return Rcpp::as>(this->names()); + } + std::vector getValues() const { return Rcpp::as>(*this); } +}; - const std::vector &getNames() const { return this->names; } - const std::vector &getValues() const { return this->values; } +using FieldColumn = std::vector; + +/** + * Stores data for input/output of a module. The class keeps track of all + * defined properties with name and import/export to 1D and 2D STL vectors. + * Also, it eases the update process of a field with an input from another + * field. + * + * It can be seen as an R-like data frame, but with less access to the members. + * Species values are stored as "columns", where column is a STL vector. + */ +class Field : private std::unordered_map { +public: + Field(){}; + + /** + * Creates a new instance of a field with fixed expected vector size. + * + * \param vec_s expected vector size of each component/column. + */ + Field(std::uint32_t vec_s) : req_vec_size(vec_s){}; + + /** + * Initializes instance with a 2D vector and according names for each columnn. + * There is no check if names were given multiple times. The order of name + * vector also defines the ordering of the output. + * + * \param vec_s expected vector size of each component/column. + * \param input 2D vector using STL semantic describing the current state of + * the field. + * \param prop_vec Name of each vector in the input. Shall match + * the count of vectors. + * + * \exception std::runtime_error If prop_vec size doesn't match input vector + * size or column vectors size doesn't match expected vector size. + */ + Field(std::uint32_t vec_s, const std::vector> &input, + const std::vector &prop_vec, bool row_major = false); + + /** + * Initializes instance with a 1D continious memory vector and according names + * for each columnn. There is no check if names were given multiple times. The + * order of name vector also defines the ordering of the output. + * + * \param vec_s expected vector size of each component/column. + * \param input 1D vector using STL semantic describing the current state of + * the field storing each column starting at index *i times requested vector + * size*. + * \param prop_vec Name of each vector in the input. Shall match the + * count of vectors. + * + * \exception std::runtime_error If prop_vec size doesn't match input vector + * size or column vectors size doesn't match expected vector size. + */ + Field(std::uint32_t vec_s, const std::vector &input, + const std::vector &prop_vec); + + Field(const SEXP &s_init) { fromSEXP(s_init); } + + Field &operator=(const Field &rhs) { + this->req_vec_size = rhs.req_vec_size; + this->props = rhs.props; + + this->clear(); + + for (const auto &pair : rhs) { + this->insert(pair); + } + return *this; + } + + /** + * Read in a (previously exported) 1D vector. It has to have the same + * dimensions as the current column count times the requested vector size of + * this instance. Import order is given by the species name vector. + * + * \param cont_field 1D field as vector. + * + * \exception std::runtime_error Input vector does not match the expected + * size. + */ + Field &operator=(const FieldColumn &cont_field); + + /** + * Read in a (previously exported) 2D vector. It has to have the same + * dimensions as the current column count of this instance and each vector + * must have the size of the requested vector size. Import order is given by + * the species name vector. + * + * \param cont_field 2D field as vector. + * + * \exception std::runtime_error Input vector has more or less elements than + * the instance or a column vector does not match expected vector size. + */ + Field &operator=(const std::vector &cont_field); + + Field &operator=(const SEXP &s_rhs) { + fromSEXP(s_rhs); + return *this; + } + + /** + * Returns a reference to the column vector with given name. Creates a new + * vector if prop was not found. The prop name will be added to the end of the + * list. + * + * \param key Name of the prop. + * + * \return Reference to the column vector. + */ + FieldColumn &operator[](const std::string &key); + + /** + * Returns the names of all species defined in the instance. + * + * \return Vector containing all species names in output order. + */ + auto GetProps() const -> std::vector { return this->props; } + + /** + * Return the requested vector size. + * + * \return Requested vector size set in the instanciation of the object. + */ + auto GetRequestedVecSize() const -> std::uint32_t { + return this->req_vec_size; + }; + + /** + * Updates all species with values from another field. If one element of the + * input field doesn't match the names of the calling instance it will get + * skipped. + * + * \param input Field to update the current instance's columns. + */ + void update(const Field &input); + + /** + * Builds a new 1D vector with the current state of the instance. The output + * order is given by the given species name vector set earlier and/or added + * values using the [] operator. + * + * \return 1D STL vector stored each column one after another. + */ + auto AsVector() const -> FieldColumn; + + /** + * Builds a new 2D vector with the current state of the instance. The output + * order is given by the given species name vector set earlier and/or added + * values using the [] operator. + * + * \return 2D STL vector stored each column one after anothe in a new vector + * element. + */ + auto As2DVector() const -> std::vector; + + SEXP asSEXP() const; + + std::size_t cols() const { return this->props.size(); } private: - std::size_t container_size{0}; + std::uint32_t req_vec_size{0}; - std::vector names{}; - std::vector values{}; + std::vector props{{}}; + + void fromSEXP(const SEXP &s_rhs); }; + } // namespace poet #endif // DATASTRUCTURES_H_ diff --git a/include/poet/DiffusionModule.hpp b/include/poet/DiffusionModule.hpp index e9b9c61f5..8ea379027 100644 --- a/include/poet/DiffusionModule.hpp +++ b/include/poet/DiffusionModule.hpp @@ -21,7 +21,7 @@ #ifndef DIFFUSION_MODULE_H #define DIFFUSION_MODULE_H -#include "Field.hpp" +#include "DataStructures.hpp" #include "SimParams.hpp" #include "poet/SimParams.hpp" #include @@ -94,7 +94,8 @@ private: enum { DIM_1D = 1, DIM_2D }; - void initialize(const poet::DiffusionParams &args); + void initialize(const poet::DiffusionParams &args, + std::uint32_t n_grid_cells); uint8_t dim; diff --git a/include/poet/Field.hpp b/include/poet/Field.hpp deleted file mode 100644 index cb4e46996..000000000 --- a/include/poet/Field.hpp +++ /dev/null @@ -1,159 +0,0 @@ -#ifndef FIELD_H_ -#define FIELD_H_ - -#include -#include -#include -#include -#include -#include - -namespace poet { - -using FieldColumn = std::vector; - -/** - * Stores data for input/output of a module. The class keeps track of all - * defined properties with name and import/export to 1D and 2D STL vectors. - * Also, it eases the update process of a field with an input from another - * field. - * - * It can be seen as an R-like data frame, but with less access to the members. - * Species values are stored as "columns", where column is a STL vector. - */ -class Field : private std::unordered_map { -public: - /** - * Creates a new instance of a field with fixed expected vector size. - * - * \param vec_s expected vector size of each component/column. \ - */ - Field(uint32_t vec_s) : req_vec_size(vec_s){}; - - /** - * Initializes instance with a 2D vector and according names for each columnn. - * There is no check if names were given multiple times. The order of name - * vector also defines the ordering of the output. - * - * \param input 2D vector using STL semantic describing the current state of - * the field. - * \param prop_vec Name of each vector in the input. Shall match - * the count of vectors. - * - * \exception std::runtime_error If prop_vec size doesn't match input vector - * size or column vectors size doesn't match expected vector size. - */ - void InitFromVec(const std::vector> &input, - const std::vector &prop_vec); - - /** - * Initializes instance with a 1D continious memory vector and according names - * for each columnn. There is no check if names were given multiple times. The - * order of name vector also defines the ordering of the output. - * - * \param input 1D vector using STL semantic describing the current state of - * the field storing each column starting at index *i times requested vector - * size*. - * \param prop_vec Name of each vector in the input. Shall match the - * count of vectors. - * - * \exception std::runtime_error If prop_vec size doesn't match input vector - * size or column vectors size doesn't match expected vector size. - */ - void InitFromVec(const std::vector &input, - const std::vector &prop_vec); - - Field &operator=(const Field &rhs) { - this->req_vec_size = rhs.req_vec_size; - this->props = rhs.props; - return *this; - } - - /** - * Returns a reference to the column vector with given name. Creates a new - * vector if prop was not found. The prop name will be added to the end of the - * list. - * - * \param key Name of the prop. - * - * \return Reference to the column vector. - */ - FieldColumn &operator[](const std::string &key); - - /** - * Returns the names of all species defined in the instance. - * - * \return Vector containing all species names in output order. - */ - auto GetProps() const -> std::vector { return this->props; } - - /** - * Return the requested vector size. - * - * \return Requested vector size set in the instanciation of the object. - */ - auto GetRequestedVecSize() const -> std::uint32_t { - return this->req_vec_size; - }; - - /** - * Updates all species with values from another field. If one element of the - * input field doesn't match the names of the calling instance it will get - * skipped. - * - * \param input Field to update the current instance's columns. - */ - void UpdateFromField(const Field &input); - - /** - * Builds a new 1D vector with the current state of the instance. The output - * order is given by the given species name vector set earlier and/or added - * values using the [] operator. - * - * \return 1D STL vector stored each column one after another. - */ - auto AsVector() const -> FieldColumn; - - /** - * Builds a new 2D vector with the current state of the instance. The output - * order is given by the given species name vector set earlier and/or added - * values using the [] operator. - * - * \return 2D STL vector stored each column one after anothe in a new vector - * element. - */ - auto As2DVector() const -> std::vector; - - /** - * Read in a (previously exported) 1D vector. It has to have the same - * dimensions as the current column count times the requested vector size of - * this instance. Import order is given by the species name vector. - * - * \param cont_field 1D field as vector. - * - * \exception std::runtime_error Input vector does not match the expected - * size. - */ - void SetFromVector(const FieldColumn &cont_field); - - /** - * Read in a (previously exported) 2D vector. It has to have the same - * dimensions as the current column count of this instance and each vector - * must have the size of the requested vector size. Import order is given by - * the species name vector. - * - * \param cont_field 2D field as vector. - * - * \exception std::runtime_error Input vector has more or less elements than - * the instance or a column vector does not match expected vector size. - */ - void SetFromVector(const std::vector &cont_field); - -private: - std::uint32_t req_vec_size; - - std::vector props; -}; - -} // namespace poet -#endif // FIELD_H_ diff --git a/include/poet/Interpolation.hpp b/include/poet/Interpolation.hpp index 5f85872ca..71b9f9a01 100644 --- a/include/poet/Interpolation.hpp +++ b/include/poet/Interpolation.hpp @@ -1,4 +1,4 @@ -// Time-stamp: "Last modified 2023-08-09 12:52:37 mluebke" +// Time-stamp: "Last modified 2023-08-16 16:49:31 mluebke" #ifndef INTERPOLATION_H_ #define INTERPOLATION_H_ @@ -9,6 +9,7 @@ #include "LookupKey.hpp" #include "poet/DHT_Wrapper.hpp" #include "poet/Rounding.hpp" +#include "poet/SimParams.hpp" #include #include #include @@ -46,7 +47,7 @@ public: std::vector> in_values; std::vector> out_values; - std::uint32_t getSize() const { + std::uint32_t getMemSize() const { std::uint32_t sum{0}; for (const auto &results : out_values) { sum += results.size() * sizeof(double); @@ -65,9 +66,9 @@ public: void writeLocationToPHT(LookupKey key, DHT_Location location); const PHT_Result &query(const LookupKey &key, - const std::uint32_t min_entries_needed, - const std::uint32_t input_count, - const std::uint32_t output_count); + std::uint32_t min_entries_needed, + std::uint32_t input_count, + std::uint32_t output_count); std::uint64_t getLocations(const LookupKey &key); @@ -163,7 +164,9 @@ public: std::uint64_t size_per_process, std::uint32_t min_entries_needed, DHT_Wrapper &dht, const NamedVector &interp_key_signifs, - const std::vector &dht_key_indices); + const std::vector &dht_key_indices, + const std::vector &out_names, + const ChemistryParams::Chem_Hook_Functions &hooks); enum result_status { RES_OK, INSUFFICIENT_DATA, NOT_NEEDED }; @@ -258,6 +261,10 @@ private: return out_key; } + + const ChemistryParams::Chem_Hook_Functions &hooks; + const std::vector &out_names; + const std::vector dht_names; }; } // namespace poet diff --git a/include/poet/LookupKey.hpp b/include/poet/LookupKey.hpp index c3366735f..fbb064c53 100644 --- a/include/poet/LookupKey.hpp +++ b/include/poet/LookupKey.hpp @@ -1,4 +1,4 @@ -// Time-stamp: "Last modified 2023-07-26 10:20:10 mluebke" +// Time-stamp: "Last modified 2023-08-11 10:12:52 mluebke" #ifndef LOOKUPKEY_H_ #define LOOKUPKEY_H_ @@ -27,8 +27,6 @@ union Lookup_Keyelement { class LookupKey : public std::vector { public: - explicit LookupKey(size_type count); - using std::vector::vector; std::vector to_double() const; diff --git a/include/poet/RInsidePOET.hpp b/include/poet/RInsidePOET.hpp index fb2102685..6467e5760 100644 --- a/include/poet/RInsidePOET.hpp +++ b/include/poet/RInsidePOET.hpp @@ -1,11 +1,11 @@ #ifndef RPOET_H_ #define RPOET_H_ -#include "DataStructures.hpp" - #include #include #include +#include +#include #include #include #include @@ -28,34 +28,32 @@ public: this->parseEval("'" + R_name + "' %in% names(" + where + ")")); } - template inline T getSTL(const std::string &R_name) { - return Rcpp::as(RInside::operator[](R_name)); - } - - template - poet::NamedVector wrapNamedVector(const std::string &R_name) { - std::size_t name_size = this->parseEval("length(names(" + R_name + "))"); - if (name_size == 0) { - throw std::runtime_error("wrapNamedVector: " + R_name + - " is not a named vector!"); - } - - auto names = Rcpp::as>( - this->parseEval("names(" + R_name + ")")); - auto values = Rcpp::as(this->parseEval(R_name)); - - poet::NamedVector ret_map; - - for (std::size_t i = 0; i < names.size(); i++) { - ret_map.insert( - std::make_pair(names[i], static_cast(values[i]))); - } - - return ret_map; - } - private: RInsidePOET() : RInside(){}; }; +template class RHookFunction { +public: + RHookFunction() {} + RHookFunction(RInside &R, const std::string &f_name) { + try { + this->func = Rcpp::Function(Rcpp::as(R.parseEval(f_name.c_str()))); + } catch (const std::exception &e) { + } + } + + template T operator()(Args... args) const { + if (func.has_value()) { + return (Rcpp::as(this->func.value()(args...))); + } else { + throw std::exception(); + } + } + + bool isValid() const { return this->func.has_value(); } + +private: + std::optional func; +}; + #endif // RPOET_H_ diff --git a/include/poet/SimParams.hpp b/include/poet/SimParams.hpp index a73669865..9928bcdba 100644 --- a/include/poet/SimParams.hpp +++ b/include/poet/SimParams.hpp @@ -22,13 +22,14 @@ #define PARSER_H #include +#include #include #include #include #include -#include "Macros.hpp" #include "DataStructures.hpp" +#include "Macros.hpp" #include "RInsidePOET.hpp" #include "argh.hpp" // Argument handler https://github.com/adishavit/argh #include @@ -124,6 +125,13 @@ struct ChemistryParams { std::uint32_t interp_min_entries; NamedVector pht_signifs; + struct Chem_Hook_Functions { + RHookFunction dht_fill; + RHookFunction> dht_fuzz; + RHookFunction> interp_pre; + RHookFunction interp_post; + } hooks; + void initFromR(RInsidePOET &R); }; @@ -240,8 +248,8 @@ public: private: std::list validateOptions(argh::parser cmdl); - const std::set flaglist{ - "ignore-result", "dht", "P", "progress", "interp"}; + const std::set flaglist{"ignore-result", "dht", "P", "progress", + "interp"}; const std::set paramlist{ "work-package-size", "dht-strategy", "dht-size", "dht-snaps", diff --git a/src/ChemistryModule/ChemistryModule.cpp b/src/ChemistryModule/ChemistryModule.cpp index 6f194e8bc..3060e8617 100644 --- a/src/ChemistryModule/ChemistryModule.cpp +++ b/src/ChemistryModule/ChemistryModule.cpp @@ -291,7 +291,6 @@ void poet::ChemistryModule::initializeField(const Field &trans_field) { if (is_master) { this->n_cells = trans_field.GetRequestedVecSize(); - chem_field = Field(n_cells); std::vector> phreeqc_dump(this->nxyz); this->getDumpedField(phreeqc_dump); @@ -305,7 +304,7 @@ void poet::ChemistryModule::initializeField(const Field &trans_field) { const auto tmp_buffer{init_vec}; this->unshuffleField(tmp_buffer, n_cells, prop_count, 1, init_vec); - chem_field.InitFromVec(init_vec, prop_names); + this->chem_field = Field(n_cells, init_vec, prop_names); return; } @@ -325,7 +324,7 @@ void poet::ChemistryModule::initializeField(const Field &trans_field) { initial_values.push_back(std::move(init)); } - chem_field.InitFromVec(initial_values, prop_names); + this->chem_field = Field(n_cells, initial_values, prop_names); } else { ChemBCast(base_totals.data(), 2, MPI_DOUBLE); } @@ -366,11 +365,9 @@ void poet::ChemistryModule::initializeDHT( auto map_copy = key_species; if (key_species.empty()) { - const auto &all_species = this->prop_names; - for (std::size_t i = 0; i < all_species.size(); i++) { - map_copy.insert(std::make_pair(all_species[i], - DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT)); - } + std::vector default_signif( + this->prop_names.size(), DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT); + map_copy = NamedVector(this->prop_names, default_signif); } auto key_indices = parseDHTSpeciesVec(key_species, this->prop_names); @@ -381,8 +378,9 @@ void poet::ChemistryModule::initializeDHT( const std::uint64_t dht_size = size_mb * MB_FACTOR; - this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices, - this->prop_count); + this->dht = + new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices, + this->prop_names, params.hooks, this->prop_count); this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1)); } } @@ -488,8 +486,8 @@ void poet::ChemistryModule::initializeInterp( const uint64_t pht_size = size_mb * MB_FACTOR; interp = std::make_unique( - bucket_size, pht_size, min_entries, *(this->dht), map_copy, - key_indices); + bucket_size, pht_size, min_entries, *(this->dht), map_copy, key_indices, + this->prop_names, this->params.hooks); interp->setInterpolationFunction(inverseDistanceWeighting); } diff --git a/src/ChemistryModule/MasterFunctions.cpp b/src/ChemistryModule/MasterFunctions.cpp index 77bfa6023..51b656da3 100644 --- a/src/ChemistryModule/MasterFunctions.cpp +++ b/src/ChemistryModule/MasterFunctions.cpp @@ -342,7 +342,7 @@ void poet::ChemistryModule::MasterRunSequential() { std::vector out_vec{shuffled_field}; unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec); - chem_field.SetFromVector(out_vec); + chem_field = out_vec; } void poet::ChemistryModule::MasterRunParallel() { @@ -423,7 +423,7 @@ void poet::ChemistryModule::MasterRunParallel() { std::vector out_vec{mpi_buffer}; unshuffleField(mpi_buffer, this->n_cells, this->prop_count, wp_sizes_vector.size(), out_vec); - chem_field.SetFromVector(out_vec); + chem_field = out_vec; /* do master stuff */ diff --git a/src/ChemistryModule/SurrogateModels/DHT_Wrapper.cpp b/src/ChemistryModule/SurrogateModels/DHT_Wrapper.cpp index bc3e8cfd5..d1ee112c2 100644 --- a/src/ChemistryModule/SurrogateModels/DHT_Wrapper.cpp +++ b/src/ChemistryModule/SurrogateModels/DHT_Wrapper.cpp @@ -1,4 +1,4 @@ -// Time-stamp: "Last modified 2023-08-09 14:05:01 mluebke" +// Time-stamp: "Last modified 2023-08-16 16:44:17 mluebke" /* ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of @@ -27,7 +27,10 @@ #include "poet/Rounding.hpp" #include "poet/enums.hpp" +#include "poet/RInsidePOET.hpp" + #include +#include #include #include #include @@ -43,10 +46,12 @@ namespace poet { DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size, const NamedVector &key_species, const std::vector &key_indices, + const std::vector &_output_names, + const ChemistryParams::Chem_Hook_Functions &_hooks, uint32_t data_count) : key_count(key_indices.size()), data_count(data_count), input_key_elements(key_indices), communicator(dht_comm), - key_species(key_species) { + key_species(key_species), output_names(_output_names), hooks(_hooks) { // initialize DHT object // key size = count of key elements + timestep uint32_t key_size = (key_count + 1) * sizeof(Lookup_Keyelement); @@ -128,12 +133,13 @@ void DHT_Wrapper::fillDHT(const WorkPackage &work_package) { // check if calcite or dolomite is absent and present, resp.n and vice // versa in input/output. If this is the case -> Do not write to DHT! // HACK: hardcoded, should be fixed! - if ((work_package.input[i][7] == 0) != (work_package.output[i][7] == 0)) { - continue; - } + if (hooks.dht_fill.isValid()) { + NamedVector old_values(output_names, work_package.input[i]); + NamedVector new_values(output_names, work_package.output[i]); - if ((work_package.input[i][9] == 0) != (work_package.output[i][9] == 0)) { - continue; + if (hooks.dht_fill(old_values, new_values)) { + continue; + } } uint32_t proc, index; @@ -228,6 +234,38 @@ void DHT_Wrapper::printStatistics() { } } +LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector &cell, + double dt) { + const auto c_zero_val = std::pow(10, AQUEOUS_EXP); + + NamedVector input_nv(this->output_names, cell); + + const std::vector eval_vec = hooks.dht_fuzz(input_nv); + assert(eval_vec.size() == this->key_count); + LookupKey vecFuzz(this->key_count + 1, {.0}); + + DHT_Rounder rounder; + + int totals_i = 0; + // introduce fuzzing to allow more hits in DHT + // loop over every variable of grid cell + for (std::uint32_t i = 0; i < eval_vec.size(); i++) { + double curr_key = eval_vec[i]; + if (curr_key != 0) { + if (this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL) { + curr_key -= base_totals[totals_i++]; + } + vecFuzz[i] = + rounder.round(curr_key, dht_signif_vector[i], + this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL); + } + } + // add timestep to the end of the key as double value + vecFuzz[this->key_count].fp_element = dt; + + return vecFuzz; +} + LookupKey DHT_Wrapper::fuzzForDHT(const std::vector &cell, double dt) { const auto c_zero_val = std::pow(10, AQUEOUS_EXP); diff --git a/src/ChemistryModule/SurrogateModels/InterpolationModule.cpp b/src/ChemistryModule/SurrogateModels/InterpolationModule.cpp index 1de57554b..9384731b2 100644 --- a/src/ChemistryModule/SurrogateModels/InterpolationModule.cpp +++ b/src/ChemistryModule/SurrogateModels/InterpolationModule.cpp @@ -1,10 +1,13 @@ -// Time-stamp: "Last modified 2023-08-11 11:01:11 delucia" +// Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke" #include "poet/DHT_Wrapper.hpp" +#include "poet/DataStructures.hpp" #include "poet/HashFunctions.hpp" #include "poet/Interpolation.hpp" #include "poet/LookupKey.hpp" #include "poet/Rounding.hpp" +#include +#include #include #include #include @@ -27,9 +30,13 @@ InterpolationModule::InterpolationModule( std::uint32_t entries_per_bucket, std::uint64_t size_per_process, std::uint32_t min_entries_needed, DHT_Wrapper &dht, const NamedVector &interp_key_signifs, - const std::vector &dht_key_indices) + const std::vector &dht_key_indices, + const std::vector &_out_names, + const ChemistryParams::Chem_Hook_Functions &_hooks) : dht_instance(dht), key_signifs(interp_key_signifs), - key_indices(dht_key_indices), min_entries_needed(min_entries_needed) { + key_indices(dht_key_indices), min_entries_needed(min_entries_needed), + dht_names(dht.getKeySpecies().getNames()), out_names(_out_names), + hooks(_hooks) { initPHT(this->key_signifs.size(), entries_per_bucket, size_per_process, dht.getCommunicator()); @@ -63,45 +70,41 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) { const auto dht_results = this->dht_instance.getDHTResults(); - for (int i = 0; i < work_package.size; i++) { - if (work_package.mapping[i] != CHEM_PQC) { - interp_result.status[i] = NOT_NEEDED; + for (int wp_i = 0; wp_i < work_package.size; wp_i++) { + if (work_package.mapping[wp_i] != CHEM_PQC) { + interp_result.status[wp_i] = NOT_NEEDED; continue; } - const auto rounded_key = roundKey(dht_results.keys[i]); + const auto rounded_key = roundKey(dht_results.keys[wp_i]); auto pht_result = pht->query(rounded_key, this->min_entries_needed, dht_instance.getInputCount(), dht_instance.getOutputCount()); - int pht_i = 0; + if (pht_result.size < this->min_entries_needed) { + interp_result.status[wp_i] = INSUFFICIENT_DATA; + continue; + } + + if (hooks.interp_pre.isValid()) { + NamedVector nv_in(this->out_names, work_package.input[wp_i]); + + auto rm_indices = hooks.interp_pre(nv_in, pht_result.in_values); + + pht_result.size -= rm_indices.size(); - while (pht_i < pht_result.size) { if (pht_result.size < this->min_entries_needed) { - break; - } - - auto in_it = pht_result.in_values.begin() + pht_i; - auto out_it = pht_result.out_values.begin() + pht_i; - - bool same_sig_calcite = (pht_result.in_values[pht_i][7] == 0) == - (work_package.input[i][7] == 0); - bool same_sig_dolomite = (pht_result.in_values[pht_i][8] == 0) == - (work_package.input[i][9] == 0); - if (!same_sig_calcite || !same_sig_dolomite) { - pht_result.size -= 1; - pht_result.in_values.erase(in_it); - pht_result.out_values.erase(out_it); + interp_result.status[wp_i] = INSUFFICIENT_DATA; continue; } - pht_i += 1; - } - - if (pht_result.size < this->min_entries_needed) { - interp_result.status[i] = INSUFFICIENT_DATA; - continue; + for (const auto &index : rm_indices) { + pht_result.in_values.erase( + std::next(pht_result.in_values.begin(), index - 1)); + pht_result.out_values.erase( + std::next(pht_result.out_values.begin(), index - 1)); + } } #ifdef POET_PHT_ADD @@ -109,20 +112,17 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) { #endif double start_fc = MPI_Wtime(); - // mean water - // double mean_water = 0; - // for (int out_i = 0; out_i < pht_result.size; out_i++) { - // mean_water += pht_result.out_values[out_i][0]; - // } - // mean_water /= pht_result.size; - work_package.output[i] = - f_interpolate(dht_instance.getKeyElements(), work_package.input[i], + work_package.output[wp_i] = + f_interpolate(dht_instance.getKeyElements(), work_package.input[wp_i], pht_result.in_values, pht_result.out_values); - if (work_package.output[i][7] < 0 || work_package.output[i][9] < 0) { - interp_result.status[i] = INSUFFICIENT_DATA; - continue; + if (hooks.interp_post.isValid()) { + NamedVector nv_result(this->out_names, work_package.output[wp_i]); + if (hooks.interp_post(nv_result)) { + interp_result.status[wp_i] = INSUFFICIENT_DATA; + continue; + } } // interp_result.results[i][0] = mean_water; @@ -130,8 +130,8 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) { this->interpolations++; - work_package.mapping[i] = CHEM_INTERP; - interp_result.status[i] = RES_OK; + work_package.mapping[wp_i] = CHEM_INTERP; + interp_result.status[wp_i] = RES_OK; } } diff --git a/src/ChemistryModule/SurrogateModels/ProximityHashTable.cpp b/src/ChemistryModule/SurrogateModels/ProximityHashTable.cpp index 285ed21d3..3ea44430d 100644 --- a/src/ChemistryModule/SurrogateModels/ProximityHashTable.cpp +++ b/src/ChemistryModule/SurrogateModels/ProximityHashTable.cpp @@ -1,4 +1,4 @@ -// Time-stamp: "Last modified 2023-08-09 13:32:11 mluebke" +// Time-stamp: "Last modified 2023-08-15 14:50:59 mluebke" #include "poet/DHT_Wrapper.hpp" #include "poet/HashFunctions.hpp" @@ -203,7 +203,7 @@ void ProximityHashTable::Cache::operator()(const LookupKey &key, if (this->free_mem < 0) { const LookupKey &to_del = this->lru_queue.back(); const auto elem_d = this->find(to_del); - this->free_mem += elem_d->second.getSize(); + this->free_mem += elem_d->second.getMemSize(); this->erase(to_del); this->keyfinder.erase(to_del); this->lru_queue.pop_back(); @@ -212,7 +212,7 @@ void ProximityHashTable::Cache::operator()(const LookupKey &key, this->insert({key, val}); this->lru_queue.emplace_front(key); this->keyfinder[key] = lru_queue.begin(); - this->free_mem -= val.getSize(); + this->free_mem -= val.getMemSize(); return; } diff --git a/src/ChemistryModule/WorkerFunctions.cpp b/src/ChemistryModule/WorkerFunctions.cpp index e13641941..624a78b4f 100644 --- a/src/ChemistryModule/WorkerFunctions.cpp +++ b/src/ChemistryModule/WorkerFunctions.cpp @@ -1,4 +1,4 @@ -// Time-stamp: "Last modified 2023-08-15 17:33:00 mluebke" +// Time-stamp: "Last modified 2023-08-16 14:50:04 mluebke" #include "poet/ChemistryModule.hpp" #include "poet/DHT_Wrapper.hpp" @@ -52,7 +52,7 @@ void poet::ChemistryModule::WorkerLoop() { break; } case CHEM_INIT_SPECIES: { - Field dummy{0}; + Field dummy; initializeField(dummy); break; } @@ -220,8 +220,8 @@ void poet::ChemistryModule::WorkerPostIter(MPI_Status &prope_status, if (this->dht_enabled) { dht_hits.push_back(dht->getHits()); dht_evictions.push_back(dht->getEvictions()); - dht->resetCounter(); - + dht->resetCounter(); + if (this->dht_snaps_type == DHT_SNAPS_ITEREND) { WorkerWriteDHTDump(iteration); } @@ -238,8 +238,9 @@ void poet::ChemistryModule::WorkerPostIter(MPI_Status &prope_status, interp->dumpPHTState(out.str()); } } -} + RInsidePOET::getInstance().parseEvalQ("gc()"); +} void poet::ChemistryModule::WorkerPostSim(uint32_t iteration) { if (this->dht_enabled && this->dht_snaps_type >= DHT_SNAPS_ITEREND) { @@ -248,7 +249,7 @@ void poet::ChemistryModule::WorkerPostSim(uint32_t iteration) { if (this->interp_enabled && this->dht_snaps_type >= DHT_SNAPS_ITEREND) { std::stringstream out; out << this->dht_file_out_dir << "/iter_" << std::setfill('0') - << std::setw(this->file_pad) << iteration << ".pht"; + << std::setw(this->file_pad) << iteration << ".pht"; interp->dumpPHTState(out.str()); } } @@ -265,7 +266,7 @@ void poet::ChemistryModule::WorkerWriteDHTDump(uint32_t iteration) { std::cout << "CPP: Worker: Successfully written DHT to file " << out.str() << "\n"; } - + void poet::ChemistryModule::WorkerReadDHTDump( const std::string &dht_input_file) { int res = dht->fileToTable((char *)dht_input_file.c_str()); diff --git a/src/DataStructures/Field.cpp b/src/DataStructures/Field.cpp index 13b143c1e..de364a599 100644 --- a/src/DataStructures/Field.cpp +++ b/src/DataStructures/Field.cpp @@ -1,31 +1,53 @@ -#include "poet/Field.hpp" +#include "poet/DataStructures.hpp" +#include +#include +#include #include #include #include #include #include -void poet::Field::InitFromVec(const std::vector> &input, - const std::vector &prop_vec) { - if (prop_vec.size() != input.size()) { - throw std::runtime_error("Prop vector shall name all elements."); - } - - auto name_it = prop_vec.begin(); - - for (const auto &in_vec : input) { - if (in_vec.size() != req_vec_size) { - throw std::runtime_error( - "Input vector doesn't match expected vector size."); +poet::Field::Field(std::uint32_t vec_s, + const std::vector> &input, + const std::vector &prop_vec, bool row_major) + : req_vec_size(vec_s), props(prop_vec) { + if (row_major) { + if (this->props.size() != input[0].size()) { + throw std::runtime_error("Prop vector shall name all elements."); } - this->insert({*(name_it++), in_vec}); - } - this->props = prop_vec; + const std::size_t n_input = input.size(); + + for (std::size_t i = 0; i < this->props.size(); i++) { + std::vector curr_col(n_input); + + for (std::size_t j = 0; j < n_input; j++) { + curr_col[j] = input[j][i]; + } + + this->insert({this->props[i], std::move(curr_col)}); + } + } else { + if (this->props.size() != input.size()) { + throw std::runtime_error("Prop vector shall name all elements."); + } + + auto name_it = this->props.begin(); + + for (const auto &in_vec : input) { + if (in_vec.size() != req_vec_size) { + throw std::runtime_error( + "Input vector doesn't match expected vector size."); + } + this->insert({*(name_it++), in_vec}); + } + } } -void poet::Field::InitFromVec(const std::vector &input, - const std::vector &prop_vec) { +poet::Field::Field(std::uint32_t vec_s, const std::vector &input, + const std::vector &prop_vec) + : Field(vec_s) { const uint32_t expected_size = prop_vec.size() * req_vec_size; if (expected_size != input.size()) { @@ -62,47 +84,7 @@ auto poet::Field::AsVector() const -> poet::FieldColumn { return output; } -void poet::Field::SetFromVector(const poet::FieldColumn &cont_field) { - if (cont_field.size() != this->size() * req_vec_size) { - throw std::runtime_error( - "Field::SetFromVector: vector does not match expected size"); - } - - uint32_t vec_p = 0; - for (const auto &elem : props) { - const auto start = cont_field.begin() + vec_p; - const auto end = start + req_vec_size; - - const auto map_it = this->find(elem); - - map_it->second = FieldColumn(start, end); - - vec_p += req_vec_size; - } -} - -void poet::Field::SetFromVector( - const std::vector &cont_field) { - if (cont_field.size() != this->size()) { - throw std::runtime_error( - "Input field contains more or less elements than this container."); - } - - auto in_field_it = cont_field.begin(); - - for (const auto &elem : props) { - if (in_field_it->size() != req_vec_size) { - throw std::runtime_error( - "One vector contains more or less elements than expected."); - } - - const auto map_it = this->find(elem); - - map_it->second = *(in_field_it++); - } -} - -void poet::Field::UpdateFromField(const poet::Field &input) { +void poet::Field::update(const poet::Field &input) { for (const auto &input_elem : input) { auto it_self = this->find(input_elem.first); if (it_self == this->end()) { @@ -132,3 +114,118 @@ poet::FieldColumn &poet::Field::operator[](const std::string &key) { return std::unordered_map::operator[](key); } + +SEXP poet::Field::asSEXP() const { + const std::size_t cols = this->props.size(); + + SEXP s_names = PROTECT(Rf_allocVector(STRSXP, cols)); + SEXP s_output = PROTECT(Rf_allocVector(VECSXP, cols)); + + for (std::size_t prop_i = 0; prop_i < this->props.size(); prop_i++) { + const auto &name = this->props[prop_i]; + SEXP s_values = PROTECT(Rf_allocVector(REALSXP, this->req_vec_size)); + const auto values = this->find(name)->second; + + SET_STRING_ELT(s_names, prop_i, Rf_mkChar(name.c_str())); + + for (std::size_t i = 0; i < this->req_vec_size; i++) { + REAL(s_values)[i] = values[i]; + } + + SET_VECTOR_ELT(s_output, prop_i, s_values); + + UNPROTECT(1); + } + + SEXP s_rownames = PROTECT(Rf_allocVector(INTSXP, this->req_vec_size)); + for (std::size_t i = 0; i < this->req_vec_size; i++) { + INTEGER(s_rownames)[i] = static_cast(i + 1); + } + + Rf_setAttrib(s_output, R_ClassSymbol, + Rf_ScalarString(Rf_mkChar("data.frame"))); + Rf_setAttrib(s_output, R_NamesSymbol, s_names); + Rf_setAttrib(s_output, R_RowNamesSymbol, s_rownames); + + UNPROTECT(3); + + return s_output; +} + +poet::Field &poet::Field::operator=(const FieldColumn &cont_field) { + if (cont_field.size() != this->size() * req_vec_size) { + throw std::runtime_error( + "Field::SetFromVector: vector does not match expected size"); + } + + uint32_t vec_p = 0; + for (const auto &elem : props) { + const auto start = cont_field.begin() + vec_p; + const auto end = start + req_vec_size; + + const auto map_it = this->find(elem); + + map_it->second = FieldColumn(start, end); + + vec_p += req_vec_size; + } + + return *this; +} + +poet::Field & +poet::Field::operator=(const std::vector &cont_field) { + if (cont_field.size() != this->size()) { + throw std::runtime_error( + "Input field contains more or less elements than this container."); + } + + auto in_field_it = cont_field.begin(); + + for (const auto &elem : props) { + if (in_field_it->size() != req_vec_size) { + throw std::runtime_error( + "One vector contains more or less elements than expected."); + } + + const auto map_it = this->find(elem); + + map_it->second = *(in_field_it++); + } + return *this; +} + +void poet::Field::fromSEXP(const SEXP &s_rhs) { + this->clear(); + + SEXP s_vec = PROTECT(Rf_coerceVector(s_rhs, VECSXP)); + SEXP s_names = PROTECT(Rf_getAttrib(s_vec, R_NamesSymbol)); + + std::size_t cols = static_cast(Rf_length(s_vec)); + + this->props.clear(); + this->props.reserve(cols); + + for (std::size_t i = 0; i < cols; i++) { + const std::string prop_name(CHAR(STRING_ELT(s_names, i))); + this->props.push_back(prop_name); + + SEXP s_values = PROTECT(VECTOR_ELT(s_vec, i)); + + if (i == 0) { + this->req_vec_size = static_cast(Rf_length(s_values)); + } + + FieldColumn input(this->req_vec_size); + + for (std::size_t j = 0; j < this->req_vec_size; j++) { + input[j] = static_cast(REAL(s_values)[j]); + } + + UNPROTECT(1); + + this->insert({prop_name, input}); + } + + UNPROTECT(2); +} diff --git a/src/DiffusionModule.cpp b/src/DiffusionModule.cpp index 576785461..e02c42e03 100644 --- a/src/DiffusionModule.cpp +++ b/src/DiffusionModule.cpp @@ -18,15 +18,15 @@ ** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ -#include -#include -#include -#include #include #include #include #include #include +#include +#include +#include +#include #include #include @@ -58,17 +58,18 @@ inline const char *convert_bc_to_config_file(uint8_t in) { DiffusionModule::DiffusionModule(const poet::DiffusionParams &diffu_args, const poet::GridParams &grid_params) - : t_field{grid_params.total_n}, n_cells_per_prop(grid_params.total_n) { + : n_cells_per_prop(grid_params.total_n) { this->diff_input.setGridCellN(grid_params.n_cells[0], grid_params.n_cells[1]); this->diff_input.setDomainSize(grid_params.s_cells[0], grid_params.s_cells[1]); this->dim = grid_params.dim; - this->initialize(diffu_args); + this->initialize(diffu_args, grid_params.total_n); } -void DiffusionModule::initialize(const poet::DiffusionParams &args) { +void DiffusionModule::initialize(const poet::DiffusionParams &args, + std::uint32_t n_grid_cells) { // const poet::DiffusionParams args(this->R); // name of props @@ -87,7 +88,7 @@ void DiffusionModule::initialize(const poet::DiffusionParams &args) { this->alpha.push_back(args.alpha[this->prop_names[i]]); const double val = args.initial_t[prop_names[i]]; - std::vector init_val(t_field.GetRequestedVecSize(), val); + std::vector init_val(n_grid_cells, val); initial_values.push_back(std::move(init_val)); if (this->dim == this->DIM_2D) { @@ -100,7 +101,7 @@ void DiffusionModule::initialize(const poet::DiffusionParams &args) { } } - t_field.InitFromVec(initial_values, prop_names); + this->t_field = Field(n_grid_cells, initial_values, prop_names); // apply boundary conditions to each ghost node uint8_t bc_size = (this->dim == this->DIM_1D ? 2 : 4); @@ -156,30 +157,30 @@ void DiffusionModule::simulate(double dt) { MSG_NOENDL("DiffusionModule::simulate(): Starting diffusion ..."); std::cout << std::flush; - + std::vector> field_2d = t_field.As2DVector(); this->diff_input.setTimestep(dt); for (int i = 0; i < field_2d.size(); i++) { - std::vector in_alpha(this->n_cells_per_prop, this->alpha[i]); - this->diff_input.setBoundaryCondition(this->bc_vec[i]); - - if (this->dim == this->DIM_1D) { - tug::diffusion::BTCS_1D(this->diff_input, field_2d[i].data(), - in_alpha.data()); - } else { - tug::diffusion::ADI_2D(this->diff_input, field_2d[i].data(), - in_alpha.data()); - } + std::vector in_alpha(this->n_cells_per_prop, this->alpha[i]); + this->diff_input.setBoundaryCondition(this->bc_vec[i]); + + if (this->dim == this->DIM_1D) { + tug::diffusion::BTCS_1D(this->diff_input, field_2d[i].data(), + in_alpha.data()); + } else { + tug::diffusion::ADI_2D(this->diff_input, field_2d[i].data(), + in_alpha.data()); + } } - - t_field.SetFromVector(field_2d); - - std::cout << " done!" << std::endl; - + + t_field = field_2d; + + std::cout << " done!\n"; + sim_a_transport = MPI_Wtime(); - + transport_t += sim_a_transport - sim_b_transport; } diff --git a/src/SimParams.cpp b/src/SimParams.cpp index 761de8f3f..fd77dde20 100644 --- a/src/SimParams.cpp +++ b/src/SimParams.cpp @@ -83,14 +83,22 @@ void poet::ChemistryParams::initFromR(RInsidePOET &R) { Rcpp::as(R.parseEval("mysetup$chemistry$input_script")); if (R.checkIfExists("dht_species", "mysetup$chemistry")) { - this->dht_signifs = - R.wrapNamedVector("mysetup$chemistry$dht_species"); + this->dht_signifs = Rcpp::as>( + R.parseEval(("mysetup$chemistry$dht_species"))); } if (R.checkIfExists("pht_species", "mysetup$chemistry")) { - this->pht_signifs = - R.wrapNamedVector("mysetup$chemistry$pht_species"); + this->pht_signifs = Rcpp::as>( + R.parseEval(("mysetup$chemistry$pht_species"))); } + this->hooks.dht_fill = + RHookFunction(R, "mysetup$chemistry$hooks$dht_fill"); + this->hooks.dht_fuzz = + RHookFunction>(R, "mysetup$chemistry$hooks$dht_fuzz"); + this->hooks.interp_pre = RHookFunction>( + R, "mysetup$chemistry$hooks$interp_pre_func"); + this->hooks.interp_post = + RHookFunction(R, "mysetup$chemistry$hooks$interp_post_func"); } SimParams::SimParams(int world_rank_, int world_size_) { diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 9d3a5ac43..32078d388 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -5,6 +5,10 @@ file(GLOB test_SRC add_executable(testPOET ${test_SRC}) target_link_libraries(testPOET doctest poet_lib) +get_filename_component(TEST_RInsideSourceFile "RInsidePOET_funcs.R" REALPATH) +configure_file(testDataStructures.hpp.in testDataStructures.hpp) +target_include_directories(testPOET PRIVATE "${CMAKE_CURRENT_BINARY_DIR}") + add_custom_target(check COMMAND $ DEPENDS testPOET diff --git a/test/RInsidePOET_funcs.R b/test/RInsidePOET_funcs.R new file mode 100644 index 000000000..668b05169 --- /dev/null +++ b/test/RInsidePOET_funcs.R @@ -0,0 +1,22 @@ +simple_named_vec <- function(input) { + input["Ca"] <- 0 + + return(input) +} + +extended_named_vec <- function(input, additional) { + return(input["H"] + additional["H"]) +} + +bool_named_vec <- function(input) { + return(input["H"] == 0) +} + +simple_field <- function(field) { + field$Na <- 0 + return(field) +} + +extended_field <- function(field, additional) { + return(field + additional) +} diff --git a/test/testDataStructures.hpp.in b/test/testDataStructures.hpp.in new file mode 100644 index 000000000..c1d10271a --- /dev/null +++ b/test/testDataStructures.hpp.in @@ -0,0 +1,6 @@ +#ifndef TESTRINSIDEPOET_H_ +#define TESTRINSIDEPOET_H_ + +inline const char *RInside_source_file = "@TEST_RInsideSourceFile@"; + +#endif // TESTRINSIDEPOET_H_ diff --git a/test/testDataStructures.cpp b/test/testField.cpp similarity index 71% rename from test/testDataStructures.cpp rename to test/testField.cpp index c3046f275..8ee1bf5d6 100644 --- a/test/testDataStructures.cpp +++ b/test/testField.cpp @@ -1,12 +1,17 @@ +#include #include #include +#include #include #include #include -#include +#include +#include #include #include +#include "testDataStructures.hpp" + using namespace poet; #define CHECK_AND_FAIL_LOOP(val1, val2) \ @@ -19,14 +24,17 @@ TEST_CASE("Field") { constexpr uint32_t VEC_COUNT = 3; constexpr double INIT_VAL = 1; - Field dut(VEC_SIZE); - std::vector names = {"C", "Ca", "Na"}; std::vector init_values(names.size(), FieldColumn(VEC_SIZE, INIT_VAL)); + RInsidePOET &R = RInsidePOET::getInstance(); + + R["sourcefile"] = RInside_source_file; + R.parseEval("source(sourcefile)"); + SUBCASE("Initialize field with 2D vector") { - dut.InitFromVec(init_values, names); + Field dut(VEC_SIZE, init_values, names); auto props = dut.GetProps(); @@ -43,7 +51,7 @@ TEST_CASE("Field") { SUBCASE("Initialize field with 2D vector") { std::vector init_values(VEC_SIZE * VEC_COUNT, 1); - dut.InitFromVec(init_values, names); + Field dut(VEC_SIZE, init_values, names); auto props = dut.GetProps(); @@ -58,7 +66,7 @@ TEST_CASE("Field") { } } - dut.InitFromVec(init_values, names); + Field dut(VEC_SIZE, init_values, names); SUBCASE("Return vector") { std::vector output = dut.AsVector(); @@ -66,6 +74,37 @@ TEST_CASE("Field") { CHECK(output.size() == VEC_SIZE * VEC_COUNT); } + SUBCASE("SEXP return") { + R["test"] = dut.asSEXP(); + + Field field_constructed = R.parseEval("test"); + + Rcpp::DataFrame R_df = R.parseEval("test"); + Rcpp::CharacterVector R_names = R_df.names(); + + for (std::size_t i = 0; i < VEC_COUNT; i++) { + const auto r_values = Rcpp::as>(R_df[i]); + CHECK_EQ(r_values, dut[names[i]]); + CHECK_EQ(R_names[i], names[i]); + CHECK_EQ(field_constructed[names[i]], dut[names[i]]); + } + } + + SUBCASE("Apply R function (set Na to zero)") { + RHookFunction to_call(R, "simple_field"); + Field field_proc = to_call(dut.asSEXP()); + + CHECK_EQ(field_proc["Na"], FieldColumn(dut.GetRequestedVecSize(), 0)); + } + + SUBCASE("Apply R function (add two fields)") { + RHookFunction to_call(R, "extended_field"); + Field field_proc = to_call(dut.asSEXP(), dut.asSEXP()); + + CHECK_EQ(field_proc["Na"], + FieldColumn(dut.GetRequestedVecSize(), INIT_VAL + INIT_VAL)); + } + constexpr double NEW_VAL = 2.; std::vector new_val_vec(VEC_SIZE, NEW_VAL); @@ -102,7 +141,7 @@ TEST_CASE("Field") { } } - dut.SetFromVector(new_field); + dut = new_field; SUBCASE("SetFromVector return correct field vector") { auto ret_vec = dut.AsVector(); @@ -147,19 +186,18 @@ TEST_CASE("Field") { // reset field names = dut.GetProps(); - dut.SetFromVector( - std::vector(names.size(), FieldColumn(VEC_SIZE, INIT_VAL))); + + dut = std::vector(names.size(), FieldColumn(VEC_SIZE, INIT_VAL)); constexpr double SOME_OTHER_VAL = -0.5; - Field some_other_field(VEC_SIZE); std::vector some_other_props = {"Na", "Cl"}; std::vector some_other_values(VEC_SIZE * some_other_props.size(), SOME_OTHER_VAL); - some_other_field.InitFromVec(some_other_values, some_other_props); + Field some_other_field(VEC_SIZE, some_other_values, some_other_props); SUBCASE("Update existing field from another field") { - dut.UpdateFromField(some_other_field); + dut.update(some_other_field); auto ret_vec = dut.As2DVector(); auto ret_it = ret_vec.begin(); diff --git a/test/testNamedVector.cpp b/test/testNamedVector.cpp new file mode 100644 index 000000000..713f2729b --- /dev/null +++ b/test/testNamedVector.cpp @@ -0,0 +1,63 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "testDataStructures.hpp" + +TEST_CASE("NamedVector") { + RInsidePOET &R = RInsidePOET::getInstance(); + + R["sourcefile"] = RInside_source_file; + R.parseEval("source(sourcefile)"); + + const std::vector names{{"H", "O", "Ca"}}; + const std::vector values{{0.1, 0.2, 0.3}}; + + poet::NamedVector nv(names, values); + + SUBCASE("SEXP conversion") { + R["test"] = nv; + + const Rcpp::NumericVector R_value = R.parseEval("test"); + const Rcpp::CharacterVector R_names = R_value.names(); + + const poet::NamedVector nv_constructed = R["test"]; + + for (std::size_t i = 0; i < R_value.size(); i++) { + CHECK_EQ(R_value[i], values[i]); + CHECK_EQ(R_names[i], names[i]); + CHECK_EQ(nv_constructed[i], values[i]); + CHECK_EQ(nv_constructed.getNames()[i], names[i]); + } + } + + SUBCASE("Apply R function (set to zero)") { + RHookFunction> to_call(R, "simple_named_vec"); + nv = to_call(nv); + + CHECK_EQ(nv[2], 0); + } + + SUBCASE("Apply R function (second NamedVector)") { + RHookFunction> to_call(R, "extended_named_vec"); + + const std::vector names{{"C", "H", "Mg"}}; + const std::vector values{{0, 1, 2}}; + + poet::NamedVector second_nv(names, values); + + nv = to_call(nv, second_nv); + + CHECK_EQ(nv[0], 1.1); + } + + SUBCASE("Apply R function (check if zero)") { + RHookFunction to_call(R, "bool_named_vec"); + + CHECK_FALSE(to_call(nv)); + } +}