Fixes2

.bcr/presubmit.yml

@@ -2,7 +2,7 @@ matrix:
   platform:
     - debian11
     - ubuntu2004
-    - macos
+    - macos_arm64
   bazel:
     - 7.x
     - 8.x

.github/workflows/ci.yml

@@ -32,7 +32,6 @@ jobs:
             libeigen3-dev \
             libgtest-dev
 
-
       - name: Set up Bazel 9
         uses: bazel-contrib/setup-bazel@0.19.0
         with:
@@ -54,16 +53,17 @@ jobs:
           reporter: github-pr-review
           flags: '--extensions=h,hpp,c,cpp,cc,cu,hh,ipp --filter=-build/include_order,-whitespace/indent_namespace --recursive include'
 
-      - name: Verify
+      - name: Build (strict warnings)
         run: |
-          bazel run verify
+          bazel build //... \
+            --copt=-Wall --copt=-Wextra --copt=-Wpedantic \
+            --copt=-Wdouble-promotion --copt=-Werror
 
-      - name: Build
+      - name: Test
         run: |
-          bazel run simple
-          bazel run constrained_simple
-          bazel run constrained_simple2
-          bazel run debug
+          bazel test //... \
+            --copt=-Wall --copt=-Wextra --copt=-Wpedantic \
+            --copt=-Wdouble-promotion --copt=-Werror
 
       - name: Configure and Build with CMake
         run: |
@@ -72,3 +72,29 @@ jobs:
           cmake .. -DCMAKE_BUILD_TYPE=Release
           cmake --build .
           ctest --output-on-failure
+
+  sanitizers:
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: Checkout repository
+        uses: actions/checkout@v4
+
+      - name: Install System Dependencies
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y g++ libeigen3-dev
+
+      - name: Set up Bazel 9
+        uses: bazel-contrib/setup-bazel@0.19.0
+        with:
+          bazelisk-cache: true
+          disk-cache: ${{ github.workflow }}-san
+          repository-cache: true
+
+      - name: Test with ASAN + UBSAN
+        run: |
+          bazel test //... \
+            --copt=-fsanitize=address,undefined \
+            --copt=-fno-omit-frame-pointer \
+            --linkopt=-fsanitize=address,undefined

BUILD

@@ -1,5 +1,5 @@
 load("//:generator.bzl", "build_example", "build_test")
-load("@rules_cc//cc:defs.bzl", "cc_library")
+load("@rules_cc//cc:defs.bzl", "cc_library", "cc_test")
 
 # Shared header library used by the SVM examples.  Exposes the
 # embedded Iris versicolor-vs-virginica dataset and a classification
@@ -23,6 +23,16 @@ build_test("hager_zhang_test")
 build_test("augmented_lagrangian_test")
 build_test("trust_region_newton_test")
 
+cc_test(
+    name = "float_scalar_test",
+    srcs = ["src/test/float_scalar_test.cc"],
+    copts = ["-std=c++17", "-Wall", "-Wextra", "-Wdouble-promotion", "-Werror"],
+    deps = [
+        "//include:cppoptlib",
+        "@eigen//:eigen",
+    ],
+)
+
 build_example("svm_primal_lbfgs", extra_deps = [":iris_data"])
 build_example("svm_primal_al", extra_deps = [":iris_data"])
 build_example("svm_dual_lbfgsb", extra_deps = [":iris_data"])

include/cppoptlib/function_base.h

