@@ -66,10 +66,10 @@ namespace riccati {
6666 * standard discrete steps of the solver.
6767 */
6868template <typename SolverInfo, typename Scalar, typename Vec, typename YScalar>
69- inline auto osc_evolve (SolverInfo &&info, Scalar xi, Scalar xf,
70- YScalar yi, YScalar dyi ,
71- Scalar eps, Scalar epsilon_h, Scalar init_stepsize ,
72- Vec &&x_eval, bool hard_stop = false ) {
69+ inline auto osc_evolve (SolverInfo &&info, Scalar xi, Scalar xf, YScalar yi,
70+ YScalar dyi, Scalar eps, Scalar epsilon_h ,
71+ Scalar init_stepsize, Vec &&x_eval ,
72+ bool hard_stop = false ) {
7373 int sign = init_stepsize > 0 ? 1 : -1 ;
7474 using complex_t = promote_complex_t <Scalar>;
7575 using vectorc_t = vector_t <complex_t >;
@@ -92,12 +92,14 @@ inline auto osc_evolve(SolverInfo &&info, Scalar xi, Scalar xf,
9292 constexpr bool dense_output = compile_size_v<Vec> != 0 ;
9393 auto osc_ret = osc_step<dense_output>(info, omega_n, gamma_n, xi,
9494 init_stepsize, yi, dyi, eps);
95- auto make_osc_ret = [](auto && tup) {
96- return std::make_tuple (std::get<0 >(tup), std::get<1 >(tup), std::get<2 >(tup), std::get<3 >(tup));
95+ auto make_osc_ret = [](auto &&tup) {
96+ return std::make_tuple (std::get<0 >(tup), std::get<1 >(tup), std::get<2 >(tup),
97+ std::get<3 >(tup));
9798 };
9899 if (std::get<0 >(osc_ret) == 0 ) {
99- return std::make_tuple (false , xi, init_stepsize, make_osc_ret (osc_ret), vectorc_t (0 ),
100- vectorc_t (0 ), static_cast <Eigen::Index>(0 ),
100+ return std::make_tuple (false , xi, init_stepsize, make_osc_ret (osc_ret),
101+ vectorc_t (0 ), vectorc_t (0 ),
102+ static_cast <Eigen::Index>(0 ),
101103 static_cast <Eigen::Index>(0 ));
102104 } else {
103105 Eigen::Index dense_size = 0 ;
@@ -140,7 +142,8 @@ inline auto osc_evolve(SolverInfo &&info, Scalar xi, Scalar xf,
140142 // o and g written here
141143 auto h_next = choose_osc_stepsize (info, x_next, hosc_ini, epsilon_h);
142144 // TODO: CANNOT RETURN ARENA MATRIX FROM OSC_RET
143- return std::make_tuple (true , x_next, std::get<0 >(h_next), make_osc_ret (osc_ret), std::move (yeval),
145+ return std::make_tuple (true , x_next, std::get<0 >(h_next),
146+ make_osc_ret (osc_ret), std::move (yeval),
144147 std::move (dyeval), dense_start, dense_size);
145148 }
146149}
@@ -199,10 +202,10 @@ inline auto osc_evolve(SolverInfo &&info, Scalar xi, Scalar xf,
199202 * standard discrete steps of the solver.
200203 */
201204template <typename SolverInfo, typename Scalar, typename Vec, typename YScalar>
202- inline auto nonosc_evolve (SolverInfo &&info, Scalar xi, Scalar xf,
203- YScalar yi, YScalar dyi ,
204- Scalar eps, Scalar epsilon_h, Scalar init_stepsize ,
205- Vec &&x_eval, bool hard_stop = false ) {
205+ inline auto nonosc_evolve (SolverInfo &&info, Scalar xi, Scalar xf, YScalar yi,
206+ YScalar dyi, Scalar eps, Scalar epsilon_h ,
207+ Scalar init_stepsize, Vec &&x_eval ,
208+ bool hard_stop = false ) {
206209 using complex_t = promote_complex_t <Scalar>;
207210 using vectorc_t = vector_t <complex_t >;
208211 vectorc_t yeval;
@@ -244,12 +247,12 @@ inline auto nonosc_evolve(SolverInfo &&info, Scalar xi, Scalar xf,
244247 if (dense_size != 0 ) {
245248 auto x_eval_map = x_eval.segment (dense_start, dense_size);
246249 auto xi_scaled
247- = (xi + init_stepsize / 2
248- + (init_stepsize / 2 )
249- * std::get<2 >(info.chebyshev_ [std::get<6 >(nonosc_ret)])
250- .array ())
251- .matrix ()
252- .eval ();
250+ = (xi + init_stepsize / 2
251+ + (init_stepsize / 2 )
252+ * std::get<2 >(info.chebyshev_ [std::get<6 >(nonosc_ret)])
253+ .array ())
254+ .matrix ()
255+ .eval ();
253256 auto Linterp = interpolate (xi_scaled, x_eval_map, info.alloc_ );
254257 yeval = Linterp * std::get<4 >(nonosc_ret);
255258 dyeval = Linterp * std::get<5 >(nonosc_ret);
@@ -325,25 +328,26 @@ inline auto nonosc_evolve(SolverInfo &&info, Scalar xi, Scalar xf,
325328 * equations.
326329 * 6. @ref Eigen::Matrix<std::complex<Scalar>, -1, 1>: A vector containing the
327330 * interpolated solution at the specified `x_eval` The function returns these
328- * 7. An array of logging information for the solver process. This can be added with the
329- * info() from the `SolverInfo` object to get the total log information.
330- * vectors encapsulated in a standard tuple, providing comprehensive information
331- * about the solution process, including where the solution was evaluated, the
332- * values and derivatives of the solution at those points, success status of
333- * each step, type of each step, and phase angles where applicable.
331+ * 7. An array of logging information for the solver process. This can be added
332+ * with the info() from the `SolverInfo` object to get the total log
333+ * information. vectors encapsulated in a standard tuple, providing
334+ * comprehensive information about the solution process, including where the
335+ * solution was evaluated, the values and derivatives of the solution at those
336+ * points, success status of each step, type of each step, and phase angles
337+ * where applicable.
334338 */
335339template <typename SolverInfo, typename Scalar, typename Vec, typename YScalar>
336- inline auto evolve (SolverInfo &info, Scalar xi, Scalar xf,
337- YScalar yi, YScalar dyi,
338- Scalar eps, Scalar epsilon_h, Scalar init_stepsize,
339- Vec &&x_eval, bool hard_stop = false ,
340+ inline auto evolve (SolverInfo &info, Scalar xi, Scalar xf, YScalar yi,
341+ YScalar dyi, Scalar eps, Scalar epsilon_h,
342+ Scalar init_stepsize, Vec &&x_eval, bool hard_stop = false ,
340343 LogLevel log_level = riccati::LogLevel::ERROR) {
341344 using vectord_t = vector_t <Scalar>;
342345 Scalar direction = init_stepsize > 0 ? 1 : -1 ;
343346 if (init_stepsize * (xf - xi) < 0 ) {
344347 throw std::domain_error (
345348 " Direction of integration does not match stepsize sign,"
346- " adjust the direction or stepsize so that integration happens from xi to xf." );
349+ " adjust the direction or stepsize so that integration happens from xi "
350+ " to xf." );
347351 }
348352 // Check that yeval and x_eval are right size
349353 constexpr bool dense_output = compile_size_v<Vec> != 0 ;
@@ -432,8 +436,8 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
432436 * std::min (std::min (static_cast <Scalar>(1e8 ),
433437 static_cast <Scalar>(std::abs (omega_i / domega_i))),
434438 std::abs (gamma_i / dgamma_i));
435- info.logger ().template log <LogLevel::INFO>(log_level,
436- " [setup] " , std::pair (" hslo_ini" , hslo_ini));
439+ info.logger ().template log <LogLevel::INFO>(log_level, " [setup] " ,
440+ std::pair (" hslo_ini" , hslo_ini));
437441 if (hard_stop) {
438442 if (direction * (xi + hosc_ini) > direction * xf) {
439443 hosc_ini = xf - xi;
@@ -460,9 +464,10 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
460464 && direction * xcurrent < direction * xf) {
461465 if (unlikely (log_level >= LogLevel::INFO)) {
462466 [&]() RICCATI_COLD_PATH {
463- info.logger ().template log <LogLevel::INFO>(log_level,
464- std::pair (" std::abs(xcurrent - xf)" , std::abs (xcurrent - xf)),
465- std::pair (" direction * xcurrent" , direction * xcurrent));
467+ info.logger ().template log <LogLevel::INFO>(
468+ log_level,
469+ std::pair (" std::abs(xcurrent - xf)" , std::abs (xcurrent - xf)),
470+ std::pair (" direction * xcurrent" , direction * xcurrent));
466471 }();
467472 }
468473 Scalar phase{0.0 };
@@ -489,22 +494,20 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
489494 if constexpr (dense_output) {
490495 std::tie (success, y, dy, err, phase, un, d_un, a_pair)
491496 = osc_step<dense_output>(info, omega_n, gamma_n, xcurrent, hosc,
492- yprev, dyprev, eps);
497+ yprev, dyprev, eps);
493498 } else {
494- std::tie (success, y, dy, err, phase)
495- = osc_step<dense_output>(info, omega_n, gamma_n, xcurrent, hosc,
496- yprev, dyprev, eps);
499+ std::tie (success, y, dy, err, phase) = osc_step<dense_output>(
500+ info, omega_n, gamma_n, xcurrent, hosc, yprev, dyprev, eps);
497501 }
498502 steptype = 1 ;
499503 solver_counts[get_idx (LogInfo::RICCSTEP)].second ++;
500504 }
501505 if (unlikely (log_level >= LogLevel::INFO)) {
502506 [&]() RICCATI_COLD_PATH {
503- info.logger ().template log <LogLevel::INFO>(log_level,
504- std::pair (" oscillatory step" , (success ? " success" : " failure" )),
505- std::pair (" x" , xcurrent),
506- std::pair (" stepsize" , hosc)
507- );
507+ info.logger ().template log <LogLevel::INFO>(
508+ log_level,
509+ std::pair (" oscillatory step" , (success ? " success" : " failure" )),
510+ std::pair (" x" , xcurrent), std::pair (" stepsize" , hosc));
508511 }();
509512 }
510513 } else {
@@ -518,14 +521,14 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
518521 steptype = 0 ;
519522 if (unlikely (log_level >= LogLevel::INFO)) {
520523 [&]() RICCATI_COLD_PATH {
521- info.logger ().template log <LogLevel::INFO>(log_level,
522- std::pair ( " nonoscillatory step " , (success ? " success " : " failures " )) ,
523- std::pair (" x " , xcurrent) ,
524- std::pair ( " stepsize " , hslo ),
525- std::pair (" y " , y. real () ),
526- std::pair (" dy" , dy.real ()),
527- std::pair (" err" , std::real (err)),
528- std::pair (" cheb_steps" , cheb_N));
524+ info.logger ().template log <LogLevel::INFO>(
525+ log_level ,
526+ std::pair (" nonoscillatory step " ,
527+ (success ? " success " : " failures " ) ),
528+ std::pair (" x " , xcurrent), std::pair ( " stepsize " , hslo ),
529+ std::pair ( " y " , y. real ()), std::pair (" dy" , dy.real ()),
530+ std::pair (" err" , std::real (err)),
531+ std::pair (" cheb_steps" , cheb_N));
529532 }();
530533 }
531534 if (!success) {
@@ -546,10 +549,9 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
546549 if (dense_size != 0 ) {
547550 if (unlikely (log_level >= LogLevel::INFO)) {
548551 [&]() RICCATI_COLD_PATH {
549- info.logger ().template log <LogLevel::INFO>(log_level,
550- " dense output" ,
551- std::pair (" x_start" , xcurrent),
552- std::pair (" x_end" , xcurrent + h));
552+ info.logger ().template log <LogLevel::INFO>(
553+ log_level, " dense output" , std::pair (" x_start" , xcurrent),
554+ std::pair (" x_end" , xcurrent + h));
553555 }();
554556 }
555557 auto x_eval_map
@@ -583,17 +585,13 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
583585 }
584586 }
585587 // Finish appending and ending conditions
586- if (unlikely (log_level >= LogLevel::INFO)) {
587- [&]() RICCATI_COLD_PATH {
588- info.logger ().template log <LogLevel::INFO>(log_level,
589- " step completed" ,
590- std::pair (" y" , y),
591- std::pair (" dy" , dy),
592- std::pair (" x_next" , xcurrent + h),
593- std::pair (" phase" , phase),
594- std::pair (" steptype" , steptype)
595- );
596- }();
588+ if (unlikely (log_level >= LogLevel::INFO)) {
589+ [&]() RICCATI_COLD_PATH {
590+ info.logger ().template log <LogLevel::INFO>(
591+ log_level, " step completed" , std::pair (" y" , y), std::pair (" dy" , dy),
592+ std::pair (" x_next" , xcurrent + h), std::pair (" phase" , phase),
593+ std::pair (" steptype" , steptype));
594+ }();
597595 }
598596 ys.push_back (y);
599597 dys.push_back (dy);
@@ -641,32 +639,31 @@ inline auto evolve(SolverInfo &info, Scalar xi, Scalar xf,
641639 }
642640 info.alloc_ .recover_memory ();
643641 }
644- if constexpr (!std::is_same_v<std::decay_t <decltype (info.logger ())>, EmptyLogger>) {
642+ if constexpr (!std::is_same_v<std::decay_t <decltype (info.logger ())>,
643+ EmptyLogger>) {
645644 if (unlikely (log_level >= LogLevel::INFO)) {
646645 info.logger ().template log <LogLevel::INFO>(log_level,
647- std::string (" x" , xcurrent));
646+ std::string (" x" , xcurrent));
648647 std::size_t riccati_steps = 0 ;
649648 std::size_t rk_steps = 0 ;
650- for (auto & success : successes) {
649+ for (auto & success : successes) {
651650 if (success) {
652651 riccati_steps++;
653652 } else {
654653 rk_steps++;
655654 }
656655 }
657- info.logger ().template log <LogLevel::INFO>(log_level,
658- std::pair (" Total Steps" , successes.size ())
659- );
660- for (auto && info_pair : solver_counts) {
656+ info.logger ().template log <LogLevel::INFO>(
657+ log_level, std::pair (" Total Steps" , successes.size ()));
658+ for (auto &&info_pair : solver_counts) {
661659 info.logger ().template log <LogLevel::INFO>(log_level, info_pair);
662660 }
663661 }
664662 }
665663 return std::make_tuple (xs, ys, dys, successes, phases, steptypes, yeval,
666- dyeval, 1.0 );
664+ dyeval, 1.0 );
667665}
668666
669-
670667} // namespace riccati
671668
672669#endif
0 commit comments