Merge pull request #2085 from SCIInstitute/cleanup_constraints

Cleanup Constraints, added comments, docs, and a roadmap

Author: akenmorris

.github/workflows/build-windows.yml

@@ -115,7 +115,7 @@ jobs:
       # Note the current convention is to use the -S and -B options here to specify source 
       # and build directories, but this is only available with CMake 3.13 and higher.  
       # The CMake binaries on the Github Actions machines are (as of this writing) 3.12
-      run: conda activate shapeworks && cmake $GITHUB_WORKSPACE -DCMAKE_CXX_FLAGS="-FS" -DCMAKE_C_FLAGS="-FS" -DCMAKE_CXX_FLAGS_RELEASE="-FS /Zm250 /Zi /GL /MD /O2 /Ob3 /DNDEBUG /EHsc" -DCMAKE_C_FLAGS_RELEASE="-FS /Zi /GL /MD /O2 /Ob3 /DNDEBUG /EHsc" -DCMAKE_SHARED_LINKER_FLAGS_RELEASE="-LTCG /DEBUG" -DCMAKE_EXE_LINKER_FLAGS_RELEASE="-LTCG /DEBUG" -DCMAKE_BUILD_TYPE=$BUILD_TYPE -DITK_DIR="C:\deps\lib\cmake\ITK-5.2" -DVTK_DIR="C:\deps\lib\cmake\vtk-9.1" -DXLNT_DIR="C:\deps" -DLIBIGL_DIR="C:\deps" -DJKQTCommonSharedLib_DIR="C:/deps/lib/cmake/JKQTCommonSharedLib" -DJKQTMathTextSharedLib_DIR="C:/deps/lib/cmake/JKQTMathTextSharedLib" -DJKQTPlotterSharedLib_DIR="C:/deps/lib/cmake/JKQTPlotterSharedLib" -DOpenVDB_DIR="C:\deps\lib\cmake\OpenVDB" -DGEOMETRYCENTRAL_DIR="C:\deps" -DACVD_DIR="C:\deps" -DBuild_Studio=ON -DGA_MEASUREMENT_ID=$GA_MEASUREMENT_ID -DGA_API_SECRET=$GA_API_SECRET
+      run: conda activate shapeworks && cmake $GITHUB_WORKSPACE -DCMAKE_CXX_FLAGS="-FS" -DCMAKE_C_FLAGS="-FS" -DCMAKE_CXX_FLAGS_RELEASE="-FS /Zm500 /Zi /GL /MD /O2 /Ob3 /DNDEBUG /EHsc" -DCMAKE_C_FLAGS_RELEASE="-FS /Zi /GL /MD /O2 /Ob3 /DNDEBUG /EHsc" -DCMAKE_SHARED_LINKER_FLAGS_RELEASE="-LTCG /DEBUG" -DCMAKE_EXE_LINKER_FLAGS_RELEASE="-LTCG /DEBUG" -DCMAKE_BUILD_TYPE=$BUILD_TYPE -DITK_DIR="C:\deps\lib\cmake\ITK-5.2" -DVTK_DIR="C:\deps\lib\cmake\vtk-9.1" -DXLNT_DIR="C:\deps" -DLIBIGL_DIR="C:\deps" -DJKQTCommonSharedLib_DIR="C:/deps/lib/cmake/JKQTCommonSharedLib" -DJKQTMathTextSharedLib_DIR="C:/deps/lib/cmake/JKQTMathTextSharedLib" -DJKQTPlotterSharedLib_DIR="C:/deps/lib/cmake/JKQTPlotterSharedLib" -DOpenVDB_DIR="C:\deps\lib\cmake\OpenVDB" -DGEOMETRYCENTRAL_DIR="C:\deps" -DACVD_DIR="C:\deps" -DBuild_Studio=ON -DGA_MEASUREMENT_ID=$GA_MEASUREMENT_ID -DGA_API_SECRET=$GA_API_SECRET
 
     - name: Build
       working-directory: "C:/build"

Libs/Optimize/Constraints/Constraint.cpp

@@ -2,58 +2,6 @@
 #include <iostream>
 namespace shapeworks {
 
-void Constraint::updateZ(const Eigen::Vector3d &pt, double C) {
-  /*double a = 2*C;
-  double c = 2*mu + 2*C*ConstraintEval(pt);
-  if(c >= 0){
-      z = 0;
-      //std::cout << "if    z: " << z << std::endl;
-  }
-  else{
-      double z1 = std::sqrt(-c/a);
-      double z2 = -z1;
-      //std::cout << "else    z: " << z1 << std::endl;
-
-      double pteval = ConstraintEval(pt);
-
-      double res1 = mu*(pteval + z1*z1) + C/2*std::abs(pteval + z1*z1)*std::abs(pteval + z1*z1);
-      double res2 = mu*(pteval + z2*z2) + C/2*std::abs(pteval + z2*z2)*std::abs(pteval + z2*z2);
-
-      if(res1 < res2){
-          z = z1;
-      }
-      else{
-          z = z2;
-      }
-  }*/
-
-  // Augmented lagrangian inequality equation: f(x) = mu*(g(x)+z^2) + C/2|g(x)+z^2|^2
-  // f'(x) = mu*g'(x) + C*y' where by substitution
-  // y = √(u^2) where by substitution
-  // u = g(x) + z^2
-  // u' = 2*z
-  //
-  // Then we compute y'
-  // y' = (dy / du) (du / dx)
-  //    = (1/2)*(2 * u) / √(u^2) (du / dx)
-  //    = u * u' / | u |
-  //    = sgn(u) * u'
-  //
-  // So we substitute
-  // f'(x) = 2*mu*z + 2*C*sgn(g(x)+z^2)*z
-
-  // z iterative update as explained above
-  /*double update = 1000;
-  size_t count = 1000;
-  while(update < 0.1|| count > 0){
-      std::cout << "pt: " << pt.transpose() << " count: " << count << " z: " << z << std::endl;
-      update = 2*mu*z + 2*C*sgn(ConstraintEval(pt)+z*z)*z;
-      z = z + update;
-      count--;
-  }*/
-  // std::cout << "z: " << z << std::endl;
-}
-
 void Constraint::updateMu(const Eigen::Vector3d &pt, double C, size_t index) {
   double eval = -constraintEval(pt);
   double maxterm = mus_[index] + C * eval;

Libs/Optimize/Constraints/Constraint.h

@@ -10,36 +10,48 @@
 
 namespace shapeworks {
 
+/**
+ * \class Constraint
+ * \ingroup Group-Constraints
+ *
+ *  This class is the general constraint class. Each instance represents a single constraint, either cutting-plane, sphere or free-form. They all inherit from this class. This class containts all the infrastructure to handle gradients and evaluations, which is shared among all constraint types.
+ *  NOTE: Not actually using the augmented lagrangian. We are using quadratic penalty and not lagrangian because it works better.
+ *
+ */
+
 class Constraint {
  public:
+  /// Returns if pt in vnl_vector format is violated by the constraint
   bool isViolated(const vnl_vector<double> &pt) const { return isViolated(Eigen::Vector3d(pt[0], pt[1], pt[2])); }
+  /// Returns if pt in Eigen format is violated by the constraint
   virtual bool isViolated(const Eigen::Vector3d &pt) const = 0;
+  /// Prints the constraint neatly
   virtual void print() const = 0;
 
   // For augmented lagrangian
-  void setZ(double inz) { z_ = inz; }
-  double getZ() { return z_; }
+  /// Initializes mu
   void setMus(std::vector<double> inmu) { mus_ = inmu; }
+  /// Gets mu
   std::vector<double> getMus() { return mus_; }
-  void setLambda(double inLambda) { lambda_ = inLambda; }
-  double getLambda() { return lambda_; }
 
+  /// Returns the gradient of the constraint
   virtual Eigen::Vector3d constraintGradient(const Eigen::Vector3d &pt) const = 0;
+  /// Returns the evaluation on the constraint, i.e. the signed distance to the constraint boundary
   virtual double constraintEval(const Eigen::Vector3d &pt) const = 0;
 
-  void updateZ(const Eigen::Vector3d &pt, double C);
-
+  /// Updates the value of mu according to the augmented lagrangian update
   void updateMu(const Eigen::Vector3d &pt, double C, size_t index);
 
+  /// Computes the lagrangian gradient based on lagrangian inequality equations. NOTE: Not actually lagrangian. We are using quadratic penalty and not lagrangian because it works better.
   Eigen::Vector3d lagragianGradient(const Eigen::Vector3d &pt, double C, size_t index) const;
 
  protected:
+  /// Returns the sign of the double
   int sgn(double val) { return (double(0) < val) - (val < double(0)); }
 
   // For augmented lagrangian
+  /// Mu is the lagrangian momentum term
   std::vector<double> mus_;
-  double z_;
-  double lambda_;
 };
 
 }  // namespace shapeworks

Libs/Optimize/Constraints/ConstraintType.h

@@ -4,7 +4,6 @@
 
 namespace shapeworks {
   enum class ConstraintType : char {
-    Sphere = 'S',
     CuttingPlane = 'C',
     FreeForm = 'F'
   };