Skip to content

Commit

Permalink
Merge pull request #700 from beomki-yeo/convergence-test
Browse files Browse the repository at this point in the history 
Print ridders results and make all tracks converged
  • Loading branch information
@beomki-yeo
beomki-yeo authored Mar 16, 2024
2 parents ec3c836 + d68f494 commit 86a36f4
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 32 deletions.
35 changes: 23 additions & 12 deletions tests/integration_tests/cpu/propagator/jacobian_validation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ const std::array<scalar, 5u> h_sizes_rect{3e0f, 3e0f, 2e-2f, 1e-3f, 1e-3f};
const std::array<scalar, 5u> h_sizes_wire{3e0f, 3e0f, 2e-2f, 1e-3f, 1e-3f};

// Ridders' algorithm setup
constexpr const unsigned int Nt = 50u;
constexpr const unsigned int Nt = 100u;
const std::array<scalar, 5u> safe{5.0f, 5.0f, 5.0f, 5.0f, 5.0f};
const std::array<scalar, 5u> con{1.2f, 1.2f, 1.2f, 1.2f, 1.2f};
constexpr const scalar big = std::numeric_limits<scalar>::max();
Expand Down Expand Up @@ -94,6 +94,18 @@ std::uniform_int_distribution<int> rand_bool(0, 1);
std::uniform_real_distribution<scalar> rand_gamma(0.f,
2.f * constant<scalar>::pi);

// Correlation factor in the range of [-10%, 10%]
constexpr scalar min_corr = -0.1f;
constexpr scalar max_corr = 0.1f;

// Values for sampling standard deviation
const std::array<scalar, 6u> stddevs_sampling{50.f * unit<scalar>::um,
50.f * unit<scalar>::um,
1.f * unit<scalar>::mrad,
1.f * unit<scalar>::mrad,
0.01f,
1.f * unit<scalar>::ns};

// surface types
using rect_type = rectangle2D;
using wire_type = line_square;
Expand Down Expand Up @@ -244,23 +256,22 @@ bound_covariance_type get_random_initial_covariance(const scalar ini_qop) {
bound_covariance_type ini_cov =
matrix_operator().template zero<e_bound_size, e_bound_size>();

// Correlation factor in the range of [-10%, 10%]
scalar min_corr = -0.1f;
scalar max_corr = 0.1f;
// Random correction factor
std::uniform_real_distribution<scalar> rand_corr(min_corr, max_corr);

// Distribution for sampling standard deviations
std::normal_distribution<scalar> rand_l0(0.f * unit<scalar>::um,
50.f * unit<scalar>::um);
stddevs_sampling[0u]);
std::normal_distribution<scalar> rand_l1(0.f * unit<scalar>::um,
50.f * unit<scalar>::um);
stddevs_sampling[1u]);
std::normal_distribution<scalar> rand_phi(0.f * unit<scalar>::mrad,
1.f * unit<scalar>::mrad);
stddevs_sampling[2u]);
std::normal_distribution<scalar> rand_theta(0.f * unit<scalar>::mrad,
1.f * unit<scalar>::mrad);
std::normal_distribution<scalar> rand_qop(0.f, 0.01f * ini_qop);
stddevs_sampling[3u]);
std::normal_distribution<scalar> rand_qop(0.f,
stddevs_sampling[4u] * ini_qop);
std::normal_distribution<scalar> rand_time(0.f * unit<scalar>::ns,
1.f * unit<scalar>::ns);
stddevs_sampling[5u]);

