Fixed mpfr support

Author: rjc8237

src/optimization/CMakeLists.txt

@@ -31,8 +31,7 @@ target_link_libraries(MetricOptimizationLib PUBLIC
   igl::predicates
   spdlog::spdlog
   pybind11::module
-  ${Boost_FILESYSTEM_LIBRARY}
-  ${Boost_SYSTEM_LIBRARY}
+  ${MPFR_LIBRARIES}
 )
 target_compile_definitions(MetricOptimizationLib PUBLIC
   SPDLOG_ACTIVE_LEVEL=SPDLOG_LEVEL_DEBUG

src/optimization/common.hh

@@ -501,6 +501,34 @@ convert_eigen_to_std_vector(const VectorX& vector_eigen,
   }
 }
 
+/// Convert vector of scalars to doubles
+inline Eigen::Matrix<double, Eigen::Dynamic, 1>
+convert_scalar_to_double_vector(const VectorX& vector_scalar)
+{
+  int num_entries = vector_scalar.size();
+  Eigen::Matrix<double, Eigen::Dynamic, 1> vector_double(num_entries);
+  for (int i = 0; i < num_entries; ++i)
+  {
+    vector_double[i] = (double) (vector_scalar[i]);
+  }
+
+  return vector_double;
+}
+
+/// Convert vector of scalars to doubles
+inline std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>>
+convert_scalar_to_double_vector(const std::vector<VectorX>& vector_scalar)
+{
+  int num_entries = vector_scalar.size();
+  std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>> vector_double(num_entries);
+  for (int i = 0; i < num_entries; ++i)
+  {
+    vector_double[i] = convert_scalar_to_double_vector(vector_scalar[i]);
+  }
+
+  return vector_double;
+}
+
 /// *****************
 /// Vector Generation
 /// *****************

src/optimization/energies.cpp

@@ -577,7 +577,7 @@ root_mean_square_error(const VectorX &x, const VectorX &x0)
         error += diff * diff; 
     }
 
-    return std::sqrt(error / num_var);
+    return sqrt(error / num_var);
 }
 
 Scalar
@@ -596,7 +596,7 @@ relative_root_mean_square_error(const VectorX &x, const VectorX &x0)
         denom += (x0[i] * x0[i]);
     }
 
-    return std::sqrt(error / (num_var * denom));
+    return sqrt(error / (num_var * denom));
 }
 
 Scalar
@@ -613,7 +613,7 @@ root_mean_square_relative_error(const VectorX &x, const VectorX &x0)
         error += (rel_err * rel_err);
     }
 
-    return std::sqrt(error / num_var);
+    return sqrt(error / num_var);
 
 }
 

src/optimization/explicit_optimization.cpp

@@ -8,6 +8,8 @@
 #include "nonlinear_optimization.hh"
 #include "projection.hh"
 #include "shear.hh"
+#include <Eigen/SparseLU>
+#include <Eigen/SparseQR>
 
 /// FIXME Do cleaning pass
 
