Merge pull request #702 from beomki-yeo/record-convergence-iteration

Record the number of iterations from Ridders algorithm

Author: beomki-yeo

tests/integration_tests/cpu/propagator/jacobian_validation.cpp

@@ -571,7 +571,9 @@ bound_track_parameters<transform3_type>::covariance_type directly_differentiate(
     const scalar detector_length, const field_t& field,
     const scalar overstep_tolerance, const scalar on_surface_tolerance,
     const scalar rk_tolerance, const scalar constraint_step,
-    const std::array<scalar, 5u> hs, std::array<bool, 25>& convergence) {
+    const std::array<scalar, 5u> hs,
+    std::array<unsigned int, 5u>& num_iterations,
+    std::array<bool, 25>& convergence) {
 
     // Return Jacobian
     bound_covariance_type differentiated_jacobian;
@@ -610,6 +612,7 @@ bound_track_parameters<transform3_type>::covariance_type directly_differentiate(
             ridder.run(nvec1, nvec2, delta, p, i, differentiated_jacobian);
 
             if (ridder.is_complete()) {
+                num_iterations[i] = p;
                 break;
             }
         }
@@ -684,6 +687,7 @@ void evaluate_jacobian_difference(
     std::ofstream& file, scalar& ref_rel_diff, bool use_field_gradient,
     bool do_inspect, const bool use_precal_values = false,
     [[maybe_unused]] bound_covariance_type precal_diff_jacobi = {},
+    [[maybe_unused]] std::array<unsigned int, 5u> precal_num_iterations = {},
     [[maybe_unused]] std::array<bool, 25u> precal_convergence = {}) {
 
     const auto phi0 = track.phi();
@@ -735,22 +739,29 @@ void evaluate_jacobian_difference(
          << "," << final_param.phi() << "," << final_param.theta() << ","
          << final_param.qop() << ",";
 
-    std::array<bool, 25u> convergence;
     bound_covariance_type differentiated_jacobian;
+    std::array<unsigned int, 5u> num_iterations;
+    std::array<bool, 25u> convergence;
 
     if (use_precal_values) {
-        convergence = precal_convergence;
         differentiated_jacobian = precal_diff_jacobi;
+        num_iterations = precal_num_iterations;
+        convergence = precal_convergence;
     } else {
         differentiated_jacobian = directly_differentiate<propagator_t, field_t>(
             trk_count, reference_param, det, detector_length, field,
             overstep_tolerance, on_surface_tolerance, rk_tolerance_dis,
-            constraint_step, hs, convergence);
+            constraint_step, hs, num_iterations, convergence);
     }
 
     bool total_convergence =
         (std::count(convergence.begin(), convergence.end(), false) == 0);
 
+    // Ridders number of iterations
+    for (unsigned int i = 0; i < 5u; i++) {
+        file << num_iterations[i] << ",";
+    }
+
     // Convergence
     file << total_convergence << ",";
     for (unsigned int i = 0; i < 25u; i++) {
@@ -1098,7 +1109,7 @@ void evaluate_jacobian_difference_helix(
      * ****************************/
 
     bound_covariance_type differentiated_jacobian;
-
+    std::array<unsigned int, 5u> num_iterations;
     std::array<bool, 25u> convergence;
 
     for (unsigned int i = 0u; i < 5u; i++) {
@@ -1128,6 +1139,7 @@ void evaluate_jacobian_difference_helix(
             ridder.run(nvec1, nvec2, delta, p, i, differentiated_jacobian);
 
             if (ridder.is_complete()) {
+                num_iterations[i] = p;
                 break;
             }
         }
@@ -1151,6 +1163,11 @@ void evaluate_jacobian_difference_helix(
          << getter::element(bound_vec, e_bound_theta, 0u) << ","
          << getter::element(bound_vec, e_bound_qoverp, 0u) << ",";
 
+    // Ridders number of iterations
+    for (unsigned int i = 0; i < 5u; i++) {
+        file << num_iterations[i] << ",";
+    }
+
     // Convergence
     file << total_convergence << ",";
     for (unsigned int i = 0; i < 25u; i++) {
@@ -1224,6 +1241,13 @@ void setup_csv_header_jacobian(std::ofstream& file) {
     // Final Parameter at the destination surface
     file << "l0_F,l1_F,phi_F,theta_F,qop_F,";
 
+    // Number of iterations to complete the numerical differentiation
+    file << "num_iterations_l0,"
+         << "num_iterations_l1,"
+         << "num_iterations_phi,"
+         << "num_iterations_theta,"
+         << "num_iterations_qop,";
+
     // Convergence
     file << "total_convergence,";
     file << "dl0dl0_C,dl0dl1_C,dl0dphi_C,dl0dtheta_C,dl0dqop_C,";
@@ -1767,13 +1791,14 @@ int main(int argc, char** argv) {
 
             if (rk_tolerance_iterate_mode) {
 
+                std::array<unsigned int, 5u> num_iterations;
                 std::array<bool, 25u> convergence;
                 auto differentiated_jacobian =
                     directly_differentiate<inhom_field_rect_propagator_t>(
                         track_count, rect_bparam, rect_det_w_mat,
                         detector_length, inhom_bfield, overstep_tol,
                         on_surface_tol, rk_tol_dis, constraint_step_size,
-                        h_sizes_rect, convergence);
+                        h_sizes_rect, num_iterations, convergence);
 
                 for (std::size_t i = 0u; i < log10_tols.size(); i++) {
 
@@ -1784,7 +1809,7 @@ int main(int argc, char** argv) {
                         std::pow(10.f, log10_tols[i]), rk_tol_dis,
                         constraint_step_size, h_sizes_rect, rect_files[i],
                         ref_rel_diff, true, do_inspect, true,
-                        differentiated_jacobian, convergence);
+                        differentiated_jacobian, num_iterations, convergence);
 
                     dqopdqop_rel_diffs_rect[i].push_back(ref_rel_diff);
                 }
@@ -1848,13 +1873,14 @@ int main(int argc, char** argv) {
 
             if (rk_tolerance_iterate_mode) {
 
+                std::array<unsigned int, 5u> num_iterations;
                 std::array<bool, 25u> convergence;
                 auto differentiated_jacobian =
                     directly_differentiate<inhom_field_wire_propagator_t>(
                         track_count, wire_bparam, wire_det_w_mat,
                         detector_length, inhom_bfield, overstep_tol,
                         on_surface_tol, rk_tol_dis, constraint_step_size,
-                        h_sizes_wire, convergence);
+                        h_sizes_wire, num_iterations, convergence);
 
                 for (std::size_t i = 0u; i < log10_tols.size(); i++) {
                     // Wire Inhomogeneous field with Material
@@ -1864,7 +1890,7 @@ int main(int argc, char** argv) {
                         std::pow(10.f, log10_tols[i]), rk_tol_dis,
                         constraint_step_size, h_sizes_wire, wire_files[i],
                         ref_rel_diff, true, do_inspect, true,
-                        differentiated_jacobian, convergence);
+                        differentiated_jacobian, num_iterations, convergence);
 
                     dqopdqop_rel_diffs_wire[i].push_back(ref_rel_diff);
                 }