@@ -100,8 +100,9 @@ struct FunctionCRTP : public FunctionInterface<TScalar, TMode, TDimension> {
   using MatrixType = Eigen::Matrix<TScalar, Dimension, Dimension>;
   static constexpr DifferentiabilityMode Differentiability = TMode;
 
-  virtual TScalar operator()(const VectorType& x, VectorType* grad = nullptr,
-                             MatrixType* hess = nullptr) const override {
+  virtual TScalar operator()(
+      const VectorType& x, [[maybe_unused]] VectorType* grad = nullptr,
+      [[maybe_unused]] MatrixType* hess = nullptr) const override {
     if constexpr (TMode == DifferentiabilityMode::None) {
       return static_cast<const Derived*>(this)->operator()(x);
     } else if constexpr (TMode == DifferentiabilityMode::First) {

include/cppoptlib/function_expressions.h

@@ -54,8 +54,9 @@ struct ConstExpression : public FunctionInterface<TScalar, TMode, TDimension> {
 
   TScalar c;
   explicit ConstExpression(TScalar c_) : c(c_) {}
-  virtual TScalar operator()(const VectorType& x, VectorType* grad = nullptr,
-                             MatrixType* hess = nullptr) const override {
+  virtual TScalar operator()(
+      const VectorType& x, VectorType* grad = nullptr,
+      [[maybe_unused]] MatrixType* hess = nullptr) const override {
     (void)x;
     if (grad) {
       *grad = VectorType::Zero(x.size());
@@ -164,8 +165,9 @@ struct SubExpression
   F f;
   G g;
   SubExpression(const F& f_, const G& g_) : f(f_), g(g_) {}
-  virtual ScalarType operator()(const VectorType& x, VectorType* grad = nullptr,
-                                MatrixType* hess = nullptr) const override {
+  virtual ScalarType operator()(
+      const VectorType& x, VectorType* grad = nullptr,
+      [[maybe_unused]] MatrixType* hess = nullptr) const override {
     if constexpr (TMode == DifferentiabilityMode::None) {
       return f(x) - g(x);
     } else if constexpr (TMode == DifferentiabilityMode::First) {
@@ -214,8 +216,9 @@ struct MulExpression : public FunctionInterface<TScalar, TMode, F::Dimension> {
   TScalar c;
   F f;
   MulExpression(const TScalar& c_, const F& f_) : c(c_), f(f_) {}
-  virtual ScalarType operator()(const VectorType& x, VectorType* grad = nullptr,
-                                MatrixType* hess = nullptr) const override {
+  virtual ScalarType operator()(
+      const VectorType& x, VectorType* grad = nullptr,
+      [[maybe_unused]] MatrixType* hess = nullptr) const override {
     if (c == TScalar(0)) {
       if (grad) {
         *grad = VectorType::Zero(x.size());
@@ -278,8 +281,9 @@ struct ProdExpression
 
   ProdExpression(const F& f_, const G& g_) : f(f_), g(g_) {}
 
-  virtual ScalarType operator()(const VectorType& x, VectorType* grad = nullptr,
-                                MatrixType* hess = nullptr) const override {
+  virtual ScalarType operator()(
+      const VectorType& x, VectorType* grad = nullptr,
+      [[maybe_unused]] MatrixType* hess = nullptr) const override {
     if constexpr (TMode == DifferentiabilityMode::None) {
       return f(x) * g(x);
     } else if constexpr (TMode == DifferentiabilityMode::First) {
@@ -329,8 +333,9 @@ struct MinZeroExpression
   F f;
   explicit MinZeroExpression(const F& f_) : f(f_) {}
 
-  virtual ScalarType operator()(const VectorType& x, VectorType* grad = nullptr,
-                                MatrixType* hess = nullptr) const override {
+  virtual ScalarType operator()(
+      const VectorType& x, VectorType* grad = nullptr,
+      [[maybe_unused]] MatrixType* hess = nullptr) const override {
     if constexpr (Differentiability == DifferentiabilityMode::None) {
       ScalarType val = f(x);
       // Use ConstExpression to return zero (with zero derivatives) if inactive.
@@ -372,8 +377,9 @@ struct MaxZeroExpression
   F f;
   explicit MaxZeroExpression(const F& f_) : f(f_) {}
 
-  virtual ScalarType operator()(const VectorType& x, VectorType* grad = nullptr,
-                                MatrixType* hess = nullptr) const override {
+  virtual ScalarType operator()(
+      const VectorType& x, VectorType* grad = nullptr,
+      [[maybe_unused]] MatrixType* hess = nullptr) const override {
     if constexpr (Differentiability == DifferentiabilityMode::None) {
       ScalarType val = f(x);
       return (val <= 0)

include/cppoptlib/linesearch/armijo.h

@@ -45,19 +45,19 @@ class Armijo {
   static ScalarType Search(const VectorType& x,
                            const VectorType& search_direction,
                            const FunctionType& function,
-                           const ScalarType alpha_init = 1.0) {
+                           const ScalarType alpha_init = ScalarType(1)) {
     constexpr ScalarType c = ScalarType(0.2);
     constexpr ScalarType rho = ScalarType(0.9);
     ScalarType alpha = alpha_init;
     VectorType gradient;
     const ScalarType f_in = function(x, &gradient);
-    ScalarType f = function(x + alpha * search_direction);
+    ScalarType f = function((x + alpha * search_direction).eval());
     const ScalarType cache = c * gradient.dot(search_direction);
     constexpr ScalarType alpha_min = ScalarType(1e-8);
 
     while (f > f_in + alpha * cache && alpha > alpha_min) {
       alpha *= rho;
-      f = function(x + alpha * search_direction);
+      f = function((x + alpha * search_direction).eval());
     }
 
     return alpha;
@@ -84,19 +84,19 @@ class Armijo<FunctionType, 2> {
                            const FunctionType& function) {
     constexpr ScalarType c = ScalarType(0.2);
     constexpr ScalarType rho = ScalarType(0.9);
-    ScalarType alpha = 1.0;
+    ScalarType alpha = ScalarType(1);
     VectorType gradient;
     MatrixType hessian;
     const ScalarType f_in = function(x, &gradient, &hessian);
-    ScalarType f = function(x + alpha * search_direction);
+    ScalarType f = function((x + alpha * search_direction).eval());
     const ScalarType cache =
         c * gradient.dot(search_direction) + ScalarType(0.5) * c * c *
                                                  search_direction.transpose() *
                                                  hessian * search_direction;
 
     while (f > f_in + alpha * cache) {
       alpha *= rho;
-      f = function(x + alpha * search_direction);
+      f = function((x + alpha * search_direction).eval());
     }
     return alpha;
   }

include/cppoptlib/solver/lbfgsb.h

@@ -208,7 +208,7 @@ class Lbfgsb
 
     // STEP 6: Only update if positive curvature (s'*y > 0)
     const ScalarType sTy = new_s.dot(new_y);
-    if (sTy > 1e-7 * new_y.squaredNorm()) {
+    if (sTy > ScalarType(1e-7) * new_y.squaredNorm()) {
       if (y_history_.cols() < m) {
         y_history_.conservativeResize(dim_, y_history_.cols() + 1);
         s_history_.conservativeResize(dim_, s_history_.cols() + 1);
@@ -439,7 +439,7 @@ class Lbfgsb
     assert(du.rows() == n);
 
     for (unsigned int i = 0; i < n; i++) {
-      if (std::abs(du(i)) < 1e-7) {
+      if (std::abs(du(i)) < ScalarType(1e-7)) {
         continue;
       } else if (du(i) > 0) {
         alphastar = std::min<ScalarType>(

include/cppoptlib/utils/derivatives.h

@@ -255,12 +255,12 @@ template <class FunctionType>
 bool IsGradientCorrect(const FunctionType& function,
                        const typename FunctionType::VectorType& x0,
                        int accuracy = 3) {
-  constexpr float tolerance = 1e-2;
-
   using ScalarType = typename FunctionType::ScalarType;
   using VectorType = typename FunctionType::VectorType;
   using index_t = typename VectorType::Index;
 
+  constexpr ScalarType tolerance = ScalarType(1e-2);
+
   const index_t D = x0.rows();
   VectorType actual_gradient;
   function(x0, &actual_gradient);
@@ -283,13 +283,13 @@ template <class FunctionType>
 bool IsHessianCorrect(const FunctionType& function,
                       const typename FunctionType::VectorType& x0,
                       int accuracy = 3) {
-  constexpr float tolerance = 1e-1;
-
   using ScalarType = typename FunctionType::ScalarType;
   using MatrixType = typename FunctionType::MatrixType;
   using VectorType = typename FunctionType::VectorType;
   using index_t = typename VectorType::Index;
 
+  constexpr ScalarType tolerance = ScalarType(1e-1);
+
   const index_t D = x0.rows();
 
   MatrixType actual_hessian;

src/examples/linear_regression.cc

@@ -72,6 +72,7 @@ int main() {
 
   const auto initial_state = cppoptlib::function::FunctionState(x);
   auto [solution, solver_state] = solver.Minimize(f, initial_state);
+  (void)solver_state;
   std::cout << "argmin " << solution.x.transpose() << std::endl;
 
   // Or model it as a augmented Lagrangian
@@ -97,6 +98,7 @@ int main() {
   // Run the agumented solver.
   auto [aug_solution, aug_solver_state] = aug_solver.Minimize(l_state);
   std::cout << "argmin " << aug_solution.x.transpose() << std::endl;
+  (void)aug_solver_state;
 
   return 0;
 }

src/test/augmented_lagrangian_test.cc

@@ -507,6 +507,7 @@ TEST(AugmentedLagrangianKKT, EqualityOnlyQuadratic) {
   x0 << 5.0, 5.0;
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 1, 0, 1.0);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   // Primal optimum.
   EXPECT_NEAR(1.0, solution.x[0], kkt_primal_tolerance);
@@ -558,6 +559,7 @@ TEST(AugmentedLagrangianKKT, InequalityActiveRecoversMultiplier) {
   x0 << 5.0, 5.0;
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 0, 1, 1.0);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   EXPECT_NEAR(1.0, solution.x[0], kkt_primal_tolerance);
   EXPECT_NEAR(0.0, solution.x[1], kkt_primal_tolerance);
@@ -601,6 +603,7 @@ TEST(AugmentedLagrangianKKT, BothEqualityAndInequalityActive) {
   x0 << 1.0, 1.0;
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 1, 1, 1.0);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   EXPECT_NEAR(0.5, solution.x[0], kkt_primal_tolerance);
   EXPECT_NEAR(1.5, solution.x[1], kkt_primal_tolerance);
@@ -643,6 +646,7 @@ TEST(AugmentedLagrangianOuter, FeasibleStartConvergesImmediately) {
   x0 << 0.0, 0.0;
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 1, 0, 1.0);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   EXPECT_NEAR(0.0, solution.x[0], kkt_primal_tolerance);
   EXPECT_NEAR(0.0, solution.x[1], kkt_primal_tolerance);
@@ -676,6 +680,7 @@ TEST(AugmentedLagrangianOuter, NoConstraintsIsUnconstrained) {
   x0 << 5.0, 5.0;
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 0, 0, 1.0);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   EXPECT_NEAR(0.0, solution.x[0], kkt_primal_tolerance);
   EXPECT_NEAR(0.0, solution.x[1], kkt_primal_tolerance);
@@ -711,6 +716,7 @@ TEST(AugmentedLagrangianOuter, PenaltyHoldsFlatOnFeasibleProblem) {
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 1, 0,
                                                              initial_penalty);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   // The conditional schedule should never fire on a trivially feasible
   // problem.  The penalty is pinned exactly at its initial value.
@@ -748,6 +754,7 @@ TEST(AugmentedLagrangianOuter, PenaltyGrowthCanBeDisabled) {
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 1, 0,
                                                              initial_penalty);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   EXPECT_EQ(initial_penalty, solution.penalty_state.penalty);
 }
@@ -781,6 +788,7 @@ TEST(AugmentedLagrangianOuter, PenaltyGrowsOnlyWhileViolationLags) {
   x0 << 5.0, 5.0;
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 1, 0, 1.0);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   constexpr double penalty_upper_bound = 1e4;
   EXPECT_LE(solution.penalty_state.penalty, penalty_upper_bound);
@@ -1051,6 +1059,7 @@ TEST(AugmentedLagrangianNonConvex, Hs024TriangleEscapesSpuriousOrigin) {
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(
       x0, /*num_eq=*/0, /*num_ineq=*/3, /*penalty=*/0.0);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   const double f_final = objective(solution.x);
 
@@ -1124,6 +1133,7 @@ TEST(AugmentedLagrangianNonConvex, Hs029EllipseEscapesOrigin) {
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(
       x0, /*num_eq=*/0, /*num_ineq=*/1, /*penalty=*/0.0);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   const double f_final = objective(solution.x);
   // True optimum for `min -x0*x1 s.t. 48 - x0^2 - 2 x1^2 >= 0` is at
@@ -1171,6 +1181,7 @@ TEST(AugmentedLagrangianOuter, KktStationarityReportedOnFinishedState) {
   x0 << 5.0, 5.0;
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 1, 0, 1.0);
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   ASSERT_EQ(cppoptlib::solver::Status::Finished, progress.status);
   // The outer-loop KKT tolerance defaults match the primal one to a
@@ -1197,7 +1208,7 @@ TEST(AugmentedLagrangianOuter, KktStationarityReportedOnFinishedState) {
 // `progress.status == Finished` fails.
 TEST(AugmentedLagrangianBoxInterface, BoxPinnedOptimumStopsOnKkt) {
   using cppoptlib::function::FunctionExpr;
-  using VectorType = Eigen::Matrix<double, 2, 1>;
+  using VectorType [[maybe_unused]] = Eigen::Matrix<double, 2, 1>;
 
   class Rosenbrock : public cppoptlib::function::FunctionCRTP<
                          Rosenbrock, double,
@@ -1261,6 +1272,7 @@ TEST(AugmentedLagrangianBoxInterface, BoxPinnedOptimumStopsOnKkt) {
   cppoptlib::solver::AugmentedLagrangeState<double, 2> state(x0, 0, 2, 0.0);
 
   auto [solution, progress] = solver.Minimize(state);
+  (void)progress;
 
   // The outer loop must reach `Finished` status on its own -- an
   // `IterationLimit` exit would signal regression of the projected-