@@ -279,7 +281,9 @@ bool compute_domain_coordinate_energy_with_gradient(
         constraint_codomain_matrix.transpose() * constraint_penner_jacobian.transpose();
 
     // Solve for the component of the gradient corresponding to the implicit metric coordinates
-    Eigen::SparseLU<Eigen::SparseMatrix<Scalar>> solver;
+    //Eigen::SparseLU<MatrixX> solver;
+    //Eigen::SparseLU<Eigen::SparseMatrix<Scalar>, Eigen::COLAMDOrdering<int>> solver;
+    Eigen::SparseQR<Eigen::SparseMatrix<Scalar>, Eigen::COLAMDOrdering<int>> solver;
     solver.compute(constraint_codomain_jacobian);
     VectorX energy_implicit_gradient =
         -constraint_domain_jacobian * solver.solve(energy_codomain_gradient);

src/optimization/implicit_optimization.cpp

@@ -9,6 +9,7 @@
 #include "logging.hh"
 #include "nonlinear_optimization.hh"
 #include "projection.hh"
+#include <Eigen/SparseQR>
 
 /// FIXME Do cleaning pass
 
@@ -347,7 +348,9 @@ void compute_projected_newton_descent_direction(
     // Solve for the optimal descent direction
     igl::Timer timer;
     timer.start();
-    Eigen::SparseLU<Eigen::SparseMatrix<Scalar>> solver;
+    // TODO Make solver a global variable or function so easy to replace
+    //Eigen::SparseLU<Eigen::SparseMatrix<Scalar>> solver;
+    Eigen::SparseQR<Eigen::SparseMatrix<Scalar>, Eigen::COLAMDOrdering<int>> solver;
     solver.compute(hessian_lagrangian);
     VectorX solution = solver.solve(rhs);
     double time = timer.getElapsedTime();

src/optimization/layout.cpp

@@ -140,7 +140,7 @@ bool check_uv(
     }
 
     // Check length consistency
-    double uv_length_error = compute_uv_length_error(F, uv, F_uv);
+    Scalar uv_length_error = compute_uv_length_error(F, uv, F_uv);
     if (!float_equal(uv_length_error, 0.0)) {
         spdlog::warn("Inconsistent uv length error {} across edges", uv_length_error);
     }

src/optimization/refinement.cpp

@@ -191,24 +191,24 @@ RefinementMesh::get_face_iterator(Index face_index) const
   return FaceIterator(*this, face_index);
 }
 
-VectorX
+Eigen::VectorXd
 RefinementMesh::get_vertex(
 	RefinementMesh::Index vertex_index
 ) const {
-	return m_V.row(vertex_index);
+	return m_V.row(vertex_index).transpose();
 }
 
-VectorX
+Eigen::VectorXd
 RefinementMesh::get_uv_vertex(
 	RefinementMesh::Index vertex_index
 ) const {
-	return m_uv.row(vertex_index);
+	return m_uv.row(vertex_index).transpose();
 }
 
 void
 RefinementMesh::get_face_vertices(
 	RefinementMesh::Index face_index,
-	std::vector<VectorX>& vertices
+	std::vector<Eigen::VectorXd>& vertices
 ) const {
 	vertices.clear();
 	for (auto iter = get_face_iterator(face_index); !iter.done(); ++iter)
@@ -222,7 +222,7 @@ RefinementMesh::get_face_vertices(
 void
 RefinementMesh::get_face_uv_vertices(
 	RefinementMesh::Index face_index,
-	std::vector<VectorX>& uv_vertices
+	std::vector<Eigen::VectorXd>& uv_vertices
 ) const {
 	uv_vertices.clear();
 	for (auto iter = get_face_iterator(face_index); !iter.done(); ++iter)
@@ -744,7 +744,7 @@ RefinementMesh::refine_mesh()
 		// Skip boundary loops
 		if (is_bnd_loop[fijk]) continue;
 
-		std::array<VectorX, 3> triangle;
+		std::array<Eigen::VectorXd, 3> triangle;
 		Index hij = h[fijk];
 		for (int l = 0; l < 3; ++l)
 		{
@@ -1115,7 +1115,7 @@ RefinementMesh::is_self_overlapping_face(
 	if (is_bnd_loop[face_index]) return true;
 
 	// Get vertices of the face
-	std::vector<VectorX> uv_vertices, vertices;
+	std::vector<Eigen::VectorXd> uv_vertices, vertices;
 	get_face_uv_vertices(face_index, uv_vertices);
 	get_face_vertices(face_index, vertices);
 
@@ -1151,8 +1151,8 @@ RefinementMesh::triangulate_face(
 	// Get uv vertices of the face
 	std::vector<int> vertex_indices;
 	std::vector<int> uv_vertex_indices;
-	std::vector<VectorX> vertices;
-	std::vector<VectorX> uv_vertices;
+	std::vector<Eigen::VectorXd> vertices;
+	std::vector<Eigen::VectorXd> uv_vertices;
 	for (auto iter = get_face_iterator(face_index); !iter.done(); ++iter)
 	{
 		Index hij = *iter;

src/optimization/refinement.hh

@@ -172,7 +172,7 @@ public:
 	///
 	/// @param[in] vertex_index: index of the vertex in the mesh
 	/// @return vertex position
-	VectorX
+	Eigen::VectorXd
 	get_vertex(
 		Index vertex_index
 	) const;
@@ -181,7 +181,7 @@ public:
 	///
 	/// @param[in] vertex_index: index of the vertex in the mesh
 	/// @return parametric domain coordinates of the vertex
-	VectorX
+	Eigen::VectorXd
 	get_uv_vertex(
 		Index vertex_index
 	) const;
@@ -193,7 +193,7 @@ public:
 	void
 	get_face_vertices(
 		Index face_index,
-		std::vector<VectorX>& vertices
+		std::vector<Eigen::VectorXd>& vertices
 	) const;
 
 	/// Get the parametric domain coordinates of the vertices of a face
@@ -203,7 +203,7 @@ public:
 	void
 	get_face_uv_vertices(
 		Index face_index,
-		std::vector<VectorX>& uv_vertices
+		std::vector<Eigen::VectorXd>& uv_vertices
 	) const;
 
 	/// Clear the data of the mesh

src/optimization/translation.cpp

@@ -5,6 +5,7 @@
 #include "reparametrization.hh"
 #include "shear.hh"
 #include <Eigen/SparseLU>
+#include <Eigen/SparseQR>
 
 /// FIXME Do cleaning pass
 
@@ -132,7 +133,8 @@ compute_as_symmetric_as_possible_translations(
     lagrangian_matrix.cols(),
     right_hand_side.size()
   );
-  Eigen::SparseLU<Eigen::SparseMatrix<Scalar>> solver;
+  //Eigen::SparseLU<Eigen::SparseMatrix<Scalar>, Eigen::COLAMDOrdering<int>> solver;
+  Eigen::SparseQR<Eigen::SparseMatrix<Scalar>, Eigen::COLAMDOrdering<int>> solver;
   solver.compute(lagrangian_matrix);
   VectorX lagrangian_solution = solver.solve(-right_hand_side);
 

src/optimization/triangulation.cpp

@@ -23,7 +23,7 @@ namespace CurvatureMetric {
 
 Scalar
 compute_face_area(
-	const std::array<VectorX, 3>& vertices
+	const std::array<Eigen::VectorXd, 3>& vertices
 ) {
 	// Get edge lengths for triangle
 	Scalar li = (vertices[1] - vertices[0]).norm();
@@ -36,7 +36,7 @@ compute_face_area(
 
 bool
 is_inverted_triangle(
-	const std::array<VectorX, 3>& vertices
+	const std::array<Eigen::VectorXd, 3>& vertices
 ) {
 	// Build matrix of triangle homogenous coordinates
 	Eigen::Matrix<Scalar, 3, 3> tri_homogenous_coords;
@@ -45,15 +45,15 @@ is_inverted_triangle(
 	tri_homogenous_coords.col(2) << vertices[2][0], vertices[2][1], 1.0;
 
 	// Triangle is flipped iff the determinant is negative
-	double det = tri_homogenous_coords.determinant();
+	Scalar det = tri_homogenous_coords.determinant();
 	return (det < 0.0);
 }
 
 
 bool
 is_self_overlapping_polygon(
-	const std::vector<VectorX>& uv_vertices,
-	const std::vector<VectorX>& vertices,
+	const std::vector<Eigen::VectorXd>& uv_vertices,
+	const std::vector<Eigen::VectorXd>& vertices,
 	std::vector<std::vector<bool>>& is_self_overlapping_subpolygon,
 	std::vector<std::vector<int>>& splitting_vertices,
 	std::vector<std::vector<Scalar>>& min_face_areas 
@@ -108,8 +108,8 @@ is_self_overlapping_polygon(
 			for (int k = (i + 1) % face_size; k != j; k = (k + 1) % face_size)
 			{
 				// Check if triangle T_ikj is positively oriented
-				std::array<VectorX, 3> uv_triangle = { uv_vertices[i], uv_vertices[k], uv_vertices[j] };
-				std::array<VectorX, 3> triangle = { vertices[i], vertices[k], vertices[j] };
+				std::array<Eigen::VectorXd, 3> uv_triangle = { uv_vertices[i], uv_vertices[k], uv_vertices[j] };
+				std::array<Eigen::VectorXd, 3> triangle = { vertices[i], vertices[k], vertices[j] };
 				if (is_inverted_triangle(uv_triangle)) continue;
 
 				// Check if the two subpolygons (i, k) and (k, j) are self overlapping
@@ -266,8 +266,8 @@ triangulate_self_overlapping_polygon(
 // Helper function to view the triangulated face
 void
 view_triangulation(
-	const std::vector<VectorX>& uv_vertices,
-	const std::vector<VectorX>& vertices,
+	const std::vector<Eigen::VectorXd>& uv_vertices,
+	const std::vector<Eigen::VectorXd>& vertices,
 	const std::vector<std::vector<bool>>& is_self_overlapping_subpolygon,
 	const std::vector<std::vector<int>>& splitting_vertices,
 	const std::vector<std::vector<Scalar>>& min_face_areas,