From 8f4049f58dc7d4bde5fabfcc7ae812a47fb3ce16 Mon Sep 17 00:00:00 2001
From: Kyle <kyle.c.klenk@gmail.com>
Date: Fri, 30 Jun 2023 11:29:19 -0600
Subject: [PATCH] Add Actor to find s2
---
.gitignore | 1 +
environment_setup/Dockerfile | 1 +
find_s1.py | 49 ++++-----
find_s2.py | 9 +-
legacy_scripts/find_s1_example.py | 101 +++++++++++++++++
main.cpp | 174 +++++++++++++++++++++++++++---
6 files changed, 287 insertions(+), 48 deletions(-)
create mode 100644 legacy_scripts/find_s1_example.py
diff --git a/.gitignore b/.gitignore
index 8b90690..36d6dd3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,4 @@ meshes/
output8.txt
.vscode/*
simulations/logs/*
+.DS_Store
diff --git a/environment_setup/Dockerfile b/environment_setup/Dockerfile
index 1727d45..c2edf64 100755
--- a/environment_setup/Dockerfile
+++ b/environment_setup/Dockerfile
@@ -70,6 +70,7 @@ RUN cmake --build _build
WORKDIR /opt/openCARP/external/carputils
RUN python3 -m pip install --user -r requirements.txt
ENV PATH=$PATH:/opt/openCARP/external/carputils/bin
+ENV PATH=$PATH:/opt/openCARP/_build/bin
ENV PYTHONPATH=$PYTHONPATH:/opt/openCARP/external/carputils/
# Create a settings.yaml - add your own settings.yaml
diff --git a/find_s1.py b/find_s1.py
index 7aed094..77f5104 100644
--- a/find_s1.py
+++ b/find_s1.py
@@ -1,8 +1,7 @@
#!/usr/bin/env python
import os
-import sys
-NAME = 'Find_S1'
+NAME = 'Find_S2'
AUTHORS = 'Joyce Reimer, Kyle Klenk, Raymond Spiteri'
GUIinclude = True
@@ -14,31 +13,32 @@ from carputils.carpio import txt
import numpy as np
from numpy import array as nplist
-def parser_for_s1():
+def parser_for_s2():
parser = tools.standard_parser()
- group = parser.add_argument_group('Experiment for finding S1 threshold')
+ group = parser.add_argument_group('experiment specific options')
group.add_argument('--duration',
type = float,
default = 20.,
help = 'Duration of simulation in [ms] (default: 20.)')
group.add_argument('--S1-strength',
- type = float,
- default = 20.,
- help = 'pick transmembrane current stimulus strength in [uA/cm^2] = [pA/pF] considering fixed Cm (default: 20.)')
+ type = float,
+ default = 20.,
+ help = 'pick S1 transmembrane current stimulus strength in [uA/cm^2] = [pA/pF] considering fixed Cm (default: 20.)')
group.add_argument('--S1-dur',
type = float,
default = 2.,
help = 'pick transmembrane current stimulus duration in [ms] (default: 2.)')
return parser
+
def jobID(args):
"""
Generate name of top level output directory.
"""
today = date.today()
- return 'simulations/output/s1_attempt_{}'.format(int(args.S1_strength))
+ return 'simulations/output/s1_{}_attempt_{}'.format(int(args.duration), int(args.S1_strength))
-@tools.carpexample(parser_for_s1, jobID)
+@tools.carpexample(parser_for_s2, jobID)
def run(args, job):
#===========================================
# 1: Define the mesh: Units are mm
@@ -74,23 +74,20 @@ def run(args, job):
# Define the geometry of the stimulus at one end of the block
# Units are um
stim = ['-num_stim', 1,
- '-stimulus[0].name', 'S1',
- '-stimulus[0].stimtype', 0,
- '-stimulus[0].strength', args.S1_strength,
- '-stimulus[0].duration', args.S1_dur ]
-
- electrode = mesh.block_boundary_condition(geom, 'stimulus', 0, 'x', True)
+ '-stimulus[0].name', 'S1', #index 0 because only one type of electrode
+ '-stimulus[0].stimtype', 0, #0 means intracellular stim
+ '-stimulus[0].strength', args.S1_strength, #uA/Cm^2
+ '-stimulus[0].duration', args.S1_dur]
+ #electrode = mesh.block_boundary_condition(geom, 'stimulus', 0, 'x', True) #True refers to a BC between tissue and bath
+ E1_lower_bound = nplist([ 0.0, 0.0, 0.0])
+ E1_upper_bound = nplist([ 2.0, 0.3, 0.2])
+ electrode = mesh.block_region(geom, 'stimulus', 0, E1_lower_bound, E1_upper_bound, False) # introudce a stimulus in the region defined by the lower and upper bounds
- # #electrode = mesh.block_boundary_condition(geom, 'stimulus', 0, 'x', True) #True refers to a BC between tissue and bath
- # E1_lower_bound = nplist([ 0.0, 0.0, 0.0])
- # E1_upper_bound = nplist([ 2.0, 0.3, 0.2])
- # electrode = mesh.block_region(geom, 'stimulus', 0, E1_lower_bound, E1_upper_bound, False) # introudce a stimulus in the region defined by the lower and upper bounds
-
- # #electrode = mesh.block_boundary_condition(geom, 'stimulus', 0, 'x', True) #True refers to a BC between tissue and bath
- # E1_lower_bound = nplist([ 0.0, 0.0, 0.0])
- # E1_upper_bound = nplist([ 2.0, 0.3, 0.2])
- # electrode = mesh.block_region(geom, 'stimulus', 0, E1_lower_bound, E1_upper_bound, False) # introudce a stimulus in the region defined by the lower and upper bounds
+ #electrode = mesh.block_boundary_condition(geom, 'stimulus', 0, 'x', True) #True refers to a BC between tissue and bath
+ E1_lower_bound = nplist([ 0.0, 0.0, 0.0])
+ E1_upper_bound = nplist([ 2.0, 0.3, 0.2])
+ electrode = mesh.block_region(geom, 'stimulus', 0, E1_lower_bound, E1_upper_bound, False) # introudce a stimulus in the region defined by the lower and upper bounds
#===========================================
# 4: Defining simulator options:
@@ -109,6 +106,4 @@ def run(args, job):
# Run simulation
job.carp(cmd) #cmd includes the params
-
-
-run()
+run()
\ No newline at end of file
diff --git a/find_s2.py b/find_s2.py
index c522555..5d28d79 100644
--- a/find_s2.py
+++ b/find_s2.py
@@ -13,7 +13,7 @@ from carputils.carpio import txt
import numpy as np
from numpy import array as nplist
-def parser_for_s2(s1_strength):
+def parser_for_s2():
parser = tools.standard_parser()
group = parser.add_argument_group('experiment specific options')
group.add_argument('--duration',
@@ -22,7 +22,7 @@ def parser_for_s2(s1_strength):
help = 'Duration of simulation in [ms] (default: 20.)')
group.add_argument('--S1-strength',
type = float,
- default = s1_strength,
+ default = 172,
help = 'pick S1 transmembrane current stimulus strength in [uA/cm^2] = [pA/pF] considering fixed Cm (default: 20.)')
group.add_argument('--S1-dur',
type = float,
@@ -52,11 +52,12 @@ def jobID(args):
Generate name of top level output directory.
"""
today = date.today()
- return 'output/s2_{}_basic_{}'.format(today.isoformat(), args.S2_strength)
+ return 'simulations/output/s2_{}_attempt_{}'.format(int(args.S2_start), int(args.S2_strength))
+
@tools.carpexample(parser_for_s2, jobID)
def run(args, job):
- #===========================================
+ #===========================================
# 1: Define the mesh: Units are mm
#===========================================
diff --git a/legacy_scripts/find_s1_example.py b/legacy_scripts/find_s1_example.py
new file mode 100644
index 0000000..7594bd5
--- /dev/null
+++ b/legacy_scripts/find_s1_example.py
@@ -0,0 +1,101 @@
+#!/usr/bin/env python
+
+import os
+import sys
+NAME = 'Find_S1'
+AUTHORS = 'Joyce Reimer, Kyle Klenk, Raymond Spiteri'
+GUIinclude = True
+
+from datetime import date
+from carputils import settings
+from carputils import tools
+from carputils import mesh
+from carputils.carpio import txt
+import numpy as np
+from numpy import array as nplist
+
+def parser_for_s1():
+ parser = tools.standard_parser()
+ group = parser.add_argument_group('Experiment for finding S1 threshold')
+ group.add_argument('--duration',
+ type = float,
+ default = 20.,
+ help = 'Duration of simulation in [ms] (default: 20.)')
+ group.add_argument('--S1-strength',
+ type = float,
+ default = 20.,
+ help = 'pick transmembrane current stimulus strength in [uA/cm^2] = [pA/pF] considering fixed Cm (default: 20.)')
+ group.add_argument('--S1-dur',
+ type = float,
+ default = 2.,
+ help = 'pick transmembrane current stimulus duration in [ms] (default: 2.)')
+ return parser
+
+def jobID(args):
+ """
+ Generate name of top level output directory.
+ """
+ today = date.today()
+ return 'simulations/output/s1_attempt_{}'.format(int(args.S1_strength))
+
+@tools.carpexample(parser_for_s1, jobID)
+def run(args, job):
+ #===========================================
+ # 1: Define the mesh: Units are mm
+ #===========================================
+
+ # Create a block with a thin z dimension
+ x = 10.0 #(mm)
+ y = 1.0 #(mm)
+ z = 0.2 #(mm)
+ geom = mesh.Block(size=(x,y,z))
+
+ # Set fibre angle to 0, sheet angle to 0
+ geom.set_fibres(0, 0, 90, 90) #endo, epi, sheet endo, sheet epicard
+
+ # Finilaize the creation of the mesh
+ meshname = mesh.generate(geom)
+
+ #===========================================
+ # 2: Define the ionic models & conductivities
+ #===========================================
+
+ # Query for element labels
+ _, etags, _ = txt.read(meshname + '.elem')
+ etags = np.unique(etags)
+
+ IntraTags = etags[etags != 0].tolist()
+ ExtraTags = etags.tolist().copy()
+
+ #===========================================
+ # 3: Configure the simulation
+ #===========================================
+
+ # Define the geometry of the stimulus at one end of the block
+ # Units are um
+ stim = ['-num_stim', 1,
+ '-stimulus[0].name', 'S1',
+ '-stimulus[0].stimtype', 0,
+ '-stimulus[0].strength', args.S1_strength,
+ '-stimulus[0].duration', args.S1_dur ]
+
+ electrode = mesh.block_boundary_condition(geom, 'stimulus', 0, 'x', True)
+
+ #===========================================
+ # 4: Defining simulator options:
+ #===========================================
+
+ # Get basic command line, including solver options
+ cmd = tools.carp_cmd(os.path.join('simulations/input/basic.par')) #cmd is a carp command
+
+ cmd += ['-simID', job.ID,
+ '-meshname', meshname,
+ '-dt', 25,
+ '-tend', args.duration]
+ cmd += tools.gen_physics_opts(ExtraTags=ExtraTags, IntraTags=IntraTags)
+ cmd += stim + electrode
+
+ # Run simulation
+ job.carp(cmd) #cmd includes the params
+
+run()
diff --git a/main.cpp b/main.cpp
index 50d264b..6a3e784 100644
--- a/main.cpp
+++ b/main.cpp
@@ -17,10 +17,15 @@ using namespace caf;
class config : public actor_system_config {
public:
int duration = 20;
- int s1_min;
- int s1_max;
+ int s1_min = 0;
+ int s1_max = 100;
int tol = 1;
int s1_duration = 1;
+ int s2_min = 153000;
+ int s2_max = 157000;
+ int s2_interval_min = 220;
+ int s2_interval_max = 220;
+ int s2_interval_increment = 5;
config() {
opt_group{custom_options_, "global"}
@@ -28,7 +33,12 @@ class config : public actor_system_config {
.add(s1_min, "S1-min", "minimum S1 strength for bisection [uA/cm^2]")
.add(s1_max, "S1-max", "maximum S1 strength for bisection [uA/cm^2]")
.add(s1_duration, "S1-dur", "Duration of S1 stimulus in [ms]")
- .add(tol, "tol", "Tolerance of Bisection Method");
+ .add(tol, "tol", "Tolerance of Bisection Method")
+ .add(s2_min, "S2-min", "minimum S2 strength for bisection [uA/cm^2]")
+ .add(s2_max, "S2-max", "maximum S2 strength for bisection [uA/cm^2]")
+ .add(s2_interval_min, "S2-interval-min", "minimum S2 interval for bisection [ms]")
+ .add(s2_interval_max, "S2-interval-max", "maximum S2 interval for bisection [ms]")
+ .add(s2_interval_increment, "S2-interval-inc", "increment of S2 interval for bisection [ms]");
}
};
@@ -54,7 +64,7 @@ std::string execute_command(std::stringstream &command) {
-bool verify_ap(int s1_strength) {
+bool verify_ap(int s1_strength, int start_index, bool find_s2=false) {
int num_nodes = 0;
bool ap_occured = false;
std::stringstream results_file_path, igbhead_command, igbextract_command,
@@ -62,8 +72,17 @@ bool verify_ap(int s1_strength) {
std::string command_output;
std::smatch match; // for regex matchin
+ // Set File Names
+ if (find_s2) {
+ results_file_path << "simulations/output/s2_" << start_index << "_attempt_" << s1_strength << "/vm.igb";
+ tmp_output_file_path << "simulations/output/s2_" << start_index << "_attempt_" << s1_strength << "/tmp/";
+ } else {
+ results_file_path << "simulations/output/s1_attempt_" << s1_strength << "/vm.igb";
+ tmp_output_file_path << "simulations/output/s1_attempt_" << s1_strength << "/tmp/";
+ }
+
+
// Verify the file exists
- results_file_path << "simulations/output/s1_attempt_" << s1_strength << "/vm.igb";
std::ifstream file(results_file_path.str().c_str());
if (!file.good()) throw std::runtime_error("Failed to find the results file");
@@ -78,7 +97,6 @@ bool verify_ap(int s1_strength) {
throw std::runtime_error("Failed to find the number of nodes in the mesh");
}
- tmp_output_file_path << "simulations/output/s1_attempt_" << s1_strength << "/tmp/";
if (!std::filesystem::create_directories(tmp_output_file_path.str().c_str()))
throw std::runtime_error("Failed to create temp output directory");
@@ -100,7 +118,7 @@ bool verify_ap(int s1_strength) {
tmp_output_file_in.close();
// Check if AP occured (if there is any value within the array larger than 0.0)
- if (std::any_of(node_data.begin(), node_data.end(), [](float i){return i > 0.0;})) {
+ if (std::any_of(node_data.begin() + start_index, node_data.end(), [](float i){return i > 0.0;})) {
ap_occured = true;
} else {
ap_occured = false;
@@ -123,6 +141,117 @@ int bisection(int s1_min, int s1_max) {
return s1_mid;
}
+behavior s2_actor(event_based_actor* self) {
+ return {
+ [=](int s1_strenth, int s2_min, int s2_max, int interval_min,
+ int interval_max, int interval_increment, int tol) {
+
+ int interval = interval_max;
+ int duration = 300;
+ int current_guess;
+ int s2_min_orginal = s2_min;
+ int s2_max_orginal = s2_max;
+ std::stringstream python_command, log_file;
+ std::string command_output = "";
+
+
+ // Verify S2_min does not cause an AP
+ aout(self) << "Verifying S2_min does not cause an AP\n";
+
+ python_command << "python3 /code/find_s2.py"
+ << " --S1-strength " << s1_strenth
+ << " --S2-strength " << s2_min
+ << " --S2-start " << interval
+ << " --duration " << duration;
+
+ log_file << "simulations/logs/output_s2_" << s2_min << ".txt";
+
+ command_output = execute_command(python_command);
+ if (verify_ap(s2_min, interval, true))
+ throw std::runtime_error("AP occured at S1_min - Select a new S1_min");
+
+ python_command.str("");
+ log_file.str("");
+
+ // Verify S2_max causes an AP
+ aout(self) << "Verifying S2_max causes an AP\n";
+
+ python_command << "python3 /code/find_s2.py"
+ << " --S1-strength " << s1_strenth
+ << " --S2-strength " << s2_max
+ << " --S2-start " << interval
+ << " --duration " << duration;
+
+ log_file << "simulations/logs/output_s2_" << s2_max << ".txt";
+
+ command_output = execute_command(python_command);
+ if (!verify_ap(s2_max, interval, true))
+ throw std::runtime_error("AP did not occur at S1_max - Select a new S1_max");
+
+ python_command.str("");
+ log_file.str("");
+
+ aout(self) << "Moving to bisection\n";
+
+
+
+ while((interval - interval_min) >= 0) {
+
+ current_guess = bisection(s2_min, s2_max);
+
+
+
+ // Inner threshold bisection loop
+ while(s2_max - s2_min > tol) {
+ aout(self) << "Current guess: " << current_guess << "\n";
+ aout(self) << "S2_min: " << s2_min << "\n";
+ aout(self) << "S2_max: " << s2_max << "\n";
+ aout(self) << "Interval: " << interval << "\n";
+
+ // Verify the current guess does not cause an AP
+ aout(self) << "Verifying current guess does not cause an AP\n";
+
+ python_command << "python3 /code/find_s2.py"
+ << " --S1-strength " << s1_strenth
+ << " --S2-strength " << current_guess
+ << " --S2-start " << interval
+ << " --duration " << duration;
+
+ log_file << "simulations/logs/output_s2_" << current_guess << ".txt";
+
+ command_output = execute_command(python_command);
+ if (verify_ap(current_guess, interval, true)) {
+ s2_max = current_guess;
+ } else {
+ s2_min = current_guess;
+ }
+
+ python_command.str("");
+ log_file.str("");
+
+ current_guess = bisection(s2_min, s2_max);
+ }
+
+ aout(self) << "S2 threshold found: " << current_guess << "\n";
+
+ aout(self) << "Adjusting interval\n";
+
+ interval = interval - interval_increment;
+
+ aout(self) << "New interval: " << interval << "\n";
+
+ // Reset the S2_min and S2_max values
+ s2_min = s2_min_orginal;
+ s2_max = s2_max_orginal;
+
+
+ }
+
+
+ }
+ };
+}
+
/**
* Actor that finds the S1 threshold using bisection
@@ -153,7 +282,7 @@ behavior s1_actor(event_based_actor* self) {
out_stream << command_output;
out_stream.close();
- if (verify_ap(s1_strength_min))
+ if (verify_ap(s1_strength_min, 0))
throw std::runtime_error("AP occured at S1_min - Select a new S1_min");
python_command.str("");
@@ -174,7 +303,7 @@ behavior s1_actor(event_based_actor* self) {
out_stream << command_output;
out_stream.close();
- if (!verify_ap(s1_strength_max))
+ if (!verify_ap(s1_strength_max, 0))
throw std::runtime_error("AP did not occur at S1_max - Select a new S1_max");
python_command.str("");
@@ -186,7 +315,7 @@ behavior s1_actor(event_based_actor* self) {
current_guess = bisection(s1_strength_min, s1_strength_max);
while(s1_strength_max - s1_strength_min > tol) {
- aout(self) << "Current guess: " << current_guess << "\n";
+ aout(self) << "\nNew guess: " << current_guess << "\n";
python_command << "python3 /code/find_s1.py"
<< " --duration " << duration
@@ -200,11 +329,11 @@ behavior s1_actor(event_based_actor* self) {
out_stream << command_output;
out_stream.close();
- if (verify_ap(current_guess)) {
- aout(self) << "AP occured at current guess\n";
+ if (verify_ap(current_guess, 0)) {
+ aout(self) << "AP occured at: " << current_guess << "\n";
s1_strength_max = current_guess;
} else {
- aout(self) << "AP did not occur at current guess\n";
+ aout(self) << "AP did not occur at: " << current_guess << "\n";
s1_strength_min = current_guess;
}
@@ -214,7 +343,8 @@ behavior s1_actor(event_based_actor* self) {
command_output = "";
}
- aout(self) << "Final S1_strength: " << s1_strength_max << "\n";
+ // aout(self) << "Final S1_strength: " << s1_strength_max << "\n";
+ return s1_strength_max;
}
};
@@ -232,13 +362,23 @@ void caf_main(actor_system& system, const config& cfg) {
aout(self) << "S1 min: " << cfg.s1_min << "\n";
aout(self) << "S1 max: " << cfg.s1_max << "\n";
-
// Cleanup the output directory
- for (auto& p: std::filesystem::directory_iterator("simulations/output/")) {
+ for (auto& p: std::filesystem::directory_iterator("simulations/output/"))
std::filesystem::remove_all(p);
- }
+
auto s1_solver = system.spawn(s1_actor);
+ self->request(s1_solver, infinite, cfg.duration, cfg.s1_min, cfg.s1_max, cfg.tol).receive(
+ [&](int s1_strength) {
+ aout(self) << "S1 strength: " << s1_strength << "\n";
+ auto s2_solver = system.spawn(s2_actor);
+ self->send(s2_solver, s1_strength, cfg.s2_min, cfg.s2_max, cfg.s2_interval_min,
+ cfg.s2_interval_max, cfg.s2_interval_increment, cfg.tol);
+ },
+ [&](error& err) {
+ aout(self) << "Error: " << err << "\n";
+ }
+ );
self->send(s1_solver, cfg.duration, cfg.s1_min, cfg.s1_max, cfg.tol);
}
--
GitLab