/*
// Typical stddev range taken from the figures of ATL-PHYS-PUB-2021-024 and
Expand Down Expand Up @@ -1902,9 +1913,9 @@ int main(int argc, char** argv) {
EXPECT_EQ(dqopdqop_rel_diffs_rect[i].size(), n_tracks);
EXPECT_EQ(dqopdqop_rel_diffs_wire[i].size(), n_tracks);

EXPECT_GE(statistics::mean(dqopdqop_rel_diffs_rect[i]), 1e-10f);
EXPECT_GE(statistics::mean(dqopdqop_rel_diffs_rect[i]), 1e-12f);
EXPECT_LE(statistics::mean(dqopdqop_rel_diffs_rect[i]), 1e-2f);
EXPECT_GE(statistics::mean(dqopdqop_rel_diffs_wire[i]), 1e-10f);
EXPECT_GE(statistics::mean(dqopdqop_rel_diffs_wire[i]), 1e-12f);
EXPECT_LE(statistics::mean(dqopdqop_rel_diffs_wire[i]), 1e-2f);
}
}
Expand Down
14 changes: 11 additions & 3 deletions tests/validation/root/jacobian_comparison.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ double y_min = -15;
double y_max = 10;
double y_margin = 0;
double header_size = 0.05;
std::array<float, 4> ldim{0.545075, 0.621849, 0.942404, 0.881048};
std::array<float, 4> ldim{0.54424, 0.621849, 0.942404, 0.880252};
double pad_x0 = 0.00;
double pad_x1 = 1;
double pad_y0 = 0.00;
Expand Down Expand Up @@ -171,6 +171,14 @@ TH1D* get_histogram(std::string name, const int n_labels,

auto rdf = ROOT::RDF::FromCSV(csv_name);
auto rdf_means = get_means(rdf);

// Count the number of non-convergence event
auto conv_success = rdf.Filter("total_convergence == 1").Count();
auto conv_failure = rdf.Filter("total_convergence == 0").Count();

std::cout << "Convergence events: " << *conv_success << std::endl;
std::cout << "Non-convergence events: " << *conv_failure << std::endl;

TH1D* histo = new TH1D(histo_name.c_str(), histo_name.c_str(), 3, 0, 3);
histo->GetYaxis()->SetRangeUser(y_min - y_margin, y_max + y_margin);
histo->GetYaxis()->SetLabelSize(0);
Expand Down Expand Up @@ -225,7 +233,7 @@ void draw_pad(const std::string& pad_name) {
void draw_text(const std::string& text) {

const float x1 = 1.23427;
const float y1 = 6.511;
const float y1 = 6.44408;

TLatex* ttext = new TLatex(0.f, 0.f, text.c_str());
ttext->SetTextFont(132);
Expand Down Expand Up @@ -266,7 +274,7 @@ void jacobian_comparison() {

gStyle->SetOptTitle(0);
const std::array<float, 2> cdim{1200, 500};
const std::array<int, 4> markers{kOpenCross, kOpenTriangleUp, kOpenSquare,
const std::array<int, 4> markers{kOpenCross, kOpenSquare, kOpenTriangleUp,
kOpenCircle};
const std::array<int, 4> hues{kOrange + 8, kMagenta + 1, kAzure,
kGreen + 2};
Expand Down
72 changes: 55 additions & 17 deletions tests/validation/root/rk_tolerance_comparison.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,6 @@ double title_font_size_rk_tol = 0.046;
int title_font = 132;
int label_font = 132;
int legend_font = 132;
double marker_size = 1.9;
double pad_x0 = 0.005f;
double pad_x1 = 1.f;
double pad_y0 = 0.005f;
Expand Down Expand Up @@ -189,7 +188,9 @@ void draw_graphs(const std::string header_title, const std::string geom_title,
TGraph* gr[25];
TMultiGraph* mg = new TMultiGraph();

const std::array<int, 5u> marker_styles = {7, 2, 5, 27, 32};
const std::array<int, 5u> marker_styles = {7, 2, 5, 27, 25};
const std::array<double, 5u> marker_sizes = {2.135, 2.135, 2.135, 2.135,
1.49};
const std::array<int, 5u> line_styles = {1, 3, 2, 7, 4};
const std::array<int, 5u> hues = {kOrange + 2, kPink + 5, kBlue + 2,
kCyan + 2, kGreen + 2};
Expand All @@ -211,7 +212,7 @@ void draw_graphs(const std::string header_title, const std::string geom_title,
const int m = i % 5;

gr[i]->SetMarkerStyle(marker_styles[n]);
gr[i]->SetMarkerSize(marker_size);
gr[i]->SetMarkerSize(marker_sizes[n]);
gr[i]->SetLineStyle(line_styles[m]);
gr[i]->SetMarkerColor(hues[m]);
gr[i]->SetLineColor(hues[m]);
Expand Down Expand Up @@ -247,6 +248,10 @@ void draw_graphs(const std::string header_title, const std::string geom_title,
double x_max = x_vec.back();

if (x_vec.size() > 10) {
if (x_vec.size() == 13) {
x_margin = 2;
}

x_min = x_min - x_margin;
x_max = x_max + x_margin;
mg->GetXaxis()->SetLimits(x_min, x_max);
Expand All @@ -255,13 +260,23 @@ void draw_graphs(const std::string header_title, const std::string geom_title,

mg->Draw("APL");

auto ga = new TGaxis(x_vec.front(), yaxis_min - yaxis_margin,
x_vec.back(), yaxis_min - yaxis_margin,
x_vec.front(), x_vec.back(), 405, "N");
ga->SetLabelFont(label_font);
ga->SetLabelSize(label_font_size_rk_tol);
ga->SetLabelOffset(-0.0065);
ga->Draw();
if (x_vec.size() == 11) {
auto ga = new TGaxis(x_vec.front(), yaxis_min - yaxis_margin,
x_vec.back(), yaxis_min - yaxis_margin,
x_vec.front(), x_vec.back(), 405, "N");
ga->SetLabelFont(label_font);
ga->SetLabelSize(label_font_size_rk_tol);
ga->SetLabelOffset(-0.0065);
ga->Draw();
} else if (x_vec.size() == 13) {
auto ga = new TGaxis(x_vec.front(), yaxis_min - yaxis_margin,
x_vec.back(), yaxis_min - yaxis_margin,
x_vec.front(), x_vec.back(), 304, "N");
ga->SetLabelFont(label_font);
ga->SetLabelSize(label_font_size_rk_tol);
ga->SetLabelOffset(-0.0065);
ga->Draw();
}

mg->GetYaxis()->SetLabelSize(0);
mg->GetYaxis()->SetTickLength(0);
Expand All @@ -284,6 +299,10 @@ void draw_graphs(const std::string header_title, const std::string geom_title,
double rk_title_deltaX =
(x_vec.back() - x_vec.front()) * rk_title_offset_fraction;

if (x_vec.size() == 13) {
rk_title_deltaX = -0.2;
}

legend->Draw();
draw_text(rk_title_deltaX + x_vec.front(), rk_title_y, rk_ygap,
rk_header_text_size, rk_geom_text_size, header_title, geom_title);
Expand All @@ -309,6 +328,11 @@ void draw_mean_step_size(const std::string header_title,
gr->GetXaxis()->SetTitle("log_{10}(#font[12]{#tau} [mm])");
gr->GetYaxis()->SetTitle("log_{10}(Mean of avg. step size [mm])");
gr->GetXaxis()->SetLimits(x_vec.front() - 0.5, x_vec.back() + 0.5);

if (x_vec.size() == 13) {
ymin = -4;
}

gr->GetYaxis()->SetRangeUser(ymin, ymax);
gr->GetYaxis()->SetNdivisions(505);
gr->GetXaxis()->SetLabelSize(label_font_size_step);
Expand All @@ -327,18 +351,32 @@ void draw_mean_step_size(const std::string header_title,
gr->GetYaxis()->SetLabelFont(label_font);

if (x_vec.size() > 10) {
gr->GetXaxis()->SetLimits(x_vec.front() - 1, x_vec.back() + 1);
if (x_vec.size() == 13) {
x_margin = 2;
}

gr->GetXaxis()->SetLimits(x_vec.front() - x_margin,
x_vec.back() + x_margin);
gr->GetXaxis()->SetLabelSize(0);
gr->GetXaxis()->SetTickLength(0);

gr->Draw("APL");

auto ga = new TGaxis(x_vec.front(), ymin, x_vec.back(), ymin,
x_vec.front(), x_vec.back(), 405, "N");
ga->SetLabelFont(label_font);
ga->SetLabelSize(label_font_size_step);
ga->SetLabelOffset(x_label_offset);
ga->Draw();
if (x_vec.size() == 11) {
auto ga = new TGaxis(x_vec.front(), ymin, x_vec.back(), ymin,
x_vec.front(), x_vec.back(), 405, "N");
ga->SetLabelFont(label_font);
ga->SetLabelSize(label_font_size_step);
ga->SetLabelOffset(x_label_offset);
ga->Draw();
} else if (x_vec.size() == 13) {
auto ga = new TGaxis(x_vec.front(), ymin, x_vec.back(), ymin,
x_vec.front(), x_vec.back(), 304, "N");
ga->SetLabelFont(label_font);
ga->SetLabelSize(label_font_size_step);
ga->SetLabelOffset(x_label_offset);
ga->Draw();
}

} else {
gr->Draw();
Expand Down

0 comments on commit 86a36f4

Please sign in to comment.