Merge branch 'rsachetto-master'

bsbash/build_functions.sh

@@ -21,7 +21,7 @@ LIBRARY_OUTPUT_DIRECTORY=$ROOT_DIR
 GLOBAL_FORCE_COMPILATION=""
 QUIET=''
 WAIT_ENTER=''
-WRITE_COMPILE_COMMANDS='y'
+WRITE_COMPILE_COMMANDS=''
 
 DEFAULT_BUILD_DIR="build_"
 COMPILE_COMMANDS_FILE="${ROOT_DIR}/compile_commands.json"
@@ -397,9 +397,6 @@ COMPILE_OBJECT () {
 		ANY_COMPILED="y"
 
 	fi
-
-
-
 }
 
 CHECK_HEADERS() {

example_configs/rabbit_mesh_eikonal.ini

@@ -0,0 +1,35 @@
+[main]
+simulation_time = 100
+dt = 1
+
+[domain]
+;These values are mandatory
+name=UCLA Rabbit Mesh
+;this mesh always start at 250.0
+maximum_discretization = 250.0
+main_function=initialize_grid_with_rabbit_mesh
+;These can be optional depending on the domain main_function
+mesh_file=meshes/rabheart.alg
+
+[stim_point]
+current = 1
+center_x=31875
+center_y=10875
+center_z=16375
+start = 0.0
+duration = 2.0
+main_function= stim_if_point_equal
+
+;[stim_mouse]
+;x_limit=500.0
+;start = 0.0
+;duration = 2.0
+;current = -50.0
+;main_function= stim_if_x_less_than
+
+[save_result]
+;/////mandatory/////////
+output_dir=./outputs/rabbit_eikonal
+main_function=save_as_text_or_binary
+;//////////////////
+file_prefix=T

src/config/config_parser.h

@@ -246,6 +246,7 @@ void free_batch_options(struct batch_options * options);
 void free_visualization_options(struct visualization_options * options);
 void free_conversion_options(struct conversion_options *options);
 void free_fibers_conversion_options(struct fibers_conversion_options *options);
+void free_eikonal_options(struct eikonal_options *options);
 
 void set_or_overwrite_common_data(struct config* config, const char *key, const char *value, const char *section, const char *config_file);
 

src/eikonal/eikonal_solver.c

@@ -1,3 +1,7 @@
+//
+// Created by sachetto on 06/06/24.
+// Based on https://github.com/SCIInstitute/StructuredEikonal
+
 #include <assert.h>
 #include <cuda_runtime.h>
 #include "eikonal_solver.h"
@@ -28,7 +32,6 @@ struct eikonal_solver * new_eikonal_solver(bool verbose) {
 
 void free_eikonal_solver(struct eikonal_solver *solver) {
 
-    //TODO: free 3D arrays
     FREE_3D_ARRAY(solver->speeds, solver->width, solver->height);
     FREE_3D_ARRAY(solver->answer, solver->width, solver->height);
     FREE_3D_ARRAY(solver->mask, solver->width, solver->height);
@@ -251,19 +254,26 @@ static void get_solution(struct eikonal_solver *solver) {
         for(int k = 0; k < BLOCK_LENGTH; k++) {
             for(int j = 0; j < BLOCK_LENGTH; j++) {
                 for(int i = 0; i < BLOCK_LENGTH; i++) {
-                    double d = solver->memory_struct.h_sol[baseAddr +
-                                                           k * BLOCK_LENGTH * BLOCK_LENGTH +
-                                                           j * BLOCK_LENGTH + i];
+                    double d = solver->memory_struct.h_sol[baseAddr + k * BLOCK_LENGTH * BLOCK_LENGTH + j * BLOCK_LENGTH + i];
                     if ((i + bx * BLOCK_LENGTH) < solver->width &&
                         (j + by * BLOCK_LENGTH) < solver->height &&
                         (k + bz * BLOCK_LENGTH) < solver->depth) {
-                        solver->answer[(i + bx * BLOCK_LENGTH)][(j +
-                                                                 by * BLOCK_LENGTH)][k + bz * BLOCK_LENGTH] = d;
+                        solver->answer[(i + bx * BLOCK_LENGTH)][(j + by * BLOCK_LENGTH)][k + bz * BLOCK_LENGTH] = d;
                     }
                 }
             }
         }
     }
+
+    for(int i = 0 ; i < solver->num_active_cells; i++) {
+        solver->active_cells[i]->v = solver->answer[(int)solver->active_cells[i]->center.x][(int)solver->active_cells[i]->center.y][(int)solver->active_cells[i]->center.z];
+
+        //Translating back to original space
+        solver->active_cells[i]->center.x = solver->active_cells[i]->center.x * solver->active_cells[i]->discretization.x + solver->active_cells[i]->discretization.x/2.0;
+        solver->active_cells[i]->center.y  = solver->active_cells[i]->center.y * solver->active_cells[i]->discretization.y + solver->active_cells[i]->discretization.y/2.0;
+        solver->active_cells[i]->center.z = solver->active_cells[i]->center.z * solver->active_cells[i]->discretization.z + solver->active_cells[i]->discretization.z/2.0;
+
+    }
 }
 
 void map_generator(struct eikonal_solver *solver) {
@@ -389,29 +399,6 @@ void use_seeds(struct eikonal_solver *solver) {
     check_cuda_error( cudaMemset(solver->memory_struct.d_con, 1, volSize*sizeof(bool)) );
 }
 
-void write_alg(struct eikonal_solver *solver, char *filename) {
-
-    FILE *out = fopen(filename, "w");
-
-    if (out == NULL) {
-        log_error_and_exit("Error opening file: %s\n", filename);
-    }
-
-    for (int k = 0; k < solver->depth; k++) {
-        for (int j = 0; j < solver->height; j++) {
-            for (int i = 0; i < solver->width; i++) {
-                if(solver->mask[i][j][k] == true) {
-                    real d = solver->answer[i][j][k];
-                    fprintf(out, "%lf, %lf, %lf, 0.5, 0.5, 0.5, %lf\n", i + 0.5, j + 0.5, k + 0.5, d);
-                }
-            }
-        }
-    }
-
-    fclose(out);
-
-}
-
 void solve_eikonal(struct eikonal_solver *solver) {
 
     if (solver->speeds == NULL) {

src/main_eikonal.c

@@ -1,29 +1,30 @@
-#include <stdlib.h>
+#include "3dparty/ini_parser/ini.h"
+#include "3dparty/stb_ds.h"
 #include "alg/cell/cell.h"
 #include "alg/grid/grid.h"
-#include "3dparty/ini_parser/ini.h"
+#include "config/config_parser.h"
 #include "config/domain_config.h"
 #include "eikonal/eikonal_solver.h"
-#include "3dparty/stb_ds.h"
+#include "config/stim_config.h"
 #include "logger/logger.h"
-#include "config/config_parser.h"
-
+#include "utils/file_utils.h"
+#include <stdlib.h>
 
 static int compare_coordinates(const void *a, const void *b) {
 
     struct cell_node *coord1 = *(struct  cell_node **) a;
     struct cell_node *coord2 = *(struct  cell_node **) b;
 
     if (coord1->center.y != coord2->center.y) {
-        return coord1->center.y - coord2->center.y;
+        return (int) (coord1->center.y - coord2->center.y);
     }
 
     if (coord1->center.x != coord2->center.x) {
-          return coord1->center.x - coord2->center.x;
+        return (int) (coord1->center.x - coord2->center.x);
     }
 
 
-    return coord1->center.z - coord2->center.z;
+    return (int) (coord1->center.z - coord2->center.z);
 
 }
 
@@ -32,6 +33,9 @@ int main(int argc, char *argv[]) {
     struct eikonal_options *eikonal_options = new_eikonal_options();
     parse_eikonal_options(argc, argv, eikonal_options);
     struct grid *grid = new_grid();
+    struct string_voidp_hash_entry *stimuli_configs = eikonal_options->stim_configs;
+
+    struct eikonal_solver *solver = new_eikonal_solver(false);
 
     if(ini_parse(eikonal_options->config_file, parse_eikonal_config_file, eikonal_options) < 0) {
         log_error_and_exit("Error: Can't load the config file %s\n", eikonal_options->config_file);
@@ -59,10 +63,61 @@ int main(int argc, char *argv[]) {
 
     }
 
+    if(stimuli_configs) {
+
+        size_t n = shlen(stimuli_configs);
+        struct time_info time_info = {0.0, 0.0, 0.0, 0};
+
+        if(n > 0) {
+            STIM_CONFIG_HASH_FOR_INIT_FUNCTIONS(stimuli_configs);
+        }
+
+        for(int i = 0; i < n; i++) {
+
+            struct string_voidp_hash_entry e = stimuli_configs[i];
+            log_info("Stimulus name: %s\n", e.key);
+            struct config *tmp = (struct config *)e.value;
+            print_stim_config_values(tmp);
+            log_msg(LOG_LINE_SEPARATOR);
+        }
+
+        set_spatial_stim(&time_info, stimuli_configs, grid, false);
+
+        struct config *tmp = NULL;
+
+        uint32_t n_active = grid->num_active_cells;
+        struct cell_node **ac = grid->active_cells;
+
+        for(size_t k = 0; k < n; k++) {
+
+            tmp = (struct config *)stimuli_configs[k].value;
+
+            for(uint32_t i = 0; i < n_active; i++) {
+                real value =  ((real *)(tmp->persistent_data))[i];
+                if(value != 0.0) {
+                    size_t *initial_seed = calloc(sizeof(size_t), 3);
+                    double new_pos_x = (ac[i]->center.x-ac[i]->discretization.x/2.0)/ac[i]->discretization.x;
+                    double new_pos_y = (ac[i]->center.y-ac[i]->discretization.y/2.0)/ac[i]->discretization.y;
+                    double new_pos_z = (ac[i]->center.z-ac[i]->discretization.z/2.0)/ac[i]->discretization.z;
+
+                    initial_seed[0] = new_pos_x;
+                    initial_seed[1] = new_pos_y;
+                    initial_seed[2] = new_pos_z;
+                    arrput(solver->seeds, initial_seed);
+                }
+            }
+        }
+    }
+
+    //Translating to 0
     for(int i = 0 ; i < grid->num_active_cells; i++) {
-        int new_pos_x = (grid->active_cells[i]->center.x-grid->active_cells[i]->discretization.x/2)/grid->active_cells[i]->discretization.x;
-        int new_pos_y = (grid->active_cells[i]->center.y-grid->active_cells[i]->discretization.y/2)/grid->active_cells[i]->discretization.y;
-        int new_pos_z = (grid->active_cells[i]->center.z-grid->active_cells[i]->discretization.z/2)/grid->active_cells[i]->discretization.z;
+        double new_pos_x = (grid->active_cells[i]->center.x-grid->active_cells[i]->discretization.x/2.0)/grid->active_cells[i]->discretization.x;
+        double new_pos_y = (grid->active_cells[i]->center.y-grid->active_cells[i]->discretization.y/2.0)/grid->active_cells[i]->discretization.y;
+        double new_pos_z = (grid->active_cells[i]->center.z-grid->active_cells[i]->discretization.z/2.0)/grid->active_cells[i]->discretization.z;
+
+        if(new_pos_x != floor(new_pos_x) || new_pos_y != floor(new_pos_y) || new_pos_z != floor(new_pos_z)) {
+            log_error_and_exit("The current version only accepts integer coordinates\n");
+        }
 
         grid->active_cells[i]->center.x = new_pos_x;
         grid->active_cells[i]->center.y = new_pos_y;
@@ -71,26 +126,39 @@ int main(int argc, char *argv[]) {
 
     size_t itersPerBlock = 10, type = 1;
 
-    struct eikonal_solver *solver = new_eikonal_solver(false);
-    solver->width = grid->cube_side_length.x/grid->active_cells[0]->discretization.x;
-    solver->height = grid->cube_side_length.y/grid->active_cells[0]->discretization.y;
-    solver->depth = grid->cube_side_length.z/grid->active_cells[0]->discretization.z;
+    solver->width  = (int) (grid->cube_side_length.x/grid->active_cells[0]->discretization.x);
+    solver->height = (int) (grid->cube_side_length.y/grid->active_cells[0]->discretization.y);
+    solver->depth  = (int) (grid->cube_side_length.z/grid->active_cells[0]->discretization.z);
 
     solver->solver_type = type;
     solver->iters_per_block = itersPerBlock;
 
-    size_t *initial_seed = calloc(sizeof(size_t), 3);
-    initial_seed[0] = grid->active_cells[0]->center.x;
-    initial_seed[1] = grid->active_cells[0]->center.y;
-    initial_seed[2] = grid->active_cells[0]->center.z;
-
-    arrput(solver->seeds, initial_seed);
     solver->active_cells = grid->active_cells;
     solver->num_active_cells = grid->num_active_cells;
 
     solve_eikonal(solver);
-    write_alg(solver, "test.alg");
+
+    struct config * save_result_config = eikonal_options->save_mesh_config;
+
+    if(save_result_config) {
+
+        init_config_functions(save_result_config, "./shared_libs/libdefault_save_mesh.so", "save_result");
+
+        print_save_mesh_config_values(save_result_config);
+        log_msg(LOG_LINE_SEPARATOR);
+
+        char *out_dir_name = NULL;
+        GET_PARAMETER_STRING_VALUE_OR_USE_DEFAULT(out_dir_name, save_result_config, "output_dir");
+        if(out_dir_name != NULL) {
+            create_dir(out_dir_name);
+            struct time_info ti = ZERO_TIME_INFO;
+            ((save_mesh_fn *)save_result_config->main_function)(&ti, save_result_config, grid, NULL, NULL);
+        } else {
+            log_warn("Not output dir provided. The result will not be saved!");
+        }
+    }
 
     free_eikonal_solver(solver);
+    free_eikonal_options(eikonal_options);
 
 }

src/stimuli_library/stimuli.c

@@ -11,6 +11,48 @@
 #include "../config_helpers/config_helpers.h"
 #include "../utils/utils.h"
 
+SET_SPATIAL_STIM(stim_if_point_equal) {
+
+    uint32_t n_active = the_grid->num_active_cells;
+    struct cell_node **ac = the_grid->active_cells;
+
+    if(is_purkinje) {
+        n_active = the_grid->purkinje->num_active_purkinje_cells;
+        ac = the_grid->purkinje->purkinje_cells;
+    }
+
+    ALLOCATE_STIMS();
+
+    real stim_current = 0.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real, stim_current, config, "current");
+
+    bool stim;
+    real stim_value;
+
+    real_cpu center_x = 0.0;
+    real_cpu center_y = 0.0;
+    real_cpu center_z = 0.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, center_x, config, "center_x");
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, center_y, config, "center_y");
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, center_z, config, "center_z");
+
+    uint32_t i;
+
+    OMP(parallel for private(stim, stim_value))
+    for(i = 0; i < n_active; i++) {
+        stim = ac[i]->center.x == center_x && ac[i]->center.y == center_y && ac[i]->center.z == center_z;
+
+        if(stim) {
+            stim_value = stim_current;
+        } else {
+            stim_value = 0.0;
+        }
+
+        SET_STIM_VALUE(i, stim_value);
+    }
+
+}
+
 SET_SPATIAL_STIM(stim_if_x_less_than) {
 
     uint32_t n_active = the_grid->num_active_cells;

src/tests/build.sh

@@ -3,11 +3,15 @@ if [ -n "$CUDA_FOUND" ]; then
     TEST_OPT_DEPS=gpu_utils
 fi
 
-	TESTS_DYNAMIC_DEPS="$CUDA_LIBRARIES $CRITERION_LIBRARIES dl m logger ${TEST_OPT_DEPS}"
+TESTS_DYNAMIC_DEPS="$CUDA_LIBRARIES $CRITERION_LIBRARIES dl m logger ${TEST_OPT_DEPS}"
+
+if [ -n "$AMGX_FOUND" ]; then
+    TESTS_DYNAMIC_DEPS="$TESTS_DYNAMIC_DEPS $AMGX_LIBRARIES"
+fi
 
 ##Tests
 TESTS_STATIC_DEPS="monodomain ode_solver config tinyexpr config_helpers alg graph utils sds"
-COMPILE_EXECUTABLE "TestSolvers" "test_solvers.c" "" "$TESTS_STATIC_DEPS" "$TESTS_DYNAMIC_DEPS $AMGX_LIBRARIES" "$CUDA_LIBRARY_PATH $AMGX_LIBRARY_PATH $CRITERION_LIBRARY_PATH $LIBRARY_OUTPUT_DIRECTORY" "$CRITERION_INCLUDE_PATH"
+COMPILE_EXECUTABLE "TestSolvers" "test_solvers.c" "" "$TESTS_STATIC_DEPS" "$TESTS_DYNAMIC_DEPS" "$CUDA_LIBRARY_PATH $AMGX_LIBRARY_PATH $CRITERION_LIBRARY_PATH $LIBRARY_OUTPUT_DIRECTORY" "$CRITERION_INCLUDE_PATH"
 
 TESTS_STATIC_DEPS="config alg graph utils sds"
 COMPILE_EXECUTABLE "TestMesh" "test_mesh.c" "" "$TESTS_STATIC_DEPS" "$TESTS_DYNAMIC_DEPS" "$CUDA_LIBRARY_PATH $CRITERION_LIBRARY_PATH $LIBRARY_OUTPUT_DIRECTORY" "$CRITERION_INCLUDE_PATH"