main.cpp 8.08 KiB
#include <iostream>
#include <fstream>
#include <chrono>
#include <thread>
#include <caf/all.hpp>
#include "caf/io/all.hpp"
#include <cstdio>
#include <sstream>
#include <regex>
#include <filesystem>
#include <algorithm>
#include <stdexcept>
#include <cmath>
using namespace caf;
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");
}
};
/**
* Function that takes in a command and
* returns the piped output as a string
*/
std::string execute_command(std::stringstream &command) {
std::string command_output = "";
char buffer[128];
FILE *pipe;
pipe = popen(command.str().c_str(), "r");
if (!pipe) throw std::runtime_error("popen() failed!");
while (!feof(pipe)) {
if (fgets(buffer, 128, pipe) != NULL)
command_output += buffer;
}
pclose(pipe);
return command_output;
}
bool verify_ap(int s1_strength) {
int num_nodes = 0;
bool ap_occured = false;
std::stringstream results_file_path, igbhead_command, igbextract_command,
tmp_output_file_path, tmp_output_file;
std::string command_output;
std::smatch match; // for regex matchin
// 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");
// Find the number of nodes in the mesh igbhead_command << "igbhead " << results_file_path.str().c_str();
command_output = execute_command(igbhead_command);
std::regex rgx("x dimension:\\s+(\\d+)");
if (std::regex_search(command_output, match, rgx)) {
num_nodes = std::stoi(match[1]);
} else {
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");
for (int i = 0; i < num_nodes; i++) {
tmp_output_file << tmp_output_file_path.str() << "node" << i << ".bin";
igbextract_command << "igbextract " << "-o binary -l " << i << " " <<
results_file_path.str() << " -O " << tmp_output_file.str();
command_output = execute_command(igbextract_command);
// Read in the values from the tmp_output_file
std::vector<float> node_data;
std::ifstream tmp_output_file_in(tmp_output_file.str().c_str(), std::ios::binary);
tmp_output_file_in.seekg(0, std::ios::end);
std::streampos tmp_output_file_size = tmp_output_file_in.tellg();
tmp_output_file_in.seekg(0, std::ios::beg);
node_data.resize(tmp_output_file_size / sizeof(float));
tmp_output_file_in.read(reinterpret_cast<char*>(node_data.data()), tmp_output_file_size);
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;})) {
ap_occured = true;
} else {
ap_occured = false;
break;
}
igbextract_command.str("");
tmp_output_file.str("");
}
return ap_occured;
}
/**
* Function that takes in the min and max values of S1 strength
* and returns the midpoint of the two values
*/
int bisection(int s1_min, int s1_max) {
int s1_mid = (s1_max - s1_min) / 2.0 + s1_min;
return s1_mid;
}
/**
* Actor that finds the S1 threshold using bisection
*/
behavior s1_actor(event_based_actor* self) {
return {
[=](int duration, int s1_strength_min, int s1_strength_max, int tol) {
// Initialize the variables
std::stringstream python_command, log_file;
std::string command_output = "";
int current_guess;
std::ofstream out_stream;
// Verify S1_min does not cause an AP
aout(self) << "Verifying S1_min does not cause an AP\n";
python_command << "python3 /code/find_s1.py"
<< " --duration " << duration
<< " --S1-strength " << s1_strength_min;
log_file << "simulations/logs/output" << s1_strength_min << ".txt";
// Run S1 By calling the python script
command_output = execute_command(python_command);
// Log Ouput
out_stream.open(log_file.str().c_str());
out_stream << command_output;
out_stream.close();
if (verify_ap(s1_strength_min))
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";
command_output = execute_command(python_command);
out_stream.open(log_file.str().c_str());
out_stream << command_output;
out_stream.close();
if (!verify_ap(s1_strength_max))
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";
command_output = execute_command(python_command);
out_stream.open(log_file.str().c_str());
out_stream << command_output;
out_stream.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("");
command_output = "";
}
aout(self) << "Final S1_strength: " << s1_strength_max << "\n";
}
};
}
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";
// Cleanup the output directory
for (auto& p: std::filesystem::directory_iterator("simulations/output/")) {
std::filesystem::remove_all(p);
}
auto s1_solver = system.spawn(s1_actor);
self->send(s1_solver, cfg.duration, cfg.s1_min, cfg.s1_max, cfg.tol);
}
CAF_MAIN(io::middleman)