From 6a88de5c5ddf9cc11177b8cd733b2c6534e92caa Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Thu, 4 Apr 2024 09:27:52 +0000
Subject: [PATCH] Refactor code for grid creation and result storage

---
 R_lib/init_r_lib.R    | 32 ++++++++++++++++++--------------
 R_lib/kin_r_library.R | 11 ++++++++---
 src/poet.cpp          | 11 +++++------
 3 files changed, 31 insertions(+), 23 deletions(-)

diff --git a/R_lib/init_r_lib.R b/R_lib/init_r_lib.R
index f6b433b42..8e9bef96c 100644
--- a/R_lib/init_r_lib.R
+++ b/R_lib/init_r_lib.R
@@ -3,22 +3,26 @@ has_doParallel <- require(doParallel)
 
 seq_pqc_to_grid <- function(pqc_in, grid) {
     # Convert the input DataFrame to a matrix
-    dt <- as.matrix(pqc_in)
+    pqc_in <- as.matrix(pqc_in)
 
     # Flatten the matrix into a vector
-    id_vector <- as.vector(t(grid))
+    id_vector <- as.numeric(t(grid))
 
     # Initialize an empty matrix to store the results
-    result_mat <- matrix(nrow = 0, ncol = ncol(dt))
+    # result_mat <- matrix(NA, nrow = length(id_vector), ncol = ncol(pqc_in))
+
+    row_indices <- match(id_vector, pqc_in[, "ID"])
+
+    result_mat <- pqc_in[row_indices, ]
 
     # Iterate over each ID in the vector
-    for (id_mat in id_vector) {
-        # Find the matching row in the matrix
-        matching_row <- dt[dt[, "ID"] == id_mat, ]
+    # for (i in seq_along(id_vector)) {
+    #     # Find the matching row in the matrix
+    #     # matching_row <- pqc_in[pqc_in[, "ID"] == i, ]
 
-        # Append the matching row to the result matrix
-        result_mat <- rbind(result_mat, matching_row)
-    }
+    #     # Append the matching row to the result matrix
+    #     result_mat[i, ] <- pqc_in[pqc_in[, "ID"] == i, ]
+    # }
 
     # Convert the result matrix to a data frame
     res_df <- as.data.frame(result_mat)
@@ -72,13 +76,13 @@ par_pqc_to_grid <- function(pqc_in, grid) {
 }
 
 pqc_to_grid <- function(pqc_in, grid) {
-    if (has_doParallel && has_foreach) {
-        print("Using parallel grid creation")
-        return(par_pqc_to_grid(pqc_in, grid))
-    } else {
+    # if (has_doParallel && has_foreach) {
+    #     print("Using parallel grid creation")
+    #     return(par_pqc_to_grid(pqc_in, grid))
+    # } else {
         print("Using sequential grid creation")
         return(seq_pqc_to_grid(pqc_in, grid))
-    }
+    # }
 }
 
 resolve_pqc_bound <- function(pqc_mat, transport_spec, id) {
diff --git a/R_lib/kin_r_library.R b/R_lib/kin_r_library.R
index 33e6ef1dc..b0e34ca61 100644
--- a/R_lib/kin_r_library.R
+++ b/R_lib/kin_r_library.R
@@ -15,7 +15,7 @@
 ### this program; if not, write to the Free Software Foundation, Inc., 51
 ### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 
-master_init <- function(setup, out_dir) {
+master_init <- function(setup, out_dir, init_field) {
     ## Setup the directory where we will store the results
     if (!dir.exists(out_dir)) {
         dir.create(out_dir)
@@ -23,7 +23,7 @@ master_init <- function(setup, out_dir) {
     } else {
         msgm("dir ", out_dir, " already exists, I will overwrite!")
     }
-    if (!exists("setup$store_result")) {
+    if (is.null(setup$store_result)) {
         msgm("store_result doesn't exist!")
     } else {
         msgm("store_result is ", setup$store_result)
@@ -40,11 +40,16 @@ master_init <- function(setup, out_dir) {
     }
 
     if (setup$store_result) {
+        init_field_out <- paste0(out_dir, "/iter_0.rds")
+        saveRDS(init_field, file = init_field_out)
+        msgm("Stored initial field in ", init_field_out)
         if (is.null(setup[["out_save"]])) {
             setup$out_save <- seq(1, setup$iterations)
         }
     }
 
+    setup$out_dir <- out_dir
+
     return(setup)
 }
 
@@ -63,7 +68,7 @@ master_iteration_end <- function(setup, state_T, state_C) {
     ## comprised in setup$out_save
     if (setup$store_result) {
         if (iter %in% setup$out_save) {
-            nameout <- paste0(fileout, "/iter_", sprintf(fmt = fmt, iter), ".rds")
+            nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".rds")
             saveRDS(list(
                 T = state_T, C = state_C,
                 simtime = as.integer(setup$simulation_time)
diff --git a/src/poet.cpp b/src/poet.cpp
index 253062c88..8ee5585a1 100644
--- a/src/poet.cpp
+++ b/src/poet.cpp
@@ -31,7 +31,6 @@
 #include <Rcpp.h>
 #include <Rcpp/Function.h>
 #include <Rcpp/vector/instantiation.h>
-#include <cstdio>
 #include <cstdlib>
 #include <memory>
 #include <mpi.h>
@@ -386,7 +385,8 @@ int main(int argc, char *argv[]) {
       // // if (MY_RANK == 0) { // get timestep vector from
       // // grid_init function ... //
       *global_rt_setup =
-          master_init_R.value()(*global_rt_setup, run_params.out_dir);
+          master_init_R.value()(*global_rt_setup, run_params.out_dir,
+                                init_list.getInitialGrid().asSEXP());
 
       // MDL: store all parameters
       // MSG("Calling R Function to store calling parameters");
@@ -405,13 +405,12 @@ int main(int argc, char *argv[]) {
 
       MSG("finished simulation loop");
 
-      // R["simtime"] = dSimTime;
-      // R.parseEvalQ("profiling$simtime <- simtime");
-
       R["profiling"] = profiling;
+      R["setup"] = *global_rt_setup;
 
       string r_vis_code;
-      r_vis_code = "saveRDS(profiling, file=paste0(fileout,'/timings.rds'));";
+      r_vis_code =
+          "saveRDS(profiling, file=paste0(setup$out_dir,'/timings.rds'));";
       R.parseEval(r_vis_code);
 
       MSG("Done! Results are stored as R objects into <" + run_params.out_dir +