Switch to CPM; support 3rd-order tensor jacobians

Author: zfergus

CMakeLists.txt

@@ -35,6 +35,8 @@ endif()
 list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/finitediff/")
 list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/recipes/")
 
+include(finitediff_cpm_cache)
+
 ################################################################################
 # Finite Diff Library
 ################################################################################

README.md

@@ -43,7 +43,7 @@ The library provides three main functions `finite_gradient`, `finite_jacobian`,
 
 ```c++
 void finite_gradient(
-    const Eigen::VectorXd& x,
+    const Eigen::Ref<const Eigen::VectorXd>& x,
     const std::function<double(const Eigen::VectorXd&)>& f,
     Eigen::VectorXd& grad,
     const AccuracyOrder accuracy = SECOND,
@@ -56,8 +56,8 @@ The `finite_gradient` function computes the [gradient](https://en.wikipedia.org/
 
 ```c++
 void finite_jacobian(
-    const Eigen::VectorXd& x,
-    const std::function<Eigen::VectorXd(const Eigen::VectorXd&)>& f,
+    const Eigen::Ref<const Eigen::VectorXd>& x,
+    const std::function<Eigen::MatrixXd(const Eigen::VectorXd&)>& f,
     Eigen::MatrixXd& jac,
     const AccuracyOrder accuracy = SECOND,
     const double eps = 1.0e-8);
