From 582f1f4d8b83fb21efd3cfcd9a1eb53ef5880cc2 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Sun, 16 Aug 2020 21:27:20 -0400 Subject: [PATCH] Added a new acceptor based on targets in a data file --- Makefile | 2 +- loadtxt.cpp | 4 ++- loadtxt.h | 2 +- scattering_experiment.cpp | 60 ++++++++++++++++++++++++++++++--------- 4 files changed, 52 insertions(+), 16 deletions(-) diff --git a/Makefile b/Makefile index 59ad505..37d98b0 100644 --- a/Makefile +++ b/Makefile @@ -5,7 +5,7 @@ LIB += -lgsl -pthread EXECUTABLE ?= scattering_experiment -defualt: loadtxt.o replay.o scattering_experiment.o +default: loadtxt.o replay.o scattering_experiment.o $(CXX) $(CPPFLAGS) $(CXXFLAGS) $(INC) $? -o $(EXECUTABLE) $(LIB) .cpp.o: diff --git a/loadtxt.cpp b/loadtxt.cpp index 6cc1ccb..f36d3bf 100644 --- a/loadtxt.cpp +++ b/loadtxt.cpp @@ -6,7 +6,7 @@ //TODO if cols is an empty vector, get all columns from the file //TODO error checking -Loadtxt::Loadtxt(std::string file_name, std::vector cols) +Loadtxt::Loadtxt(std::string file_name, std::vector cols, int skiprows) { std::sort(cols.begin(), cols.end()); n_cols = cols.size(); @@ -14,7 +14,9 @@ Loadtxt::Loadtxt(std::string file_name, std::vector cols) int n_rows_alloc = tmp_number_of_rows; buffer = (double*)malloc(n_cols * sizeof(double) * n_rows_alloc); std::ifstream file(file_name); + if (!file.good()) throw std::runtime_error("Input file not found (" + file_name + ")"); std::string line; + for (int row = 0; row < skiprows; row++) getline(file, line); int row = -1; while (getline(file, line)) { if (line[line.find_first_not_of(whitespaces)]=='#') continue; diff --git a/loadtxt.h b/loadtxt.h index 02ad1b3..54a1a59 100644 --- a/loadtxt.h +++ b/loadtxt.h @@ -3,7 +3,7 @@ #include class Loadtxt { public: - Loadtxt(std::string file_name, std::vector cols); + Loadtxt(std::string file_name, std::vector cols, int skiprows=0); ~Loadtxt(); std::vector> get_cols(); private: diff --git a/scattering_experiment.cpp b/scattering_experiment.cpp index 657d284..0d27d92 100644 --- a/scattering_experiment.cpp +++ b/scattering_experiment.cpp @@ -11,6 +11,8 @@ #include // We just want to print the timestamp #include +#include "loadtxt.h" + using double3 = Eigen::Vector3d; using Generator = std::default_random_engine; @@ -113,6 +115,42 @@ public: } }; +//////////////////////////////////////////////////////////////////////////////////////////////////////////// +class Match_data : public Acceptor { +public: + Match_data() : + data(Loadtxt("Peter_7D.txt.input", {2, 3, 4, 5, 6, 7}, 3).get_cols()), + n(data[0].size()), + tolerance(0.1) + { + constexpr double kms_to_kpcGyr = 1.0226911647958985; + for (unsigned i = 0; i < n; i++) { + data[3][i] *= kms_to_kpcGyr; + data[4][i] *= kms_to_kpcGyr; + data[5][i] *= kms_to_kpcGyr; + } + } + bool operator()(const double w[6]) const override + { + for (unsigned i = 0; i < n; i++) { + for (unsigned k = 0; k < 6; k++) { + if (abs((w[k] - data[k][i])/w[k]) > 0.1) break; + return true; + } + } + return false; + } + void output_information(std::ostream& stream) const override + { + stream << "# Acceptor type: match data\n"; + stream << "# tolerance = " << std::setw(13) << tolerance << '\n'; + } +private: + std::vector> data; + unsigned n; + double tolerance; +}; + class Initial_conditions { public: virtual std::array operator()(Generator& generator) const = 0; @@ -171,14 +209,14 @@ public: Box_initial_conditions(double x_max, double v_max) : x_max(x_max), v_max(v_max) { - position_distribution = std::uniform_real_distribution(-x_max, x_max); - velocity_distribution = std::uniform_real_distribution(-v_max, v_max); + position_distribution = std::uniform_real_distribution(-x_max, x_max); + velocity_distribution = std::uniform_real_distribution(-v_max, v_max); } std::array operator()(Generator& generator) const override { 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); + for (unsigned i=0; i<3; i++) result[i] = position_distribution(generator); + for (unsigned i=3; i<6; i++) result[i] = velocity_distribution(generator); return result; } void output_information(std::ostream& stream) const override @@ -191,7 +229,7 @@ public: } private: double x_max, v_max; - mutable std::uniform_real_distribution<> position_distribution, velocity_distribution; + mutable std::uniform_real_distribution position_distribution, velocity_distribution; }; class Scattering_experiment { @@ -298,17 +336,13 @@ int main() std::cout << "Hi\n"; const double t_max = 11.65551390; // Gyr - 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 - In_box_cyl accept(5, 50, target_pos, target_vel); Uniform_sphere_escape_velocity_cutoff uniform_sphere_escape_velocity_cutoff(75, 1.82777387E+07, 9.86242602E+00); - Scattering_experiment scattering_experiment("results.dat", t_max, 1024*256, uniform_sphere_escape_velocity_cutoff, Accept_all()); - //scattering_experiment.set_concurrency(1); + auto zzzzz = Match_data(); + + Scattering_experiment scattering_experiment("results.dat", t_max, 1024*1024, uniform_sphere_escape_velocity_cutoff, std::ref(zzzzz)); + scattering_experiment.set_concurrency(80); scattering_experiment.write_header(); scattering_experiment.launch();