Added scattering experiment

.gitignore

@@ -3,4 +3,5 @@ main
 libmain.so
 data
 __pycache__
-*.png
+*.png
+*.o

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)

replay.cpp

@@ -144,14 +144,64 @@ void integrate(const double y0[], const double t_max, const double stride_size,
     }
 }
 
-int main()
+
+#include <random>
+#include <chrono>
+#include <unistd.h>
+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<double> position_distribution(-x_max, x_max);
+    auto random_position_component = [&]() { return position_distribution(generator); };
+
+    std::uniform_real_distribution<double> 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;
 
 
 

scattering_experiment.cpp

@@ -0,0 +1,160 @@
+#include <array>
+#include <atomic>
+#include <fstream>
+#include <iostream>
+#include <mutex>
+#include <random>
+#include <string>
+#include <thread>
+#include <time.h> // We just want to print the timestamp
+#include <unistd.h>
+
+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<double>(-x_max, x_max);
+        velocity_distribution = std::uniform_real_distribution<double>(-v_max, v_max);
+    }
+    std::array<double,6> generate()
+    {
+        std::array<double,6> 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<double> 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<std::mutex> 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<double> 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<bunch_size; i++) {
+                auto ic = initial_conditions.generate();
+                /* INTEGRATE */
+                /* CHECK IF RESULT IS A MATCH */
+            }
+            counter += bunch_size;
+            if (tid == 0) {
+                auto now = std::chrono::system_clock::now();
+                if (now >= next_report_time) {
+                    std::lock_guard<std::mutex> 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<n_threads; tid++) {
+            threads.push_back(std::thread(&Scattering_experiment::thread_task, this, tid));
+        }
+        for (auto& t : threads) {
+            t.join();
+        }
+        std::cout << counter << '\n';
+    }
+    std::ofstream file;
+    std::mutex mtx;
+    int n_threads;
+    std::vector<std::thread> threads;
+    std::atomic<unsigned long long> 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;
+}