From 9dc7a50073fb116cc161ccbc699a6d3fb866a15e Mon Sep 17 00:00:00 2001
From: Kyle <kyle.c.klenk@gmail.com>
Date: Wed, 28 Jun 2023 20:42:48 +0000
Subject: [PATCH] Actor can use bisection to obtain the correct result for S1

---
 .gitignore |   1 +
 main.cpp   | 169 ++++++++++++++++++++++++++++++++++++++++++-----------
 2 files changed, 136 insertions(+), 34 deletions(-)

diff --git a/.gitignore b/.gitignore
index f7c4f67..8b90690 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 d002843..af06b5d 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);
 
 }
 
-- 
GitLab