@@ -310,6 +310,8 @@ inline auto nonosc_evolve(SolverInfo &&info, Scalar xi, Scalar xf,
310310 * equations.
311311 * 6. @ref Eigen::Matrix<std::complex<Scalar>, -1, 1>: A vector containing the
312312 * interpolated solution at the specified `x_eval` The function returns these
313+ * 7. An array of logging information for the solver process. This can be added with the
314+ * info() from the `SolverInfo` object to get the total log information.
313315 * vectors encapsulated in a standard tuple, providing comprehensive information
314316 * about the solution process, including where the solution was evaluated, the
315317 * values and derivatives of the solution at those points, success status of
@@ -319,15 +321,16 @@ template <typename SolverInfo, typename Scalar, typename Vec>
319321inline auto evolve (SolverInfo &info, Scalar xi, Scalar xf,
320322 std::complex <Scalar> yi, std::complex <Scalar> dyi,
321323 Scalar eps, Scalar epsilon_h, Scalar init_stepsize,
322- Vec &&x_eval, bool hard_stop = false ) {
324+ Vec &&x_eval, bool hard_stop = false ,
325+ LogLevel log_level = riccati::LogLevel::ERROR) {
323326 static_assert (std::is_floating_point<Scalar>::value,
324327 " Scalar type must be a floating-point type." );
325328 using vectord_t = vector_t <Scalar>;
326329 Scalar direction = init_stepsize > 0 ? 1 : -1 ;
327330 if (init_stepsize * (xf - xi) < 0 ) {
328331 throw std::domain_error (
329332 " Direction of integration does not match stepsize sign,"
330- " adjusting it so that integration happens from xi to xf." );
333+ " adjust the direction or stepsize so that integration happens from xi to xf." );
331334 }
332335 // Check that yeval and x_eval are right size
333336 constexpr bool dense_output = compile_size_v<Vec> != 0 ;
@@ -436,9 +439,9 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
436439 matrixc_t y_eval;
437440 matrixc_t dy_eval;
438441 std::pair<complex_t , complex_t > a_pair;
442+ std::array<std::pair<LogInfo, std::size_t >, 5 > solver_counts = info.info ();
439443 while (std::abs (xcurrent - xf) > Scalar (1e-8 )
440444 && direction * xcurrent < direction * xf) {
441- // std::cout << "t = " << xcurrent << std::endl;
442445 Scalar phase{0.0 };
443446 bool success = false ;
444447 bool steptype = true ;
@@ -447,7 +450,6 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
447450 arena_matrix<vectorc_t > un (info.alloc_ , omega_n.size (), 1 );
448451 arena_matrix<vectorc_t > d_un (info.alloc_ , omega_n.size (), 1 );
449452 if (direction * hosc > direction * hslo){
450- // && (direction * hosc * wnext / (2.0 * pi<Scalar>()) > 1.0)) {
451453 if (hard_stop) {
452454 if (direction * (xcurrent + hosc) > direction * xf) {
453455 hosc = xf - xcurrent;
@@ -463,13 +465,27 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
463465 std::tie (success, y, dy, err, phase, un, d_un, a_pair)
464466 = osc_step<dense_output>(info, omega_n, gamma_n, xcurrent, hosc,
465467 yprev, dyprev, eps);
466- // std::cout << "Attempted osc step with hosc = " << hosc << ", successful? " << success << std::endl;
467468 steptype = 1 ;
469+ solver_counts[get_idx (LogInfo::RICCSTEP)].second ++;
470+ }
471+ if (unlikely (log_level == LogLevel::INFO && !success)) {
472+ info.logger ().template log <LogLevel::INFO>(
473+ std::string (" [Non-oscillatory step][x = " ) +
474+ std::to_string (xcurrent) + std::string (" ][stepsize =" ) +
475+ std::to_string (hslo) + std::string (" ]" )
476+ );
477+ } else if (unlikely (log_level == LogLevel::INFO)) {
478+ info.logger ().template log <LogLevel::INFO>(
479+ std::string (" [Oscillatory step][x = " ) +
480+ std::to_string (xcurrent) + std::string (" ][stepsize =" ) +
481+ std::to_string (hslo) + std::string (" ]" )
482+ );
468483 }
469484 while (!success) {
470485 std::tie (success, y, dy, err, y_eval, dy_eval, cheb_N)
471486 = nonosc_step (info, xcurrent, hslo, yprev, dyprev, eps);
472- // std::cout << "Attempted nonosc step with hslo = " << hslo << ", successful? " << success << std::endl;
487+ solver_counts[get_idx (LogInfo::CHEBSTEP)].second ++;
488+ solver_counts[get_idx (LogInfo::LS)].second += cheb_N + 1 ;
473489 steptype = 0 ;
474490 if (!success) {
475491 hslo *= Scalar{0.5 };
@@ -487,6 +503,13 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
487503 std::tie (dense_start, dense_size)
488504 = get_slice (x_eval, xcurrent, (xcurrent + h));
489505 if (dense_size != 0 ) {
506+ if (unlikely (log_level == LogLevel::INFO)) {
507+ info.logger ().template log <LogLevel::INFO>(
508+ std::string (" [Dense output][x_start = " ) +
509+ std::to_string (xcurrent) + std::string (" ][x_end =" ) +
510+ std::to_string (xcurrent + h) + std::string (" ]" )
511+ );
512+ }
490513 auto x_eval_map
491514 = Eigen::Map<vectord_t >(x_eval.data () + dense_start, dense_size);
492515 auto y_eval_map
@@ -564,11 +587,29 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
564587 }
565588 info.alloc_ .recover_memory ();
566589 }
567- #ifdef RICCATI_DEBUG
568- std::cout << " Total riccati steps: " << successes.size () << std::endl;
569- #endif
590+ if constexpr (!std::is_same_v<std::decay_t <decltype (info.logger ())>, EmptyLogger>) {
591+ if (unlikely (log_level == LogLevel::INFO)) {
592+ std::size_t riccati_steps = 0 ;
593+ std::size_t rk_steps = 0 ;
594+ for (auto & success : successes) {
595+ if (success) {
596+ riccati_steps++;
597+ } else {
598+ rk_steps++;
599+ }
600+ }
601+ info.logger ().template log <LogLevel::INFO>(
602+ std::string (" Total Steps = " ) + std::to_string (successes.size ())
603+ );
604+ for (auto && info_pair : solver_counts) {
605+ info.logger ().template log <LogLevel::INFO>(
606+ std::string (to_string (info_pair.first )) + " = " + std::to_string (info_pair.second )
607+ );
608+ }
609+ }
610+ }
570611 return std::make_tuple (xs, ys, dys, successes, phases, steptypes, yeval,
571- dyeval);
612+ dyeval, 1.0 );
572613}
573614
574615} // namespace riccati
0 commit comments