Skip to content
Snippets Groups Projects
Commit 8f4049f5 authored by Kyle's avatar Kyle
Browse files

Add Actor to find s2

parent 4eebac05
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
.gitignore
+ 1
− 0
View file @ 8f4049f5
...@@ -6,3 +6,4 @@ meshes/ ...@@ -6,3 +6,4 @@ meshes/
output8.txt output8.txt
.vscode/* .vscode/*
simulations/logs/* simulations/logs/*
.DS_Store
environment_setup/Dockerfile
+ 1
− 0
View file @ 8f4049f5
...@@ -70,6 +70,7 @@ RUN cmake --build _build ...@@ -70,6 +70,7 @@ RUN cmake --build _build
WORKDIR /opt/openCARP/external/carputils WORKDIR /opt/openCARP/external/carputils
RUN python3 -m pip install --user -r requirements.txt RUN python3 -m pip install --user -r requirements.txt
ENV PATH=$PATH:/opt/openCARP/external/carputils/bin ENV PATH=$PATH:/opt/openCARP/external/carputils/bin
ENV PATH=$PATH:/opt/openCARP/_build/bin
ENV PYTHONPATH=$PYTHONPATH:/opt/openCARP/external/carputils/ ENV PYTHONPATH=$PYTHONPATH:/opt/openCARP/external/carputils/
# Create a settings.yaml - add your own settings.yaml # Create a settings.yaml - add your own settings.yaml
......
find_s1.py
+ 22
− 27
View file @ 8f4049f5
#!/usr/bin/env python #!/usr/bin/env python
import os import os
import sys NAME = 'Find_S2'
NAME = 'Find_S1'
AUTHORS = 'Joyce Reimer, Kyle Klenk, Raymond Spiteri' AUTHORS = 'Joyce Reimer, Kyle Klenk, Raymond Spiteri'
GUIinclude = True GUIinclude = True
...@@ -14,31 +13,32 @@ from carputils.carpio import txt ...@@ -14,31 +13,32 @@ from carputils.carpio import txt
import numpy as np import numpy as np
from numpy import array as nplist from numpy import array as nplist
def parser_for_s1(): def parser_for_s2():
parser = tools.standard_parser() 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', group.add_argument('--duration',
type = float, type = float,
default = 20., default = 20.,
help = 'Duration of simulation in [ms] (default: 20.)') help = 'Duration of simulation in [ms] (default: 20.)')
group.add_argument('--S1-strength', group.add_argument('--S1-strength',
type = float, type = float,
default = 20., default = 20.,
help = 'pick transmembrane current stimulus strength in [uA/cm^2] = [pA/pF] considering fixed Cm (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', group.add_argument('--S1-dur',
type = float, type = float,
default = 2., default = 2.,
help = 'pick transmembrane current stimulus duration in [ms] (default: 2.)') help = 'pick transmembrane current stimulus duration in [ms] (default: 2.)')
return parser return parser
def jobID(args): def jobID(args):
""" """
Generate name of top level output directory. Generate name of top level output directory.
""" """
today = date.today() 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): def run(args, job):
#=========================================== #===========================================
# 1: Define the mesh: Units are mm # 1: Define the mesh: Units are mm
...@@ -74,23 +74,20 @@ def run(args, job): ...@@ -74,23 +74,20 @@ def run(args, job):
# Define the geometry of the stimulus at one end of the block # Define the geometry of the stimulus at one end of the block
# Units are um # Units are um
stim = ['-num_stim', 1, stim = ['-num_stim', 1,
'-stimulus[0].name', 'S1', '-stimulus[0].name', 'S1', #index 0 because only one type of electrode
'-stimulus[0].stimtype', 0, '-stimulus[0].stimtype', 0, #0 means intracellular stim
'-stimulus[0].strength', args.S1_strength, '-stimulus[0].strength', args.S1_strength, #uA/Cm^2
'-stimulus[0].duration', args.S1_dur ] '-stimulus[0].duration', args.S1_dur]
electrode = mesh.block_boundary_condition(geom, 'stimulus', 0, 'x', True)
#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 #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_lower_bound = nplist([ 0.0, 0.0, 0.0])
# E1_upper_bound = nplist([ 2.0, 0.3, 0.2]) 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_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: # 4: Defining simulator options:
...@@ -109,6 +106,4 @@ def run(args, job): ...@@ -109,6 +106,4 @@ def run(args, job):
# Run simulation # Run simulation
job.carp(cmd) #cmd includes the params job.carp(cmd) #cmd includes the params
run()
\ No newline at end of file
run()
find_s2.py
+ 5
− 4
View file @ 8f4049f5
...@@ -13,7 +13,7 @@ from carputils.carpio import txt ...@@ -13,7 +13,7 @@ from carputils.carpio import txt
import numpy as np import numpy as np
from numpy import array as nplist from numpy import array as nplist
def parser_for_s2(s1_strength): def parser_for_s2():
parser = tools.standard_parser() parser = tools.standard_parser()
group = parser.add_argument_group('experiment specific options') group = parser.add_argument_group('experiment specific options')
group.add_argument('--duration', group.add_argument('--duration',
...@@ -22,7 +22,7 @@ def parser_for_s2(s1_strength): ...@@ -22,7 +22,7 @@ def parser_for_s2(s1_strength):
help = 'Duration of simulation in [ms] (default: 20.)') help = 'Duration of simulation in [ms] (default: 20.)')
group.add_argument('--S1-strength', group.add_argument('--S1-strength',
type = float, 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.)') help = 'pick S1 transmembrane current stimulus strength in [uA/cm^2] = [pA/pF] considering fixed Cm (default: 20.)')
group.add_argument('--S1-dur', group.add_argument('--S1-dur',
type = float, type = float,
...@@ -52,11 +52,12 @@ def jobID(args): ...@@ -52,11 +52,12 @@ def jobID(args):
Generate name of top level output directory. Generate name of top level output directory.
""" """
today = date.today() 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) @tools.carpexample(parser_for_s2, jobID)
def run(args, job): def run(args, job):
#=========================================== #===========================================
# 1: Define the mesh: Units are mm # 1: Define the mesh: Units are mm
#=========================================== #===========================================
......
legacy_scripts/find_s1_example.py 0 → 100644
+ 101
− 0
View file @ 8f4049f5
#!/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()
main.cpp
+ 157
− 17
View file @ 8f4049f5
...@@ -17,10 +17,15 @@ using namespace caf; ...@@ -17,10 +17,15 @@ using namespace caf;
class config : public actor_system_config { class config : public actor_system_config {
public: public:
int duration = 20; int duration = 20;
int s1_min; int s1_min = 0;
int s1_max; int s1_max = 100;
int tol = 1; int tol = 1;
int s1_duration = 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() { config() {
opt_group{custom_options_, "global"} opt_group{custom_options_, "global"}
...@@ -28,7 +33,12 @@ class config : public actor_system_config { ...@@ -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_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_max, "S1-max", "maximum S1 strength for bisection [uA/cm^2]")
.add(s1_duration, "S1-dur", "Duration of S1 stimulus in [ms]") .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) { ...@@ -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; int num_nodes = 0;
bool ap_occured = false; bool ap_occured = false;
std::stringstream results_file_path, igbhead_command, igbextract_command, std::stringstream results_file_path, igbhead_command, igbextract_command,
...@@ -62,8 +72,17 @@ bool verify_ap(int s1_strength) { ...@@ -62,8 +72,17 @@ bool verify_ap(int s1_strength) {
std::string command_output; std::string command_output;
std::smatch match; // for regex matchin 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 // Verify the file exists
results_file_path << "simulations/output/s1_attempt_" << s1_strength << "/vm.igb";
std::ifstream file(results_file_path.str().c_str()); std::ifstream file(results_file_path.str().c_str());
if (!file.good()) throw std::runtime_error("Failed to find the results file"); if (!file.good()) throw std::runtime_error("Failed to find the results file");
...@@ -78,7 +97,6 @@ bool verify_ap(int s1_strength) { ...@@ -78,7 +97,6 @@ bool verify_ap(int s1_strength) {
throw std::runtime_error("Failed to find the number of nodes in the mesh"); 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())) if (!std::filesystem::create_directories(tmp_output_file_path.str().c_str()))
throw std::runtime_error("Failed to create temp output directory"); throw std::runtime_error("Failed to create temp output directory");
...@@ -100,7 +118,7 @@ bool verify_ap(int s1_strength) { ...@@ -100,7 +118,7 @@ bool verify_ap(int s1_strength) {
tmp_output_file_in.close(); tmp_output_file_in.close();
// Check if AP occured (if there is any value within the array larger than 0.0) // 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; ap_occured = true;
} else { } else {
ap_occured = false; ap_occured = false;
...@@ -123,6 +141,117 @@ int bisection(int s1_min, int s1_max) { ...@@ -123,6 +141,117 @@ int bisection(int s1_min, int s1_max) {
return s1_mid; 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 * Actor that finds the S1 threshold using bisection
...@@ -153,7 +282,7 @@ behavior s1_actor(event_based_actor* self) { ...@@ -153,7 +282,7 @@ behavior s1_actor(event_based_actor* self) {
out_stream << command_output; out_stream << command_output;
out_stream.close(); 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"); throw std::runtime_error("AP occured at S1_min - Select a new S1_min");
python_command.str(""); python_command.str("");
...@@ -174,7 +303,7 @@ behavior s1_actor(event_based_actor* self) { ...@@ -174,7 +303,7 @@ behavior s1_actor(event_based_actor* self) {
out_stream << command_output; out_stream << command_output;
out_stream.close(); 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"); throw std::runtime_error("AP did not occur at S1_max - Select a new S1_max");
python_command.str(""); python_command.str("");
...@@ -186,7 +315,7 @@ behavior s1_actor(event_based_actor* self) { ...@@ -186,7 +315,7 @@ behavior s1_actor(event_based_actor* self) {
current_guess = bisection(s1_strength_min, s1_strength_max); current_guess = bisection(s1_strength_min, s1_strength_max);
while(s1_strength_max - s1_strength_min > tol) { 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" python_command << "python3 /code/find_s1.py"
<< " --duration " << duration << " --duration " << duration
...@@ -200,11 +329,11 @@ behavior s1_actor(event_based_actor* self) { ...@@ -200,11 +329,11 @@ behavior s1_actor(event_based_actor* self) {
out_stream << command_output; out_stream << command_output;
out_stream.close(); out_stream.close();
if (verify_ap(current_guess)) { if (verify_ap(current_guess, 0)) {
aout(self) << "AP occured at current guess\n"; aout(self) << "AP occured at: " << current_guess << "\n";
s1_strength_max = current_guess; s1_strength_max = current_guess;
} else { } 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; s1_strength_min = current_guess;
} }
...@@ -214,7 +343,8 @@ behavior s1_actor(event_based_actor* self) { ...@@ -214,7 +343,8 @@ behavior s1_actor(event_based_actor* self) {
command_output = ""; 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) { ...@@ -232,13 +362,23 @@ void caf_main(actor_system& system, const config& cfg) {
aout(self) << "S1 min: " << cfg.s1_min << "\n"; aout(self) << "S1 min: " << cfg.s1_min << "\n";
aout(self) << "S1 max: " << cfg.s1_max << "\n"; aout(self) << "S1 max: " << cfg.s1_max << "\n";
// Cleanup the output directory // 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); std::filesystem::remove_all(p);
}
auto s1_solver = system.spawn(s1_actor); 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); self->send(s1_solver, cfg.duration, cfg.s1_min, cfg.s1_max, cfg.tol);
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment