New project PSD option (#95)

* New project PSD option
* Update scalable_ccd

---------

Co-authored-by: Zizhou Huang <zizhou@nyu.edu>
Co-authored-by: Zachary Ferguson <zy.fergus@gmail.com>

cmake/recipes/scalable_ccd.cmake

@@ -9,4 +9,4 @@ message(STATUS "Third-party: creating target 'scalable_ccd::scalable_ccd'")
 set(SCALABLE_CCD_WITH_CUDA ${IPC_TOOLKIT_WITH_CUDA} CACHE BOOL "Enable CUDA CCD" FORCE)
 
 include(CPM)
-CPMAddPackage("gh:continuous-collision-detection/scalable-ccd#60692d7cbc03390a3a198792a3ff2b869c5da05b")
+CPMAddPackage("gh:continuous-collision-detection/scalable-ccd#318352e6aa3c1008a67c80547c6fa21bf18fb1a5")

python/src/bindings.cpp

@@ -10,6 +10,9 @@ PYBIND11_MODULE(ipctk, m)
 
     m.doc() = "IPC Toolkit";
 
+    // define these early because they are used in other definitions
+    define_eigen_ext(m);
+
     // barrier
     define_barrier(m);
     define_adaptive_stiffness(m);
@@ -91,7 +94,7 @@ PYBIND11_MODULE(ipctk, m)
 
     // utils
     define_area_gradient(m);
-    define_eigen_ext(m);
+    // define_eigen_ext(m);
     define_interval(m);
     define_intersection(m);
     define_logger(m);

python/src/potentials/barrier_potential.cpp

@@ -57,7 +57,8 @@ void define_barrier_potential(py::module_& m)
             "hessian",
             py::overload_cast<
                 const Collisions&, const CollisionMesh&, const Eigen::MatrixXd&,
-                const bool>(&BarrierPotential::Potential::hessian, py::const_),
+                const PSDProjectionMethod>(
+                &BarrierPotential::Potential::hessian, py::const_),
             R"ipc_Qu8mg5v7(
             Compute the hessian of the barrier potential.
 
@@ -71,7 +72,7 @@ void define_barrier_potential(py::module_& m)
                 The hessian of all barrier potentials (not scaled by the barrier stiffness). This will have a size of |vertices|x|vertices|.
             )ipc_Qu8mg5v7",
             py::arg("collisions"), py::arg("mesh"), py::arg("vertices"),
-            py::arg("project_hessian_to_psd") = false)
+            py::arg("project_hessian_to_psd") = PSDProjectionMethod::NONE)
         .def(
             "shape_derivative",
             py::overload_cast<
@@ -125,7 +126,8 @@ void define_barrier_potential(py::module_& m)
         .def(
             "hessian",
             py::overload_cast<
-                const Collision&, const VectorMax12d&, const bool>(
+                const Collision&, const VectorMax12d&,
+                const PSDProjectionMethod>(
                 &BarrierPotential::hessian, py::const_),
             R"ipc_Qu8mg5v7(
             Compute the hessian of the potential for a single collision.
@@ -138,7 +140,7 @@ void define_barrier_potential(py::module_& m)
                 The hessian of the potential.
             )ipc_Qu8mg5v7",
             py::arg("collision"), py::arg("x"),
-            py::arg("project_hessian_to_psd") = false)
+            py::arg("project_hessian_to_psd") = PSDProjectionMethod::NONE)
         .def(
             "shape_derivative",
             [](const DistanceBasedPotential& self, const Collision& collision,

python/src/potentials/friction_potential.cpp

@@ -72,7 +72,7 @@ void define_friction_potential(py::module_& m)
             "hessian",
             py::overload_cast<
                 const FrictionCollisions&, const CollisionMesh&,
-                const Eigen::MatrixXd&, const bool>(
+                const Eigen::MatrixXd&, const PSDProjectionMethod>(
                 &FrictionPotential::Potential::hessian, py::const_),
             R"ipc_Qu8mg5v7(
             Compute the hessian of the friction dissipative potential.
@@ -87,7 +87,7 @@ void define_friction_potential(py::module_& m)
                 The hessian of all friction dissipative potentials. This will have a size of |velocities|×|velocities|.
             )ipc_Qu8mg5v7",
             py::arg("collisions"), py::arg("mesh"), py::arg("vertices"),
-            py::arg("project_hessian_to_psd") = false)
+            py::arg("project_hessian_to_psd") = PSDProjectionMethod::NONE)
         .def(
             "force",
             py::overload_cast<
@@ -179,7 +179,8 @@ void define_friction_potential(py::module_& m)
         .def(
             "hessian",
             py::overload_cast<
-                const FrictionCollision&, const VectorMax12d&, const bool>(
+                const FrictionCollision&, const VectorMax12d&,
+                const PSDProjectionMethod>(
                 &FrictionPotential::hessian, py::const_),
             R"ipc_Qu8mg5v7(
             Compute the hessian of the potential for a single collision.
@@ -192,7 +193,7 @@ void define_friction_potential(py::module_& m)
                 The hessian of the potential.
             )ipc_Qu8mg5v7",
             py::arg("collision"), py::arg("x"),
-            py::arg("project_hessian_to_psd") = false)
+            py::arg("project_hessian_to_psd") = PSDProjectionMethod::NONE)
         .def(
             "force",
             py::overload_cast<

python/src/utils/eigen_ext.cpp

@@ -23,6 +23,15 @@ void define_eigen_ext(py::module_& m)
         )ipc_Qu8mg5v7",
         py::arg("A"), py::arg("eps") = 1e-8);
 
+    py::enum_<PSDProjectionMethod>(
+        m, "PSDProjectionMethod",
+        "Enumeration of implemented PSD projection methods.")
+        .value("NONE", PSDProjectionMethod::NONE, "No PSD projection")
+        .value(
+            "CLAMP", PSDProjectionMethod::CLAMP, "Clamp negative eigenvalues")
+        .value("ABS", PSDProjectionMethod::ABS, "Flip negative eigenvalues")
+        .export_values();
+
     m.def(
         "project_to_psd",
         &project_to_psd<
@@ -33,9 +42,10 @@ void define_eigen_ext(py::module_& m)
 
         Parameters:
             A: Symmetric matrix to project
+            method: PSD projection method
 
         Returns:
             Projected matrix
         )ipc_Qu8mg5v7",
-        py::arg("A"));
+        py::arg("A"), py::arg("method"));
 }

python/tests/find_ipctk.py

@@ -3,6 +3,17 @@
 except ImportError:
     import sys
     import pathlib
-    sys.path.append(
-        str(pathlib.Path(__file__).parents[2] / "build" / "release" / "python"))
+    repo_root = pathlib.Path(__file__).parents[2]
+    possible_paths = [
+        pathlib.Path("python").resolve(),
+        repo_root / "build" / "python",
+        repo_root / "build" / "release" / "python",
+        repo_root / "build" / "debug" / "python",
+    ]
+    for path in possible_paths:
+        if path.exists():
+            sys.path.append(str(path))
+            break
+    else:
+        raise ImportError("Could not find the ipctk module")
     import ipctk  # Try again

src/ipc/potentials/distance_based_potential.cpp

@@ -89,7 +89,7 @@ VectorMax12d DistanceBasedPotential::gradient(
 MatrixMax12d DistanceBasedPotential::hessian(
     const Collision& collision,
     const VectorMax12d& positions,
-    const bool project_hessian_to_psd) const
+    const PSDProjectionMethod project_hessian_to_psd) const
 {
     // d(x)
     const double d = collision.compute_distance(positions);
@@ -135,7 +135,7 @@ MatrixMax12d DistanceBasedPotential::hessian(
     }
 
     // Need to project entire hessian because w can be negative
-    return project_hessian_to_psd ? project_to_psd(hess) : hess;
+    return project_to_psd(hess, project_hessian_to_psd);
 }
 
 void DistanceBasedPotential::shape_derivative(
@@ -181,8 +181,7 @@ void DistanceBasedPotential::shape_derivative(
     MatrixMax12d local_hess;
     if (!collision.is_mollified()) {
         // w ∇ₓ∇ᵤf = w ∇ᵤ²f
-        local_hess =
-            hessian(collision, positions, /*project_hessian_to_psd=*/false);
+        local_hess = hessian(collision, positions, PSDProjectionMethod::NONE);
     } else {
         // d(x̄+u)
         const double d = collision.compute_distance(positions);

src/ipc/potentials/distance_based_potential.hpp

@@ -56,7 +56,8 @@ class DistanceBasedPotential : public Potential<Collisions> {
     MatrixMax12d hessian(
         const Collision& collision,
         const VectorMax12d& positions,
-        const bool project_hessian_to_psd = false) const override;
+        const PSDProjectionMethod project_hessian_to_psd =
+            PSDProjectionMethod::NONE) const override;
 
     /// @brief Compute the shape derivative of the potential for a single collision.
     /// @param[in] collision The collision.

src/ipc/potentials/friction_potential.cpp

@@ -183,7 +183,7 @@ VectorMax12d FrictionPotential::gradient(
 MatrixMax12d FrictionPotential::hessian(
     const FrictionCollision& collision,
     const VectorMax12d& velocities,
-    const bool project_hessian_to_psd) const
+    const PSDProjectionMethod project_hessian_to_psd) const
 {
     // ∇ₓ μ N(xᵗ) f₁(‖u‖)/‖u‖ T(xᵗ) u (where u = T(xᵗ)ᵀ v)
     //  = μ N T [(f₁'(‖u‖)‖u‖ − f₁(‖u‖))/‖u‖³ uuᵀ + f₁(‖u‖)/‖u‖ I] Tᵀ
@@ -215,7 +215,7 @@ MatrixMax12d FrictionPotential::hessian(
         //            = μ N T [-f₁(‖u‖)/‖u‖ uuᵀ/‖u‖² + f₁(‖u‖)/‖u‖ I] Tᵀ
         //            = μ N T [f₁(‖u‖)/‖u‖ (I - uuᵀ/‖u‖²)] Tᵀ
         //  ⟹ no PSD projection needed because f₁(‖u‖)/‖u‖ ≥ 0
-        if (project_hessian_to_psd && scale <= 0) {
+        if (project_hessian_to_psd != PSDProjectionMethod::NONE && scale <= 0) {
             hess.setZero(collision.ndof(), collision.ndof()); // -PSD = NSD ⟹ 0
         } else if (collision.dim() == 2) {
             // I - uuᵀ/‖u‖² = 1 - u²/u² = 0 ⟹ ∇²D(v) = 0
@@ -232,7 +232,7 @@ MatrixMax12d FrictionPotential::hessian(
         // ∇²D = μ N T [(f₁'(‖u‖)‖u‖ − f₁(‖u‖))/‖u‖³ uuᵀ + f₁(‖u‖)/‖u‖ I] Tᵀ
         // lim_{‖u‖→0} ∇²D = μ N T [f₁(‖u‖)/‖u‖ I] Tᵀ
         // no PSD projection needed because μ N f₁(‖ū‖)/‖ū‖ ≥ 0
-        if (project_hessian_to_psd && scale <= 0) {
+        if (project_hessian_to_psd != PSDProjectionMethod::NONE && scale <= 0) {
             hess.setZero(collision.ndof(), collision.ndof()); // -PSD = NSD ⟹ 0
         } else {
             hess = scale * f1_over_norm_u * T * T.transpose();
@@ -245,9 +245,7 @@ MatrixMax12d FrictionPotential::hessian(
         MatrixMax2d inner_hess = f2 * u * u.transpose();
         inner_hess.diagonal().array() += f1_over_norm_u;
         inner_hess *= scale; // NOTE: negative scaling will be projected out
-        if (project_hessian_to_psd) {
-            inner_hess = project_to_psd(inner_hess);
-        }
+        inner_hess = project_to_psd(inner_hess, project_hessian_to_psd);
 
         hess = T * inner_hess * T.transpose();
     }

src/ipc/potentials/friction_potential.hpp

@@ -109,7 +109,8 @@ class FrictionPotential : public Potential<FrictionCollisions> {
     MatrixMax12d hessian(
         const FrictionCollision& collision,
         const VectorMax12d& velocities,
-        const bool project_hessian_to_psd = false) const override;
+        const PSDProjectionMethod project_hessian_to_psd =
+            PSDProjectionMethod::NONE) const override;
 
     /// @brief Compute the friction force.
     /// @param collision The collision

src/ipc/potentials/potential.hpp

@@ -47,7 +47,8 @@ template <class TCollisions> class Potential {
         const TCollisions& collisions,
         const CollisionMesh& mesh,
         const Eigen::MatrixXd& X,
-        const bool project_hessian_to_psd = false) const;
+        const PSDProjectionMethod project_hessian_to_psd =
+            PSDProjectionMethod::NONE) const;
 
     // -- Single collision methods ---------------------------------------------
 
@@ -72,7 +73,8 @@ template <class TCollisions> class Potential {
     virtual MatrixMax12d hessian(
         const TCollision& collision,
         const VectorMax12d& x,
-        const bool project_hessian_to_psd = false) const = 0;
+        const PSDProjectionMethod project_hessian_to_psd =
+            PSDProjectionMethod::NONE) const = 0;
 };
 
 } // namespace ipc

src/ipc/potentials/potential.tpp

@@ -84,7 +84,7 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
     const TCollisions& collisions,
     const CollisionMesh& mesh,
     const Eigen::MatrixXd& X,
-    const bool project_hessian_to_psd) const
+    const PSDProjectionMethod project_hessian_to_psd) const
 {
     assert(X.rows() == mesh.num_vertices());
 

src/ipc/utils/eigen_ext.hpp

@@ -138,6 +138,8 @@ project_to_pd(
     const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& A,
     double eps = 1e-8);
 
+enum class PSDProjectionMethod { NONE, CLAMP, ABS };
+
 /// @brief Matrix projection onto positive semi-definite cone
 /// @param A Symmetric matrix to project
 /// @return Projected matrix
@@ -150,8 +152,8 @@ template <
     int _MaxCols>
 Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
 project_to_psd(
-    const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
-        A);
+    const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& A,
+    const PSDProjectionMethod method = PSDProjectionMethod::CLAMP);
 
 inline Eigen::Vector3d to_3D(const VectorMax3d& v)
 {

src/ipc/utils/eigen_ext.tpp

@@ -63,10 +63,14 @@ template <
     int _MaxCols>
 Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
 project_to_psd(
-    const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& A)
+    const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& A,
+    const PSDProjectionMethod method)
 {
     assert(A.isApprox(A.transpose()) && "A must be symmetric");
 
+    if (method == PSDProjectionMethod::NONE)
+        return A;
+
     // https://math.stackexchange.com/q/2776803
     Eigen::SelfAdjointEigenSolver<
         Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>>
@@ -86,11 +90,21 @@ project_to_psd(
     // Save a little time and only project the negative values
     for (int i = 0; i < A.rows(); i++) {
         if (D.diagonal()[i] < 0.0) {
-            D.diagonal()[i] = 0.0;
+            switch (method) {
+            case PSDProjectionMethod::CLAMP:
+                D.diagonal()[i] = 0.0;
+                break;
+            case PSDProjectionMethod::ABS:
+                D.diagonal()[i] = std::abs(D.diagonal()[i]);
+                break;
+            default:
+                throw std::runtime_error("Invalid type of PSD projection!");
+            }
         } else {
             break;
         }
     }
+
     return eigensolver.eigenvectors() * D
         * eigensolver.eigenvectors().transpose();
 }

tests/src/tests/potential/test_barrier_potential.cpp

@@ -507,7 +507,8 @@ TEST_CASE(
     };
     BENCHMARK("Compute barrier potential hessian with PSD projection")
     {
-        return barrier_potential.hessian(collisions, mesh, vertices, true);
+        return barrier_potential.hessian(
+            collisions, mesh, vertices, PSDProjectionMethod::CLAMP);
     };
     BENCHMARK("Compute compute_minimum_distance")
     {

tests/src/tests/utils/test_utils.cpp

@@ -29,17 +29,21 @@ TEST_CASE("Project to PSD", "[utils][project_to_psd]")
     Eigen::MatrixXd A, A_psd;
 
     A.setIdentity(3, 3);
-    A_psd = ipc::project_to_psd(A);
+    A_psd = ipc::project_to_psd(A, ipc::PSDProjectionMethod::CLAMP);
+    CHECK(A_psd.isApprox(A));
+    A_psd = ipc::project_to_psd(A, ipc::PSDProjectionMethod::ABS);
     CHECK(A_psd.isApprox(A));
 
     A *= -1;
-    A_psd = ipc::project_to_psd(A);
+    A_psd = ipc::project_to_psd(A, ipc::PSDProjectionMethod::CLAMP);
     CHECK(A_psd.isZero());
 
     A.resize(2, 2);
     A.row(0) << 2, 1;
     A.row(1) << 1, 2;
-    A_psd = ipc::project_to_psd(A);
+    A_psd = ipc::project_to_psd(A, ipc::PSDProjectionMethod::CLAMP);
+    CHECK(A_psd.isApprox(A));
+    A_psd = ipc::project_to_psd(A, ipc::PSDProjectionMethod::ABS);
     CHECK(A_psd.isApprox(A));
 }