@@ -70,7 +70,7 @@ The `finite_jacobian` function computes the [Jacobian](https://en.wikipedia.org/
 
 ```c++
 void finite_hessian(
-    const Eigen::VectorXd& x,
+    const Eigen::Ref<const Eigen::VectorXd>& x,
     const std::function<double(const Eigen::VectorXd&)>& f,
     Eigen::MatrixXd& hess,
     const AccuracyOrder accuracy = SECOND,

cmake/finitediff/finitediff_cpm_cache.cmake

@@ -0,0 +1,24 @@
+#
+# Copyright 2021 Adobe. All rights reserved.
+# This file is licensed to you under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License. You may obtain a copy
+# of the License at http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software distributed under
+# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
+# OF ANY KIND, either express or implied. See the License for the specific language
+# governing permissions and limitations under the License.
+#
+
+if(DEFINED ENV{CPM_SOURCE_CACHE})
+    set(CPM_SOURCE_CACHE_DEFAULT $ENV{CPM_SOURCE_CACHE})
+else()
+    # Set CPM cache folder if unset
+    file(REAL_PATH "~/.cache/CPM" CPM_SOURCE_CACHE_DEFAULT EXPAND_TILDE)
+endif()
+
+set(CPM_SOURCE_CACHE
+    ${CPM_SOURCE_CACHE_DEFAULT}
+    CACHE PATH "Directory to download CPM dependencies"
+)
+message(STATUS "Using CPM cache folder: ${CPM_SOURCE_CACHE}")

cmake/recipes/CPM.cmake

@@ -0,0 +1,33 @@
+set(CPM_DOWNLOAD_VERSION 0.40.2)
+
+if(CPM_SOURCE_CACHE)
+  set(CPM_DOWNLOAD_LOCATION "${CPM_SOURCE_CACHE}/cpm/CPM_${CPM_DOWNLOAD_VERSION}.cmake")
+elseif(DEFINED ENV{CPM_SOURCE_CACHE})
+  set(CPM_DOWNLOAD_LOCATION "$ENV{CPM_SOURCE_CACHE}/cpm/CPM_${CPM_DOWNLOAD_VERSION}.cmake")
+else()
+  set(CPM_DOWNLOAD_LOCATION "${CMAKE_BINARY_DIR}/cmake/CPM_${CPM_DOWNLOAD_VERSION}.cmake")
+endif()
+
+# Expand relative path. This is important if the provided path contains a tilde (~)
+get_filename_component(CPM_DOWNLOAD_LOCATION ${CPM_DOWNLOAD_LOCATION} ABSOLUTE)
+
+function(download_cpm)
+  message(STATUS "Downloading CPM.cmake to ${CPM_DOWNLOAD_LOCATION}")
+  file(DOWNLOAD
+       https://github.com/cpm-cmake/CPM.cmake/releases/download/v${CPM_DOWNLOAD_VERSION}/CPM.cmake
+       ${CPM_DOWNLOAD_LOCATION}
+  )
+endfunction()
+
+if(NOT (EXISTS ${CPM_DOWNLOAD_LOCATION}))
+  download_cpm()
+else()
+  # resume download if it previously failed
+  file(READ ${CPM_DOWNLOAD_LOCATION} check)
+  if("${check}" STREQUAL "")
+    download_cpm()
+  endif()
+  unset(check)
+endif()
+
+include(${CPM_DOWNLOAD_LOCATION})

cmake/recipes/catch2.cmake

@@ -1,25 +1,18 @@
-#
-# Copyright 2020 Adobe. All rights reserved.
-# This file is licensed to you under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License. You may obtain a copy
-# of the License at http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under
-# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
-# OF ANY KIND, either express or implied. See the License for the specific language
-# governing permissions and limitations under the License.
-#
+# Catch2 (https://github.com/catchorg/Catch2)
+# License: BSL-1.0
 if(TARGET Catch2::Catch2)
     return()
 endif()
 
 message(STATUS "Third-party: creating target 'Catch2::Catch2'")
 
-include(FetchContent)
-FetchContent_Declare(
-    catch2
-    GIT_REPOSITORY https://github.com/catchorg/Catch2.git
-    GIT_TAG v3.0.1
-    GIT_SHALLOW TRUE
-)
-FetchContent_MakeAvailable(catch2)
+option(CATCH_CONFIG_CPP17_STRING_VIEW "Enable support for std::string_view" ON)
+option(CATCH_INSTALL_DOCS "Install documentation alongside library" OFF)
+option(CATCH_INSTALL_EXTRAS "Install extras alongside library" OFF)
+
+include(CPM)
+CPMAddPackage("gh:catchorg/[email protected]")
+
+# Folder name for IDE
+set_target_properties(Catch2 PROPERTIES FOLDER "ThirdParty/Catch2")
+set_target_properties(Catch2WithMain PROPERTIES FOLDER "ThirdParty/Catch2")

cmake/recipes/eigen.cmake

@@ -1,53 +1,39 @@
-#
-# Copyright 2020 Adobe. All rights reserved.
-# This file is licensed to you under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License. You may obtain a copy
-# of the License at http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under
-# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
-# OF ANY KIND, either express or implied. See the License for the specific language
-# governing permissions and limitations under the License.
-#
+# Eigen (https://gitlab.com/libeigen/eigen)
+# License: MPL 2.0
 if(TARGET Eigen3::Eigen)
     return()
 endif()
 
 option(EIGEN_WITH_MKL "Use Eigen with MKL" OFF)
+option(EIGEN_DONT_VECTORIZE "Disable Eigen vectorization" OFF)
+option(EIGEN_MPL2_ONLY "Enable Eigen MPL2 license only" OFF)
 
-if(EIGEN_ROOT)
-    message(STATUS "Third-party: creating target 'Eigen3::Eigen' for external path: ${EIGEN_ROOT}")
-    set(EIGEN_INCLUDE_DIRS ${EIGEN_ROOT})
-else()
-    message(STATUS "Third-party: creating target 'Eigen3::Eigen'")
+message(STATUS "Third-party: creating target 'Eigen3::Eigen'")
 
-    include(FetchContent)
-    FetchContent_Declare(
-        eigen
-        GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
-        GIT_TAG tags/3.4.0
-        GIT_SHALLOW TRUE
-    )
-    FetchContent_GetProperties(eigen)
-    if(NOT eigen_POPULATED)
-        FetchContent_Populate(eigen)
-    endif()
-    set(EIGEN_INCLUDE_DIRS ${eigen_SOURCE_DIR})
-
-    install(DIRECTORY ${EIGEN_INCLUDE_DIRS}/Eigen
-        DESTINATION include
-    )
-endif()
+include(CPM)
+CPMAddPackage(
+    NAME eigen
+    GITLAB_REPOSITORY libeigen/eigen
+    GIT_TAG 3.4.0
+    DOWNLOAD_ONLY YES
+)
 
 add_library(Eigen3_Eigen INTERFACE)
 add_library(Eigen3::Eigen ALIAS Eigen3_Eigen)
 
 include(GNUInstallDirs)
 target_include_directories(Eigen3_Eigen SYSTEM INTERFACE
-    $<BUILD_INTERFACE:${EIGEN_INCLUDE_DIRS}>
+    $<BUILD_INTERFACE:${eigen_SOURCE_DIR}>
     $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
 )
-# target_compile_definitions(Eigen3_Eigen INTERFACE EIGEN_MPL2_ONLY)
+
+if(EIGEN_MPL2_ONLY)
+  target_compile_definitions(Eigen3_Eigen INTERFACE EIGEN_MPL2_ONLY)
+endif()
+
+if(EIGEN_DONT_VECTORIZE)
+  target_compile_definitions(Eigen3_Eigen INTERFACE EIGEN_DONT_VECTORIZE)
+endif()
 
 if(EIGEN_WITH_MKL)
     # TODO: Checks that, on 64bits systems, `mkl::mkl` is using the LP64 interface
@@ -58,6 +44,27 @@ if(EIGEN_WITH_MKL)
         EIGEN_USE_MKL_ALL
         EIGEN_USE_LAPACKE_STRICT
     )
+elseif(APPLE)
+    find_package(BLAS REQUIRED)
+    find_library(LAPACKE lapacke PATHS
+        "/opt/local/lib/lapack"
+        "/opt/homebrew/opt/lapack/lib"
+    )
+    if (NOT LAPACKE)
+        # BLAS should be available on macOS, but LAPACKE might not be
+        message(WARNING "LAPACKE library not found (required for EIGEN_USE_LAPACKE on macOS)! "
+            "Perhaps you need to install it (e.g., brew install lapack). "
+            "Eigen will be built without LAPACKE support.")
+    else()
+        message(STATUS "Found BLAS and LAPACKE. Enabling Eigen LAPACKE support.")
+        target_link_libraries(Eigen3_Eigen INTERFACE
+            ${BLAS_LIBRARIES} ${LAPACKE}
+        )
+        target_compile_definitions(Eigen3_Eigen INTERFACE
+            EIGEN_USE_BLAS
+            EIGEN_USE_LAPACKE_STRICT
+        )
+    endif()
 endif()
 
 # On Windows, enable natvis files to improve debugging experience
@@ -68,6 +75,6 @@ endif()
 # Install rules
 set(CMAKE_INSTALL_DEFAULT_COMPONENT_NAME eigen)
 set_target_properties(Eigen3_Eigen PROPERTIES EXPORT_NAME Eigen)
-install(DIRECTORY ${EIGEN_INCLUDE_DIRS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+install(DIRECTORY ${eigen_SOURCE_DIR} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
 install(TARGETS Eigen3_Eigen EXPORT Eigen_Targets)
-install(EXPORT Eigen_Targets DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/eigen NAMESPACE Eigen3::)
+install(EXPORT Eigen_Targets DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/eigen NAMESPACE Eigen3::)

cmake/recipes/spdlog.cmake

@@ -1,35 +1,21 @@
-#
-# Copyright 2020 Adobe. All rights reserved.
-# This file is licensed to you under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License. You may obtain a copy
-# of the License at http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under
-# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
-# OF ANY KIND, either express or implied. See the License for the specific language
-# governing permissions and limitations under the License.
-#
+# spdlog (https://github.com/gabime/spdlog)
+# License: MIT
 if(TARGET spdlog::spdlog)
     return()
 endif()
 
 message(STATUS "Third-party: creating target 'spdlog::spdlog'")
 
-include(FetchContent)
-FetchContent_Declare(
-    spdlog
-    GIT_REPOSITORY https://github.com/gabime/spdlog.git
-    GIT_TAG v1.10.0
-    GIT_SHALLOW TRUE
-)
-
 option(SPDLOG_INSTALL "Generate the install target" ON)
 set(CMAKE_INSTALL_DEFAULT_COMPONENT_NAME "spdlog")
-FetchContent_MakeAvailable(spdlog)
+
+include(CPM)
+CPMAddPackage("gh:gabime/[email protected]")
 
 set_target_properties(spdlog PROPERTIES POSITION_INDEPENDENT_CODE ON)
 
-set_target_properties(spdlog PROPERTIES FOLDER external)
+# Folder name for IDE
+set_target_properties(spdlog PROPERTIES FOLDER "ThirdParty")
 
 if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang" OR
    "${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")

src/finitediff.cpp

@@ -94,7 +94,7 @@ void finite_gradient(
 
 void finite_jacobian(
     const Eigen::Ref<const Eigen::VectorXd>& x,
-    const std::function<Eigen::VectorXd(const Eigen::VectorXd&)>& f,
+    const std::function<Eigen::MatrixXd(const Eigen::VectorXd&)>& f,
     Eigen::MatrixXd& jac,
     const AccuracyOrder accuracy,
     const double eps)
@@ -107,16 +107,24 @@ void finite_jacobian(
 
     const double denom = get_denominator(accuracy) * eps;
 
-    jac.setZero(f(x).rows(), x.rows());
+    // Call f once to get the size of the Jacobian
+    size_t f_rows, f_cols;
+    {
+        const Eigen::MatrixXd tmp = f(x);
+        f_rows = tmp.rows(), f_cols = tmp.cols();
+    }
+    jac.setZero(f_rows, f_cols * x.size());
 
+    // f: ℝ^n ↦ ℝ^{p×q} ⟹ ∇f: ℝ^n ↦ ℝ^{p×(qn)}
     Eigen::VectorXd x_mutable = x;
-    for (size_t i = 0; i < x.rows(); i++) {
+    for (size_t i = 0; i < x.size(); i++) {
         for (size_t ci = 0; ci < inner_steps; ci++) {
             x_mutable[i] += internal_coeffs[ci] * eps;
-            jac.col(i) += external_coeffs[ci] * f(x_mutable);
+            jac.middleCols(f_cols * i, f_cols) +=
+                external_coeffs[ci] * f(x_mutable);
             x_mutable[i] = x[i];
         }
-        jac.col(i) /= denom;
+        jac.middleCols(f_cols * i, f_cols) /= denom;
     }
 }
 

src/finitediff.hpp

@@ -39,27 +39,27 @@ void finite_gradient(
     const double eps = 1.0e-8);
 
 /**
- * @brief Compute the jacobian of a function using finite differences.
+ * @brief Compute the Jacobian of a function using finite differences.
  *
- * @param[in]  x         Point at which to compute the jacobian.
- * @param[in]  f         Compute the jacobian of this function.
- * @param[out] jac       Computed jacobian.
+ * @param[in]  x         Point at which to compute the Jacobian.
+ * @param[in]  f         Compute the Jacobian of this function.
+ * @param[out] jac       Computed Jacobian.
  * @param[in]  accuracy  Accuracy of the finite differences.
  * @param[in]  eps       Value of the finite difference step.
  */
 void finite_jacobian(
     const Eigen::Ref<const Eigen::VectorXd>& x,
-    const std::function<Eigen::VectorXd(const Eigen::VectorXd&)>& f,
+    const std::function<Eigen::MatrixXd(const Eigen::VectorXd&)>& f,
     Eigen::MatrixXd& jac,
     const AccuracyOrder accuracy = SECOND,
     const double eps = 1.0e-8);
 
 /**
- * @brief Compute the hessian of a function using finite differences.
+ * @brief Compute the Hessian of a function using finite differences.
  *
- * @param[in]  x         Point at which to compute the hessian.
- * @param[in]  f         Compute the hessian of this function.
- * @param[out] hess      Computed hessian.
+ * @param[in]  x         Point at which to compute the Hessian.
+ * @param[in]  f         Compute the Hessian of this function.
+ * @param[out] hess      Computed Hessian.
  * @param[in]  accuracy  Accuracy of the finite differences.
  * @param[in]  eps       Value of the finite difference step.
  */
@@ -87,10 +87,10 @@ bool compare_gradient(
     const std::string& msg = "compare_gradient ");
 
 /**
- * @brief Compare if two jacobians are close enough.
+ * @brief Compare if two Jacobians are close enough.
  *
- * @param[in] x         The first jacobian to compare.
- * @param[in] y         The second jacobian to compare against.
+ * @param[in] x         The first Jacobian to compare.
+ * @param[in] y         The second Jacobian to compare against.
  * @param[in] test_eps  Tolerance of equality.
  * @param[in] msg       Debug message header.
  *
@@ -103,10 +103,10 @@ bool compare_jacobian(
     const std::string& msg = "compare_jacobian ");
 
 /**
- * @brief Compare if two hessians are close enough.
+ * @brief Compare if two Hessians are close enough.
  *
- * @param[in] x         The first hessian to compare.
- * @param[in] y         The second hessian to compare against.
+ * @param[in] x         The first Hessian to compare.
+ * @param[in] y         The second Hessian to compare against.
  * @param[in] test_eps  Tolerance of equality.
  * @param[in] msg       Debug message header.
  *