diff --git a/Makefile b/Makefile index 1b6a71e..b93894e 100644 --- a/Makefile +++ b/Makefile @@ -1,5 +1,6 @@ CXX = g++ CXXFLAGS = -std=c++17 -Wall --coverage +LDFLAGS = -lyaml-cpp SRC_DIR = src INC_DIR = include OBJ_DIR = obj @@ -22,7 +23,7 @@ INC_DIRS = -I $(INC_DIR) # Linking step for src files $(BIN_DIR)/$(TARGET): $(OBJS) - $(CXX) $(CXXFLAGS) -o $@ $^ + $(CXX) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) # Compiling step for src files $(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp @@ -30,7 +31,7 @@ $(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp # Linking step for test files $(BIN_DIR)/$(TEST_TARGET): $(filter-out $(OBJ_DIR)/simulation.o, $(OBJS)) $(TEST_OBJS) - $(CXX) $(CXXFLAGS) -o $@ $^ + $(CXX) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) # Compiling step for test files $(TEST_OBJ_DIR)/%.o: $(TEST_DIR)/%.cpp diff --git a/configs/default.yaml b/configs/default.yaml new file mode 100644 index 0000000..f8bb1a4 --- /dev/null +++ b/configs/default.yaml @@ -0,0 +1,29 @@ +global: + nParticles: 1000 +mass: + dist: constant + val: 100000000 +x: + dist: uniform + min: 0 + max: 4 +y: + dist: uniform + min: -2 + max: 4 +z: + dist: uniform + min: 1 + max: 10 +vx: + dist: uniform + min: 0 + max: 4 +vy: + dist: normal + mu: -2 + sigma: 4 +vz: + dist: uniform + min: -10 + max: 10 diff --git a/include/environment.h b/include/environment.h index 120d017..08b4a0c 100644 --- a/include/environment.h +++ b/include/environment.h @@ -1,6 +1,7 @@ #pragma once #include #include +#include #include // Include for smart pointers #include "particle.h" @@ -10,8 +11,11 @@ class GravitationalEnvironment { // Constructors GravitationalEnvironment(const std::vector>& particlePtrs, const bool log); GravitationalEnvironment(const std::vector>& particlePtrs, const bool log, std::string logFilePrefix); + GravitationalEnvironment(const std::string configFileName, const bool log); + GravitationalEnvironment(const std::string configFileName, const bool log, std::string logFilePrefix); // Define member functions + void loadParticlesFromConfig(std::string configFileName); std::vector> getForces(const double timestep); void updateAll(const std::vector>& forces, const double timestep); void step(const double timestep); @@ -31,3 +35,4 @@ class GravitationalEnvironment { // Helper functions int getLargestLabelNumber(const std::vector& filenames, const std::string logFilePrefix); +std::map> loadConfig(const std::string& fileName); diff --git a/include/statistics.h b/include/statistics.h new file mode 100644 index 0000000..b9a5591 --- /dev/null +++ b/include/statistics.h @@ -0,0 +1,18 @@ +#pragma once +#include +#include + +// Assuming you have the random device and generator declared somewhere +extern std::random_device rd; +extern std::mt19937 GENERATOR; + +// Distribution template used to sample from random distributions +template +std::vector sampleFromDistribution(size_t n, Distribution& distribution) { + // Generate samples + std::vector samples(n); + for (size_t i = 0; i < n; ++i) { + samples[i] = distribution(GENERATOR); + } + return samples; +} diff --git a/src/environment.cpp b/src/environment.cpp index 049354d..aa1d170 100644 --- a/src/environment.cpp +++ b/src/environment.cpp @@ -1,13 +1,68 @@ #include +#include #include #include +#include +#include #include #include #include +#include #include "../include/environment.h" +#include "../include/statistics.h" namespace fs = std::filesystem; double G = 6.6743e-11; +std::string REPOPATH = std::string(std::getenv("HOOTSIM_PATH")); + + +// Load config file to a hashmap. These files are of the form (key: param, value: the distribution parameters). +std::map> loadConfig(const std::string& fileName) { + + std::map> configMap; // a 2D hashmap of the configuration + std::string fullPath = std::string(REPOPATH) + "/configs/" + fileName; // the full path to the configfile + + try { + + // Get file and parse like a map + YAML::Node config = YAML::LoadFile(fullPath); + if (config.IsMap()) { + + // Fill in 2D hashmap + for (const auto& element : config) { + std::string outerKey = element.first.as(); + const YAML::Node& innerNode = element.second; + + if (innerNode.IsMap()) { // If there is a 2nd dim for the specific element + std::map innerMap; + for (const auto& innerElement : innerNode) { + std::string innerKey = innerElement.first.as(); + std::string innerValue = innerElement.second.as(); + innerMap[innerKey] = innerValue; + } + configMap[outerKey] = innerMap; + } + } + } + } catch (const YAML::BadFile& e) { // fail to open + std::cerr << "Failed to open YAML file: " << fullPath << std::endl; + throw; + } catch (const YAML::Exception& e) { // issue parsing + std::cerr << "YAML parsing error: " << e.what() << std::endl; + throw; + } + + // Check if the map contains the key "dist" is assigned for each parameter + for (auto const& [property, propertyDist] : configMap) { + if (property != "global") { + if (propertyDist.find("dist") == propertyDist.end()) { + throw std::runtime_error("Key 'dist' not found in configuration for property " + property + "."); + } + } + } + + return configMap; +} // Get the largest number from a vector of filenames @@ -51,8 +106,7 @@ GravitationalEnvironment::GravitationalEnvironment(const std::vector lastLogFileNames; - const char* repoPath = std::getenv("HOOTSIM_PATH"); - std::string dataPath = std::string(repoPath) + "/data"; + std::string dataPath = REPOPATH + "/data"; for (const auto& entry : fs::directory_iterator(dataPath)) { if (fs::is_regular_file(entry.status())) { lastLogFileNames.push_back(entry.path().filename().string()); @@ -75,8 +129,7 @@ GravitationalEnvironment::GravitationalEnvironment(const std::vector lastLogFileNames; - const char* repoPath = std::getenv("HOOTSIM_PATH"); - std::string dataPath = std::string(repoPath) + "/data"; + std::string dataPath = REPOPATH + "/data"; for (const auto& entry : fs::directory_iterator(dataPath)) { if (fs::is_regular_file(entry.status())) { lastLogFileNames.push_back(entry.path().filename().string()); @@ -88,6 +141,114 @@ GravitationalEnvironment::GravitationalEnvironment(const std::vector lastLogFileNames; + std::string dataPath = REPOPATH + "/data"; + for (const auto& entry : fs::directory_iterator(dataPath)) { + if (fs::is_regular_file(entry.status())) { + lastLogFileNames.push_back(entry.path().filename().string()); + } + } + + // Get the largest log file number and create new log file + int lastLogNum = getLargestLabelNumber(lastLogFileNames, "run"); + logFileName = dataPath + "/run" + std::to_string(lastLogNum + 1) + ".csv"; + } + } +GravitationalEnvironment::GravitationalEnvironment(const std::string configFileName, const bool log, std::string logFilePrefix = "run") + : log(log), time(0) { + + // Get particles + loadParticlesFromConfig(configFileName); + + // Declare the number of particles + nParticles = particlePtrs.size(); + + // Create a log file if we want one + if (log == true) { + + // Get a vector of the filenames in the data directory + std::vector lastLogFileNames; + std::string dataPath = REPOPATH + "/data"; + for (const auto& entry : fs::directory_iterator(dataPath)) { + if (fs::is_regular_file(entry.status())) { + lastLogFileNames.push_back(entry.path().filename().string()); + } + } + + // Get the largest log file number and create new log file + int lastLogNum = getLargestLabelNumber(lastLogFileNames, logFilePrefix); + logFileName = dataPath + "/" + logFilePrefix + std::to_string(lastLogNum + 1) + ".csv"; + } + } + + +// Load a full environment from the configuration file +void GravitationalEnvironment::loadParticlesFromConfig(const std::string configFileName) { + + // Get configuration map + std::map> configMap = loadConfig(configFileName); + + // Grab the gloabl config params for the environment + std::map globalConfigMap = configMap.at("global"); + int nParticles = std::stoi(globalConfigMap.at("nParticles")); + + // Generate distributions for each param + std::map> envParams; + + // Iterate through the configuration and sample particles according to config prescirption + for (auto const& [property, propertyDist] : configMap) { + + if (property != "global") { // not the global configuarion field + + // The type of distribution for the property + std::string distType = propertyDist.at("dist"); + + // Sample from the distribution for the property + if (distType == "constant") { // Constant dist + envParams[property] = std::vector (nParticles, std::stod(propertyDist.at("val"))); + } else if (distType == "normal") { // Normal dist + double mu = std::stod(propertyDist.at("mu")); + double sigma = std::stod(propertyDist.at("sigma")); + std::normal_distribution<> normalDist(mu, sigma); + envParams[property] = sampleFromDistribution(nParticles, normalDist); + } else if (distType == "uniform") { // Uniform dist + double min = std::stod(propertyDist.at("min")); + double max = std::stod(propertyDist.at("max")); + std::uniform_real_distribution<> uniformDist(min, max); + envParams[property] = sampleFromDistribution(nParticles, uniformDist); + } else { + throw std::invalid_argument("Property " + property + " has an invalid distribution."); + } + } + } + + // Declare particle pointer vector + std::vector> positions(nParticles); + std::vector> velocities(nParticles); + for (int i = 0; i < nParticles; i++) { + + // Populate positions and velocities + positions[i] = {envParams.at("x")[i], envParams.at("y")[i], envParams.at("z")[i]}; + velocities[i] = {envParams.at("vx")[i], envParams.at("vy")[i], envParams.at("vz")[i]}; + double mass = envParams.at("mass")[i]; + + // Create Particle instance with pointers to elements in the positions and velocities vectors + particlePtrs.push_back(std::make_shared(&positions[i], &velocities[i], mass)); + } +} // Get the forces in the environment @@ -190,7 +351,7 @@ std::string GravitationalEnvironment::getStepLog() const { // Run a simulation void GravitationalEnvironment::simulate(const double duration, const double timestep) { - std::string logStr = getLogHeader(); + std::string logStr = getLogHeader() + ",\n"; std::cout << getLogHeader() + "\n"; // Get number of timesteps and take steps iteratively @@ -209,6 +370,7 @@ void GravitationalEnvironment::simulate(const double duration, const double time } // end state if (log == true) { + logStr += std::to_string(nTimesteps * timestep) + ","; logStr += getStepLog() + "\n"; } std::cout << nTimesteps * timestep << ",\t" << getStepLog() << "\n"; diff --git a/src/simulation.cpp b/src/simulation.cpp index aac8648..61715d5 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -6,42 +6,10 @@ #include "../include/environment.h" int main() { + GravitationalEnvironment defaultEnv("default.yaml", true); - // Parameters for simulation - float timestep = 0.5; - - // Initialize particle specs - double mass = 1E9; - std::array initial_position1 = {0, 0, 0}; - std::array initial_position2 = {4, 0, 0}; - std::array initial_velocity1 = {0, 0, 0}; - std::array initial_velocity2 = {0, 0, 0}; - - - // Define the particles - auto particle1Ptr = std::make_shared(&initial_position1, &initial_velocity1, mass); - auto particle2Ptr = std::make_shared(&initial_position2, &initial_velocity2, mass); - std::vector> particles = {particle1Ptr, particle2Ptr}; - - // Initialize an environment - GravitationalEnvironment env1(particles, true); - - // Take a step - env1.simulate(3, timestep); - - // std::cout << particle1.position[0] << "\n"; - // std::cout << particle2.position[0]; + // Simulate + defaultEnv.simulate(3, 0.5); return 0; - }; - - - - // for (int i=0; i<5; i++) { - // // Update the particle - // particle1.update(&force, timestep); - // } - - // // Define the force - // std::array force = {3, 0, 0}; diff --git a/src/statistics.cpp b/src/statistics.cpp new file mode 100644 index 0000000..4988d38 --- /dev/null +++ b/src/statistics.cpp @@ -0,0 +1,4 @@ +#include "../include/statistics.h" + +std::random_device rd; +std::mt19937 GENERATOR(rd()); diff --git a/test/test_environment.cpp b/test/test_environment.cpp index 51036df..727ec4d 100644 --- a/test/test_environment.cpp +++ b/test/test_environment.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -15,6 +16,7 @@ double _G = 6.6743e-11; // Data path const char* repoPath = std::getenv("HOOTSIM_PATH"); std::string dataPath = std::string(repoPath) + "/data"; +std::string configPath = std::string(repoPath) + "/config"; // Initialize particle specs double mass = 1E10; @@ -50,6 +52,11 @@ TEST_CASE("Environment Initialization") { CHECK(env2.logFileName == dataPath + "/funPrefix0.csv"); } +TEST_CASE("Error Handling for Load Config") { + // Test error handling when loading an invalid configuration file + CHECK_THROWS(loadConfig("invalid_config.yaml")); +} + void test_normalCase() { std::vector filenames = {"run0.csv", "run1.csv", "run10.csv"}; CHECK(getLargestLabelNumber(filenames, "run") == 10); @@ -161,5 +168,76 @@ TEST_CASE("Reset Environment") { // Reset environment env1.reset(); + // Check that it resets CHECK(env1.time == 0); } + + +////////// CONFIG FILE TESTS ////////// +TEST_CASE("Load Config File") { + + // Get the default map + std::map> defaultConfig = loadConfig("default.yaml"); + + // Check some values + CHECK(defaultConfig["vy"]["dist"] == "normal"); + CHECK(defaultConfig["vy"]["mu"] == "-2"); + CHECK(defaultConfig["vy"]["sigma"] == "4"); + CHECK(defaultConfig["global"]["nParticles"] == "1000"); +} + +void checkDefaultEnv(GravitationalEnvironment& env) { + // Check that it all looks good... that default vals are: +// mass: +// dist: constant +// val: 100000000 +// x: +// dist: uniform +// min: 0 +// max: 4 +// y: +// dist: uniform +// min: -2 +// max: 4 +// z: +// dist: uniform +// min: 1 +// max: 10 +// vx: +// dist: uniform +// min: 0 +// max: 4 +// vy: +// dist: uniform +// min: -2 +// max: 4 +// vz: +// dist: uniform +// min: -10 +// max: 10 + bool massCheck = true; + bool xCheck = true; + bool yCheck = true; + bool zCheck = true; + for (auto pPtr : env.particlePtrs) { + massCheck &= pPtr->mass == 100000000; + xCheck &= (pPtr->position[0] <= 4 & pPtr->position[0] >= 0); + yCheck &= (pPtr->position[1] <= 4 & pPtr->position[1] >= -2); + zCheck &= (pPtr->position[2] <= 10 & pPtr->position[2] >= 1); + } + CHECK(env.nParticles == 1000); + CHECK(xCheck); + CHECK(yCheck); + CHECK(zCheck); +} +TEST_CASE("Load Particles From Config") { + + // Get a defualt environment class + GravitationalEnvironment defaultEnv("default.yaml", true); + GravitationalEnvironment defaultEnv2("default.yaml", true, "prefixyprefix"); + + // Check 'em + checkDefaultEnv(defaultEnv); + checkDefaultEnv(defaultEnv2); + CHECK(defaultEnv2.logFileName == dataPath + "/prefixyprefix0.csv"); +} diff --git a/test/test_statistics.cpp b/test/test_statistics.cpp new file mode 100644 index 0000000..8c06c97 --- /dev/null +++ b/test/test_statistics.cpp @@ -0,0 +1,66 @@ +#include +#include +#include + +#include "../include/doctest.h" +#include "../include/statistics.h" + +TEST_CASE("Constant") { + // Sample + std::vector valVec(12, std::stod("6.32")); + + // Check that all the values are equal + bool allEq = true; + for (double v : valVec) { + allEq &= (v == 6.32); + } + CHECK(allEq); + CHECK(valVec.size() == 12); +} + +TEST_CASE("Normal") { + // Sample + int n = 1000; + std::normal_distribution<> normalDist(0, 1); + std::vector valVec = sampleFromDistribution(n, normalDist); + + // Check the mean is near 0 + double mean = 0; + for (double v : valVec) { + mean += v; + } + mean /= n; + CHECK(abs(mean) < 0.1); + + // Check std is near 1 + double std = 0; + for (double v : valVec) { + std += v*v; + } + std = pow(std / n, 0.5); + CHECK(abs(std - 1) < 0.1); + CHECK(valVec.size() == n); +} + +TEST_CASE("Uniform") { + // Sample + int n = 1000; + std::uniform_real_distribution<> uniformDis(-5, 5); + std::vector valVec = sampleFromDistribution(n, uniformDis); + + // Check the mean is near 0 + double mean = 0; + for (double v : valVec) { + mean += v; + } + mean /= n; + CHECK(abs(mean) < 0.25); + + // Check they're in the sampling range + bool inRange = true; + for (double v : valVec) { + inRange &= (v <= 5 & v >= -5); + } + CHECK(inRange); + CHECK(valVec.size() == n); +}