Skip to content

Commit

Permalink
Add a new branch/Fixes to the adjusted ToRORd model/Working tuneCV sc…
Browse files Browse the repository at this point in the history 
…ript with the new ToRORd model
  • Loading branch information
@bergolho
bergolho committed May 31, 2024
1 parent 33b72e9 commit 2a96e66
Show file tree
Hide file tree
Showing 14 changed files with 3,683 additions and 79 deletions.
168 changes: 99 additions & 69 deletions scripts/tuneCV/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,10 +11,6 @@
#include <vtkXMLPolyDataWriter.h>
#include <vtkSmartPointer.h>
#include <vtkDataSetMapper.h>
#include <vtkActor.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkUnstructuredGrid.h>
#include <vtkHexahedron.h>
#include <vtkSphereSource.h>
Expand All @@ -23,57 +19,87 @@
#include <vtkCellData.h>
#include <vtkLine.h>
#include <vtkPointData.h>
#include <vtkCellLocator.h>

using namespace std;

// Change your MonoAlg3D path here:
// ----------------------------------------------------------
const char MONOALG_PATH[500] = "/home/berg/Github/MonoAlg3D_C";
// ----------------------------------------------------------

const double TOLERANCE = 1.0e-03;

double calculate_conduction_velocity_from_simulation ()
{
string filename = "outputs/cable/purkinje_activation_time_map_pulse_it_0.vtp";
string filename = "outputs/cable/tissue_activation_time_map_pulse_it_0.vtu";

// Read all the data from the file
vtkSmartPointer<vtkXMLPolyDataReader> reader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
vtkSmartPointer<vtkXMLUnstructuredGridReader> reader = vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
reader->SetFileName(filename.c_str());
reader->Update();

vtkPolyData* polydata = reader->GetOutput();
uint32_t num_points = polydata->GetNumberOfPoints();
uint32_t num_lines = polydata->GetNumberOfLines();
vtkUnstructuredGrid* unstructuredGrid = reader->GetOutput();
uint32_t num_points = unstructuredGrid->GetNumberOfPoints();
uint32_t num_lines = unstructuredGrid->GetNumberOfCells();

vtkSmartPointer<vtkCellLocator> cellLocator = vtkSmartPointer<vtkCellLocator>::New();
cellLocator->SetDataSet(unstructuredGrid);
cellLocator->BuildLocator();

// Points of interest
double point_x_0[3] = {7500, 250, 250};
double point_x_1[3] = {12500, 250, 250};

// Find (closest points) Cell indexes for CV computation
vtkIdType cellId_x_0; // the cell id of the cell containing the closest point
vtkIdType cellId_x_1; // the cell id of the cell containing the closest point

double closestPoint_x_0[3];
double closestPoint_x_1[3];

int subId; // not needed
double closestPointDist2; // not needed

cellLocator->FindClosestPoint(point_x_0, closestPoint_x_0, cellId_x_0, subId, closestPointDist2);
cellLocator->FindClosestPoint(point_x_1, closestPoint_x_1, cellId_x_1, subId, closestPointDist2);

double delta_s_x = sqrt(pow(closestPoint_x_0[0]-closestPoint_x_1[0], 2)); // only changes over x-axis

cout << delta_s_x << endl;

// Read points scalar values
string array_name = "Scalars_";
vtkSmartPointer<vtkFloatArray> array = vtkFloatArray::SafeDownCast(polydata->GetPointData()->GetArray(array_name.c_str()));

double cv = 0.0;
if(array)
{
// Cell indexes for CV computation
uint32_t middle_cell = num_points/2;
uint32_t prev_cell = middle_cell-25;
uint32_t next_cell = middle_cell+25;
vtkSmartPointer<vtkFloatArray> array = vtkFloatArray::SafeDownCast(unstructuredGrid->GetCellData()->GetArray(array_name.c_str()));

double cv_x = -1.0;

if(array) {
double delta_lat_x = (array->GetValue(cellId_x_1) - array->GetValue(cellId_x_0)); // ms

double delta_s = 5000.0;
double delta_t = (array->GetValue(next_cell) - array->GetValue(prev_cell));
cv = (delta_s / delta_t)*0.001; // {m/s}
cout << delta_lat_x << endl;
cout << array->GetValue(cellId_x_0) << endl;
cout << array->GetValue(cellId_x_1) << endl;

cv_x = (delta_s_x / delta_lat_x)*0.001; // {m/s}
}
else
{
cerr << "[!] ERROR! No Scalar_value found for the points!" << endl;
else {
cerr << "[!] ERROR! No 'Scalar_value' found for the points!" << endl;
exit(EXIT_FAILURE);
}

return cv;
return cv_x;
}

// TODO: Maybe pass a pre-configured config file as an input parameter with the cellular model setup that the user will use
void write_configuration_file (const double sigma)
{
FILE *file = fopen("/home/berg/Github/MonoAlg3D_C/scripts/tuneCV/configs/cable.ini","w+");
void write_configuration_file (const double sigma) {
char filename[500];
sprintf(filename, "%s/scripts/tuneCV/configs/cable.ini", MONOALG_PATH);
FILE *file = fopen(filename,"w+");

fprintf(file,"[main]\n");
fprintf(file,"num_threads=4\n");
fprintf(file,"dt_pde=0.02\n");
fprintf(file,"dt_pde=0.01\n");
fprintf(file,"simulation_time=150.0\n");
fprintf(file,"abort_on_no_activity=false\n");
fprintf(file,"use_adaptivity=false\n");
Expand All @@ -82,80 +108,83 @@ void write_configuration_file (const double sigma)

fprintf(file,"[update_monodomain]\n");
fprintf(file,"main_function=update_monodomain_default\n");
fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_update_monodomain.so\n");
fprintf(file,"library_file=%s/shared_libs/libdefault_update_monodomain.so\n", MONOALG_PATH);
fprintf(file,"\n");

// For saving the LATs in a format that can be read for calculating the CVs
fprintf(file,"[save_result]\n");
fprintf(file,"print_rate=1\n");
fprintf(file,"output_dir=/home/berg/Github/MonoAlg3D_C/scripts/tuneCV/outputs/cable\n");
fprintf(file,"output_dir=%s/scripts/tuneCV/outputs/cable\n", MONOALG_PATH);
fprintf(file,"save_pvd=true\n");
fprintf(file,"file_prefix=V\n");
fprintf(file,"save_activation_time=true\n");
fprintf(file,"save_apd=true\n");
fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_save_mesh_purkinje.so\n");
fprintf(file,"main_function=save_purkinje_with_activation_times\n");
fprintf(file,"init_function=init_save_purkinje_with_activation_times\n");
fprintf(file,"end_function=end_save_purkinje_with_activation_times\n");
fprintf(file,"activation_threshold=-50.0\n");
fprintf(file,"save_apd=false\n");
fprintf(file,"library_file=%s/shared_libs/libdefault_save_mesh_purkinje.so\n", MONOALG_PATH);
fprintf(file,"main_function=save_tissue_with_activation_times\n");
fprintf(file,"init_function=init_save_tissue_with_activation_times\n");
fprintf(file,"end_function=end_save_tissue_with_activation_times\n");
fprintf(file,"remove_older_simulation=true\n");
fprintf(file,"\n");

fprintf(file,"[assembly_matrix]\n");
fprintf(file,"init_function=set_initial_conditions_fvm\n");
fprintf(file,"sigma_purkinje=%g\n",sigma);
fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libpurkinje_matrix_assembly.so\n");
fprintf(file,"main_function=purkinje_fibers_assembly_matrix\n");
fprintf(file,"sigma_x=%g\n",sigma);
fprintf(file,"sigma_y=%g\n",sigma);
fprintf(file,"sigma_z=%g\n",sigma);
fprintf(file,"library_file=%s/shared_libs/libdefault_matrix_assembly.so\n", MONOALG_PATH);
fprintf(file,"main_function=homogeneous_sigma_assembly_matrix\n");
fprintf(file,"\n");

fprintf(file,"[linear_system_solver]\n");
fprintf(file,"tolerance=1e-16\n");
fprintf(file,"use_preconditioner=yes\n");
fprintf(file,"max_iterations=500\n");
fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_linear_system_solver.so\n");
fprintf(file,"library_file=%s/shared_libs/libdefault_linear_system_solver.so\n", MONOALG_PATH);
fprintf(file,"main_function=conjugate_gradient\n");
fprintf(file,"\n");

fprintf(file,"[purkinje]\n");
fprintf(file,"name=Simple Purkinje\n");
//fprintf(file,"dx=100.0\n");
fprintf(file,"dx=400.0\n");
fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_purkinje.so\n");
fprintf(file,"main_function=initialize_purkinje_with_custom_mesh\n");
//fprintf(file,"network_file=/home/berg/Github/MonoAlg3D_C/networks/simple_cable_10cm.vtk\n");
fprintf(file,"network_file=/home/berg/Github/MonoAlg3D_C/networks/simple_cable_40cm.vtk\n");
fprintf(file,"[domain]\n");
fprintf(file,"name=Simple Cable\n");
fprintf(file,"start_dx=500.0\n");
fprintf(file,"start_dy=500.0\n");
fprintf(file,"start_dz=500.0\n");
fprintf(file,"cable_length=20000.0\n");
fprintf(file,"library_file=%s/shared_libs/libdefault_domains.so\n", MONOALG_PATH);
fprintf(file,"main_function=initialize_grid_with_cable_mesh\n");
fprintf(file,"\n");

fprintf(file,"[ode_solver]\n");
//fprintf(file,"adaptive=true\n");
fprintf(file,"dt=0.02\n");
fprintf(file,"dt=0.01\n");
fprintf(file,"use_gpu=no\n");
fprintf(file,"gpu_id=0\n");
//fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libtrovato_2019.so\n");
//fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libToRORd_fkatp_endo.so\n");
fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libten_tusscher_3_endo.so\n");
fprintf(file,"library_file=%s/shared_libs/libToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.so\n", MONOALG_PATH);
fprintf(file,"\n");

fprintf(file,"[purkinje_stim_his]\n");
fprintf(file,"[stim_benchmark]\n");
fprintf(file,"start = 0.0\n");
fprintf(file,"duration = 5.0\n");
//fprintf(file,"current = -40.0\n");
fprintf(file,"current = -53.0\n");
fprintf(file,"id_limit = 25\n");
fprintf(file,"main_function=stim_if_id_less_than\n");
fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_stimuli.so\n");
fprintf(file,"duration = 2.0\n");
fprintf(file,"current = -50.0\n");
fprintf(file, "min_x = 0.0\n");
fprintf(file, "max_x = 500.0\n");
fprintf(file, "min_y = 0.0\n");
fprintf(file, "max_y = 500.0\n");
fprintf(file, "min_z = 0.0\n");
fprintf(file, "max_z = 3000.0\n");
fprintf(file,"main_function=stim_x_y_z_limits\n");
fprintf(file,"library_file=%s/shared_libs/libdefault_stimuli.so\n", MONOALG_PATH);
fprintf(file,"\n");

fclose(file);
}

double calculate_error (const double cv, const double target_cv)
{
double calculate_error (const double cv, const double target_cv) {
return sqrt(pow(cv,2)-pow(target_cv,2))/target_cv;
}

int main (int argc, char *argv[])
{
if (argc-1 != 1)
{
int main (int argc, char *argv[]) {
if (argc-1 != 1) {

cerr << "=============================================================================" << endl;
cerr << "Usage:> " << argv[0] << " <target_CV>" << endl;
cerr << "=============================================================================" << endl;
Expand All @@ -175,12 +204,13 @@ int main (int argc, char *argv[])
double target_cv = atof(argv[1]);
double sigma = 0.0002;

do
{
do {
write_configuration_file(sigma);

// Run the simulation
system("/home/berg/Github/MonoAlg3D_C/bin/MonoAlg3D -c /home/berg/Github/MonoAlg3D_C/scripts/tuneCV/configs/cable.ini");
char command[500];
sprintf(command,"%s/bin/MonoAlg3D -c %s/scripts/tuneCV/configs/cable.ini", MONOALG_PATH, MONOALG_PATH);
system(command);

cv = calculate_conduction_velocity_from_simulation();
factor = pow(target_cv/cv,2);
Expand Down
5 changes: 5 additions & 0 deletions scripts/tuneCVbenchmark/.gitignore
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
/tables
build/*
outputs/*
bin/*
configs/*
12 changes: 12 additions & 0 deletions scripts/tuneCVbenchmark/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
cmake_minimum_required(VERSION 2.8)

PROJECT(tuneCV)

find_package(VTK REQUIRED)
include(${VTK_USE_FILE})

set( CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin )

add_executable(tuneCV main.cpp )

target_link_libraries(tuneCV ${VTK_LIBRARIES})
15 changes: 15 additions & 0 deletions scripts/tuneCVbenchmark/clean_project.sh
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
#!/bin/bash

BUILD_DIR="build"

if [[ "$#" -eq 1 ]]; then
BUILD_DIR=$1
fi

if [[ ! -d "${BUILD_DIR}" ]]; then
echo "Directory ${BUILD_DIR} does not exist"
exit
fi

rm -fr ${BUILD_DIR}/*
rm -fr bin/*
Loading

0 comments on commit 2a96e66

Please sign in to comment.