diff --git a/.gitignore b/.gitignore
index f7c4f67eaddd7e080cf84294eecbcc33e4cdeefa..8b9069068262e589c426c5e9623f95641efcbd33 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,3 +5,4 @@ build/
meshes/
output8.txt
.vscode/*
+simulations/logs/*
diff --git a/main.cpp b/main.cpp
index d00284313bfec72bdc6f7308d46856d127b3140c..af06b5d5aeedc73b95ff74465d090ed93c2981b9 100644
--- a/main.cpp
+++ b/main.cpp
@@ -9,15 +9,32 @@
#include <regex>
#include <filesystem>
#include <algorithm>
+#include <stdexcept>
+#include <cmath>
using namespace caf;
-// float bisection(float val_1, float val_2) {
+class config : public actor_system_config {
+ public:
+ int duration = 20;
+ int s1_min;
+ int s1_max;
+ int tol = 1;
+ int s1_duration = 1;
+
+ config() {
+ opt_group{custom_options_, "global"}
+ .add(duration, "duration", "Duration of simulation in [ms]")
+ .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");
+ }
+};
+
-// return
-// }
-bool verify_ap(float s1_strength) {
+bool verify_ap(int s1_strength) {
int num_nodes = 0;
bool ap_occured = false;
std::stringstream results_file_path, igbhead_command, igbextract_command,
@@ -101,41 +118,124 @@ bool verify_ap(float s1_strength) {
return ap_occured;
}
+int bisection(int s1_min, int s1_max) {
+ int s1_mid = (s1_max - s1_min) / 2.0 + s1_min;
+ return s1_mid;
+}
+
+
+
behavior s1_actor(event_based_actor* self) {
return {
- [=](float duration, float s1_strength, float s1_duration) {
+ [=](int duration, int s1_strength_min, int s1_strength_max, int tol) {
+ // Initialize the variables
+ std::stringstream python_command, log_file;
+ FILE *pipe;
+ char buffer[128];
+ std::string result = "";
+ int current_guess;
+
+ // Verify S1_min does not cause an AP
+ aout(self) << "Verifying S1_min does not cause an AP\n";
- // Run the S1 Script
- std::stringstream python_command;
python_command << "python3 /code/find_s1.py"
<< " --duration " << duration
- << " --S1-strength " << s1_strength
- << " --S1-dur " << s1_duration;
+ << " --S1-strength " << s1_strength_min;
- // create a pipe and run the python command
- FILE* pipe = popen(python_command.str().c_str(), "r");
-
- if (!pipe) std::cerr << "popen() failed!" << std::endl;
- char buffer[128];
- std::string result = "";
+ log_file << "simulations/logs/output" << s1_strength_min << ".txt";
+
+ pipe = popen(python_command.str().c_str(), "r");
+ if (!pipe) throw std::runtime_error("popen() failed!");
while (!feof(pipe)) {
if (fgets(buffer, 128, pipe) != NULL)
result += buffer;
}
- pclose(pipe); // close the pipe
-
- // Log output to a txt file
- std::stringstream log_file;
- log_file << "/code/simulations/output/logs/output" << s1_strength << ".txt";
+ pclose(pipe);
std::ofstream out(log_file.str().c_str());
out << result;
out.close();
- // Check if AP occured
- bool ap_occured = verify_ap(s1_strength);
- aout(self) << "AP occured: " << ap_occured << std::endl;
+ if (verify_ap(s1_strength_min)) {
+ aout(self) << "AP occured at S1_min\n";
+ throw std::runtime_error("AP occured at S1_min - Select a new S1_min");
+ }
+
+ python_command.str("");
+ log_file.str("");
+
+ // Verify S1_max causes an AP
+ aout(self) << "Verifying S1_max causes an AP\n";
+
+ python_command << "python3 /code/find_s1.py"
+ << " --duration " << duration
+ << " --S1-strength " << s1_strength_max;
+
+ log_file << "simulations/logs/output" << s1_strength_max << ".txt";
+
+ pipe = popen(python_command.str().c_str(), "r");
+ if (!pipe) throw std::runtime_error("popen() failed!");
+ while (!feof(pipe)) {
+ if (fgets(buffer, 128, pipe) != NULL)
+ result += buffer;
+ }
+ pclose(pipe);
+
+ out.open(log_file.str().c_str());
+ out << result;
+ out.close();
+
+ if (!verify_ap(s1_strength_max)) {
+ aout(self) << "AP did not occur at S1_max\n";
+ 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 on to bisection\n";
+
+ // Find new S1_strength
+ current_guess = bisection(s1_strength_min, s1_strength_max);
+
+ while(s1_strength_max - s1_strength_min > tol) {
+ aout(self) << "Current guess: " << current_guess << "\n";
+
+ python_command << "python3 /code/find_s1.py"
+ << " --duration " << duration
+ << " --S1-strength " << current_guess;
+
+ log_file << "simulations/logs/output" << current_guess << ".txt";
+
+ pipe = popen(python_command.str().c_str(), "r");
+ if (!pipe) throw std::runtime_error("popen() failed!");
+ while (!feof(pipe)) {
+ if (fgets(buffer, 128, pipe) != NULL)
+ result += buffer;
+ }
+ pclose(pipe);
+
+ out.open(log_file.str().c_str());
+ out << result;
+ out.close();
+
+ if (verify_ap(current_guess)) {
+ aout(self) << "AP occured at current guess\n";
+ s1_strength_max = current_guess;
+ } else {
+ aout(self) << "AP did not occur at current guess\n";
+ s1_strength_min = current_guess;
+ }
+
+ current_guess = bisection(s1_strength_min, s1_strength_max);
+ python_command.str("");
+ log_file.str("");
+ result = "";
+ }
+
+ aout(self) << "Final S1_strength: " << s1_strength_max << "\n";
+
}
};
@@ -144,22 +244,23 @@ behavior s1_actor(event_based_actor* self) {
-void caf_main(actor_system& system) {
+void caf_main(actor_system& system, const config& cfg) {
scoped_actor self{system};
+ // Check command line arguments
+ aout(self) << "Command Line Arguments: \n";
+ aout(self) << "Duration: " << cfg.duration << "\n";
+ aout(self) << "S1 min: " << cfg.s1_min << "\n";
+ aout(self) << "S1 max: " << cfg.s1_max << "\n";
- auto actor_1 = self->spawn(s1_actor);
- auto actor_2 = self->spawn(s1_actor);
-
- float duration = 20.0;
- float s1_strength_min = 100.0;
- float s1_duration = 2.0;
- // float s1_strength_max = 200.0;
-
+ // Cleanup the output directory
+ for (auto& p: std::filesystem::directory_iterator("simulations/output/")) {
+ std::filesystem::remove_all(p);
+ }
- self->send(actor_1, duration, s1_strength_min, s1_duration);
- // self->send(actor_2, duration, s1_strength_max, s1_duration);
+ auto s1_solver = system.spawn(s1_actor);
+ self->send(s1_solver, cfg.duration, cfg.s1_min, cfg.s1_max, cfg.tol);
}