diff --git a/photon-targeting/src/main/native/cpp/photon/casadi_meme/wrap/casadi_wrapper.cpp b/photon-targeting/src/main/native/cpp/photon/casadi_meme/wrap/casadi_wrapper.cpp index c1bf665a70..b4a9d36760 100644 --- a/photon-targeting/src/main/native/cpp/photon/casadi_meme/wrap/casadi_wrapper.cpp +++ b/photon-targeting/src/main/native/cpp/photon/casadi_meme/wrap/casadi_wrapper.cpp @@ -52,6 +52,8 @@ #include "../generate/casadi_meme_9_tags_False.h" #include "../generate/casadi_meme_9_tags_True.h" +constexpr bool VERBOSE = false; + struct Problem { int numTags; bool headingFree; @@ -180,7 +182,8 @@ void IterativeRefinement(const Solver& solver, } template -wpi::expected do_optimization( +wpi::expected +do_optimization( int nTags, casadi_meme::CameraCalibration cameraCal, // Note that casadi is column major, apparently Eigen::Matrix robot2camera, @@ -190,18 +193,16 @@ wpi::expected do_optimization( point_observations) { if (field2points.cols() != (nTags * 4) || point_observations.cols() != (nTags * 4)) { - fmt::println("Got unexpected num cols!"); - return wpi::unexpected{false}; + if constexpr (VERBOSE) fmt::println("Got unexpected num cols!"); + // TODO find a new error code + return wpi::unexpected{ + sleipnir::SolverExitCondition::kNonfiniteInitialCostOrConstraints}; } -#if 0 + if constexpr (VERBOSE) { fmt::println("----------------------------------"); - fmt::println("Camera cal {} {} {} {}", - cameraCal.fx, - cameraCal.fy, - cameraCal.cx, - cameraCal.cy - ); + fmt::println("Camera cal {} {} {} {}", cameraCal.fx, cameraCal.fy, + cameraCal.cx, cameraCal.cy); fmt::println("{} tags", nTags); fmt::println("nstate {}", nState); @@ -210,11 +211,12 @@ wpi::expected do_optimization( std::cout << "field2pt:\n" << field2points << std::endl; std::cout << "observations:\n" << point_observations << std::endl; fmt::println("---------^^^^^^^^---------"); -#endif + } auto problemOpt = createProblem(nTags, nState == 3); if (!problemOpt) { - return wpi::unexpected{false}; + return wpi::unexpected{ + sleipnir::SolverExitCondition::kNonfiniteInitialCostOrConstraints}; } ProblemState pState{robot2camera, field2points, point_observations, cameraCal, *problemOpt}; @@ -228,30 +230,48 @@ wpi::expected do_optimization( FullStateMat x_full = x_guess; auto x = x_full.template head(); - // Sleipnir's delta_I caching algo from + // Sleipnir's delta_I caching algo and Newton.cpp inspiration from // https://github.com/SleipnirGroup/Sleipnir/blob/5af8519f268a8075e245bb7cd411a81e1598f521/src/optimization/RegularizedLDLT.hpp#L163 // licensed under BSD 3-Clause + /// The value of δ from the previous iteration - 1e-4 is a sane guess /// TUNABLE double δ = 1e-4 * 2.0; + constexpr double ERROR_TOL = 1e-4; + for (int iter = 0; iter < 100; iter++) { auto iter_start = wpi::Now(); - HessianMat H = pState.calculateHessJ(x_full); + // Check for diverging iterates + if (x_full.template lpNorm() > 1e20 || + !x_full.allFinite()) { + return wpi::unexpected{sleipnir::SolverExitCondition::kDivergingIterates}; + } + GradMat g = pState.calculateGradJ(x_full); + // If our previous step found an x such grad(J) is acceptable, we're done + if (g.norm() < ERROR_TOL) { + // Done! + fmt::println("{}: Exiting due to convergence (‖∇J‖={})", iter, g.norm()); + break; + } + + HessianMat H = pState.calculateHessJ(x_full); + /// Regularization. If the Hessian inertia is already OK, don't adjust auto H_ldlt = H.ldlt(); if (H_ldlt.info() != Eigen::Success) { std::cerr << "LDLT decomp failed! H=" << std::endl << H << std::endl; - return wpi::unexpected{false}; + return wpi::unexpected{sleipnir::SolverExitCondition::kLocallyInfeasible}; } // Make sure H is positive definite (all eigenvalues are > 0) int i_reg{0}; - HessianMat H_reg = H; // keep track of our regularized H TODO - is this valid? + HessianMat H_reg = + H; // keep track of our regularized H TODO - is this valid? if ((H_ldlt.vectorD().array() <= 0.0).any()) { // If δthe Hessian wasn't regularized in a previous iteration, start at a // small value of δ. Otherwise, attempt a δ half as big as the previous @@ -269,7 +289,8 @@ wpi::expected do_optimization( if (H_ldlt.info() != Eigen::Success) { std::cerr << "LDLT decomp failed! H=" << std::endl << H << std::endl; - return wpi::unexpected{false}; + return wpi::unexpected{ + sleipnir::SolverExitCondition::kLocallyInfeasible}; } // If our eigenvalues aren't positive definite, pick a new δ for next @@ -279,7 +300,8 @@ wpi::expected do_optimization( // If the Hessian perturbation is too high, report failure if (δ > 1e20) { - return wpi::unexpected{false}; + return wpi::unexpected{ + sleipnir::SolverExitCondition::kLocallyInfeasible}; } } else { // Done! @@ -290,57 +312,75 @@ wpi::expected do_optimization( } if (i_reg == MAX_REG_STEPS) { - return wpi::unexpected{false}; + return wpi::unexpected{ + sleipnir::SolverExitCondition::kLocallyInfeasible}; } } else { // std::printf("Already regularized\n"); } // Solve for p_x, and refine our solution - auto Hsolver = H.fullPivLu(); + auto Hsolver = H_ldlt; // H.fullPivLu(); StateMat p_x = Hsolver.solve(-g); - // IterativeRefinement(Hsolver, H_reg, p_x, -g); casadi_real old_cost = pState.calculateJ(x_full); - double alpha = 1.0; + constexpr double α_max = 1.0; + double alpha = α_max; // Iterate until our chosen trial_x decreases our cost int alpha_refinement{0}; + FullStateMat trial_x = x_full; for (alpha_refinement = 0; alpha_refinement < 100; alpha_refinement++) { - FullStateMat trial_x = x_full; trial_x.template head() = x + alpha * p_x; casadi_real new_cost = pState.calculateJ(trial_x); + // If f(xₖ + αpₖˣ) isn't finite, reduce step size immediately + if (!std::isfinite(new_cost)) { + // Reduce step size + alpha *= 0.5; + continue; + } + // Make sure we see an improvement if (new_cost < old_cost) { + // Step accepted - update x x = trial_x.template head(); break; } else { alpha *= 0.5; - } - } - if (g.norm() < 1e-8) { - // Done! - break; + // Safety factor for the minimal step size + constexpr double α_min_frac = 0.05; + constexpr double γConstraint = 1e-5; + + // If our step size shrank too much, report local infesibility + if (alpha < α_min_frac * γConstraint) { + return wpi::unexpected{ + sleipnir::SolverExitCondition::kLocallyInfeasible}; + } + } } auto iter_end = wpi::Now(); - fmt::println( - "{}: {} uS, ‖∇J‖={}, α={} ({} refinement steps), {} regularization " - "steps", - iter, iter_end - iter_start, g.norm(), alpha, alpha_refinement, i_reg); - fmt::println("∇J={}", g); - fmt::println("H={}", H); - fmt::println("p_x={}", p_x); - fmt::println("|Hp_x + ∇f|₂={}", (H * p_x + g).norm()); + if constexpr (VERBOSE) { + fmt::println( + "{}: {} uS, ‖∇J‖={}, α={} ({} refinement steps), {} regularization " + "steps", + iter, iter_end - iter_start, g.norm(), alpha, alpha_refinement, + i_reg); + // fmt::println("∇J={}", g); + // fmt::println("H={}", H); + // fmt::println("p_x={}", p_x); + // fmt::println("|Hp_x + ∇f|₂={}", (H * p_x + g).norm()); + } } - fmt::println("======================"); + if constexpr (VERBOSE) fmt::println("======================"); + return x_full; } -wpi::expected +wpi::expected casadi_meme::do_optimization_heading_free( int nTags, casadi_meme::CameraCalibration cameraCal, // Note that casadi is column major, apparently @@ -353,7 +393,7 @@ casadi_meme::do_optimization_heading_free( field2points, point_observations); } -wpi::expected +wpi::expected casadi_meme::do_optimization_heading_fixed( int nTags, casadi_meme::CameraCalibration cameraCal, // Note that casadi is column major, apparently diff --git a/photon-targeting/src/main/native/include/photon/casadi_meme/wrap/casadi_wrapper.h b/photon-targeting/src/main/native/include/photon/casadi_meme/wrap/casadi_wrapper.h index 7ab89d3418..838f898fa8 100644 --- a/photon-targeting/src/main/native/include/photon/casadi_meme/wrap/casadi_wrapper.h +++ b/photon-targeting/src/main/native/include/photon/casadi_meme/wrap/casadi_wrapper.h @@ -18,6 +18,7 @@ #pragma once #include +#include #include namespace casadi_meme { @@ -39,7 +40,8 @@ using RobotStateMat = Eigen::Matrix; * to this. The number of columns in field2points and point_observations just be * exactly 4x nTags. */ -wpi::expected do_optimization_heading_free( +wpi::expected +do_optimization_heading_free( int nTags, CameraCalibration cameraCal, // Note that casadi is column major, apparently Eigen::Matrix robot2camera, @@ -52,7 +54,8 @@ wpi::expected do_optimization_heading_free( * this. The number of columns in field2points and point_observations just be * exactly 4x nTags. */ -wpi::expected do_optimization_heading_fixed( +wpi::expected +do_optimization_heading_fixed( int nTags, CameraCalibration cameraCal, // Note that casadi is column major, apparently Eigen::Matrix robot2camera, diff --git a/photon-targeting/src/main/native/jni/CasadiMemeJNI.cpp b/photon-targeting/src/main/native/jni/CasadiMemeJNI.cpp index 952f63a708..7d23cdb556 100644 --- a/photon-targeting/src/main/native/jni/CasadiMemeJNI.cpp +++ b/photon-targeting/src/main/native/jni/CasadiMemeJNI.cpp @@ -87,7 +87,8 @@ Java_org_photonvision_jni_CasadiMeme_do_1optimization std::cout << "observations:\n" << pointObservationsMat << std::endl; #endif - wpi::expected result; + wpi::expected + result; if (headingFree) { result = casadi_meme::do_optimization_heading_free( nTags, cameraCal_, robot2cameraMat, xGuessMat, field2pointsMat,