From b67c2cf2cb162a9d6d88d44152f1aa87315b78c5 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Wed, 6 May 2020 11:58:53 -0400 Subject: [PATCH] Added scattering experiment --- .gitignore | 3 +- Makefile | 10 +-- main.cpp => replay.cpp | 60 ++++++++++++-- scattering_experiment.cpp | 160 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 222 insertions(+), 11 deletions(-) rename main.cpp => replay.cpp (74%) create mode 100644 scattering_experiment.cpp diff --git a/.gitignore b/.gitignore index 60b2039..d0c6438 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ main libmain.so data __pycache__ -*.png \ No newline at end of file +*.png +*.o \ No newline at end of file diff --git a/Makefile b/Makefile index 02addc9..94d4882 100644 --- a/Makefile +++ b/Makefile @@ -1,15 +1,15 @@ OPTIMIZATION ?= 3 -CXXFLAGS += -O$(OPTIMIZATION) +CXXFLAGS += -O$(OPTIMIZATION) -pthread INC += -I/home/meiron/boost_1_72_0 LIB += -lgsl EXECUTABLE ?= main -default: - $(CXX) $(CPPFLAGS) $(CXXFLAGS) $(INC) loadtxt.cpp main.cpp -o $(EXECUTABLE) $(LIB) +defualt: loadtxt.o replay.o scattering_experiment.o + $(CXX) $(CPPFLAGS) $(CXXFLAGS) $(INC) $? -o $(EXECUTABLE) $(LIB) -lib: - $(CXX) $(CPPFLAGS) $(CXXFLAGS) -fPIC $(INC) loadtxt.cpp main.cpp -shared -o lib$(EXECUTABLE).so $(LIB) +.cpp.o: + $(CXX) $(CPPFLAGS) $(CXXFLAGS) $(INC) -c $< clean: rm -f *.o *.so $(EXECUTABLE) diff --git a/main.cpp b/replay.cpp similarity index 74% rename from main.cpp rename to replay.cpp index 16e84b3..5d4e40b 100644 --- a/main.cpp +++ b/replay.cpp @@ -144,14 +144,64 @@ void integrate(const double y0[], const double t_max, const double stride_size, } } -int main() + +#include +#include +#include +int mainXXX() { std::cout << "Hi" << std::endl; - double y[12]; - double y0[] = {80,0,0,0,80,0}; - integrate(y0, 10, 10, y); - for (int i=0; i<12; i++) std::cout << y[i] << std::endl; + double x_max = 100; + double v_max = 200; + double t_max = 10.; + + double3 target_pos = {-9.1817914423447837E+03, -2.3161972143758221E+03, 4.2321813890923086E+03}; // pc + double3 target_vel = { 1.6995805137819784E+02, 1.3476658021349482E+02, -1.3342192079702110E+02}; // km/s + + // Unit adjustment + target_pos *= 0.001; // Now in kpc + target_vel *= 1.0226911647958985; // Now in kpc/Gyr + + double tolrance = 2; + double max_dx = target_pos.norm()*tolrance; + double max_dv = target_vel.norm()*tolrance; + + const auto current_time = std::chrono::system_clock::now().time_since_epoch().count(); + const unsigned long big_prime = 840580612873; + const auto pid = getpid(); + const auto seed = current_time + big_prime*pid; // TODO add thread number + + std::default_random_engine generator(seed); + + std::uniform_real_distribution position_distribution(-x_max, x_max); + auto random_position_component = [&]() { return position_distribution(generator); }; + + std::uniform_real_distribution velocity_distribution(-v_max, v_max); + auto random_velocity_component = [&]() { return velocity_distribution(generator); }; + + for (int i=0; i<100000; i++) { + if (i%1000==0) std::cout << i << "\n"; + double y[12], y0[6]; + for (int i=0; i<3; i++) y0[i] = random_position_component(); + for (int i=3; i<6; i++) y0[i] = random_velocity_component(); + + integrate(y0, t_max, t_max, y); + double3 new_pos(y+6); + double3 new_vel(y+9); + + + if (((new_pos - target_pos).norm() < max_dx) && ((new_vel - target_vel).norm() < max_dv)) + std::cout << "*"; + } + + + + + // double y[12]; + // double y0[] = {80,0,0,0,80,0}; + // integrate(y0, 10, 10, y); + // for (int i=0; i<12; i++) std::cout << y[i] << std::endl; diff --git a/scattering_experiment.cpp b/scattering_experiment.cpp new file mode 100644 index 0000000..7d43715 --- /dev/null +++ b/scattering_experiment.cpp @@ -0,0 +1,160 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include // We just want to print the timestamp +#include + +extern "C" +void integrate(const double y0[], const double t_max, const double stride_size, double y[]); + +std::string git_version() { + std::string result; + std::ifstream git_log(".git/logs/HEAD"); + if (git_log) { + std::string prev_line, line; + getline(git_log, prev_line); + while (getline(git_log, line)) { + prev_line = line; + } + if (prev_line.length() > 51) result = prev_line.substr(41, 10); + else result = "UNKNOWN"; + } else result = "UNKNOWN"; + return result; +} + +std::string iso_time(std::chrono::system_clock::time_point time) +{ + auto time_c = std::chrono::system_clock::to_time_t(time); + auto timeinfo = gmtime(&time_c); + char time_string[256]; + strftime(time_string,256,"%FT%TZ", timeinfo); + return std::string(time_string); +} + +class Initial_conditions { +public: + Initial_conditions(unsigned long long seed, double x_max, double v_max) + : x_max(x_max), v_max(v_max) + { + generator = std::default_random_engine(seed); + position_distribution = std::uniform_real_distribution(-x_max, x_max); + velocity_distribution = std::uniform_real_distribution(-v_max, v_max); + } + std::array generate() + { + std::array result; + for (int i=0; i<3; i++) result[i] = position_distribution(generator); + for (int i=3; i<6; i++) result[i] = velocity_distribution(generator); + return result; + } +private: + double x_max, v_max; + std::default_random_engine generator; + std::uniform_real_distribution position_distribution, velocity_distribution; +}; + +class Scattering_experiment { +public: + Scattering_experiment(std::string file_name) + : bunch_size(16384), counter_target(10000000L), report_interval_seconds(3) + { + auto time = std::chrono::system_clock::now(); + file = std::ofstream(file_name/*, std::ofstream::app*/); + file << "#############################################\n"; + file << "# Cosmological replay scattering experiment #\n"; + file << "#############################################\n"; + file << "# Time: " << iso_time(time) << '\n'; + file << "# Last Git commit: " << git_version() << '\n'; + char hostname[256]; + gethostname(hostname, 256); + file << "# Hostname: " << hostname << '\n'; + const auto pid = getpid(); + file << "# PID: " << pid << '\n'; + const unsigned long time_count = time.time_since_epoch().count(); + const unsigned long big_prime = 840580612873L; + unique_id = time_count + big_prime*pid; // TODO hash it + file << "# Unique identifier: " << unique_id << '\n'; + n_threads = 4; + file << "# Concurrency: " << n_threads << '\n'; + // TODO print information about target and galaxy model + file.flush(); + } + ~Scattering_experiment() { + file.close(); + } + void output_data(int aaa) + { + std::lock_guard lock(mtx); + file << aaa << '\n'; + file.flush(); + } + void thread_task(int tid) + { + auto thread_begin_time = std::chrono::system_clock::now(); + auto report_interval = std::chrono::seconds(report_interval_seconds); + auto next_report_time = thread_begin_time + report_interval; + + // std::default_random_engine generator(unique_id + tid); + // std::uniform_real_distribution position_distribution(0, 1000); + // auto random_position_component = [&]() { return position_distribution(generator); }; + Initial_conditions initial_conditions(unique_id + tid, 100, 100); + + while (counter < counter_target) { + for (int i=0; i= next_report_time) { + std::lock_guard lock(mtx); + file << "# " << iso_time(now) << " counter = " << counter << std::endl; + next_report_time = now + report_interval; + } + } + } + } + void launch() + { + counter = 0; + for (int tid=0; tid threads; + std::atomic counter; + unsigned long long counter_target; + unsigned long unique_id; + int bunch_size; + int report_interval_seconds; + + int i; +}; + + + +int main() +{ + std::cout << "Hi\n"; + + Scattering_experiment scattering_experiment("aaa.out"); + scattering_experiment.launch(); + + std::cout << "Bye\n"; + + return 0; +} \ No newline at end of file