Added rmse computations and skip mesh options to the pipeline

Author: rjc8237

CMakeLists.txt

@@ -152,25 +152,22 @@ if(NOT USE_MULTIPRECISION)
     include(Catch2)
 
     # Build testing libraries
-    add_library(CurvatureMetricTestsLib
-      src/tests/test_area.cpp
-      src/tests/test_energies.cpp
-      src/tests/test_optimize.cpp
-    )
-    target_link_libraries(CurvatureMetricTestsLib PUBLIC
-      MetricOptimizationLib
-      Catch2::Catch2WithMain
-    )
-    set(TEST_DATA_ROOT "${PROJECT_SOURCE_DIR}/src/tests/regression/")
-    target_compile_definitions(CurvatureMetricTestsLib PUBLIC TEST_DATA_DIR=\"${TEST_DATA_ROOT}\")
+    #add_library(CurvatureMetricTestsLib
+    #)
+    #target_link_libraries(CurvatureMetricTestsLib PUBLIC
+    #  MetricOptimizationLib
+    #  Catch2::Catch2WithMain
+    #)
 
     # Build testing executable
     add_executable(CurvatureMetricTests
       src/tests/tests.cpp
     )
     target_link_libraries(CurvatureMetricTests PRIVATE
-      CurvatureMetricTestsLib
+      MetricOptimizationLib
       Catch2::Catch2WithMain
     )
+    set(TEST_DATA_ROOT "${PROJECT_SOURCE_DIR}/src/tests/regression/")
+    target_compile_definitions(CurvatureMetricTests PUBLIC TEST_DATA_DIR=\"${TEST_DATA_ROOT}\")
   endif()
 endif()

cmake/suitesparse.cmake

@@ -1,8 +1,12 @@
+if(TARGET SuiteSparse::SuiteSparseConfig)
+  return()
+endif()
+
 include(FetchContent)
 FetchContent_Declare(
     suitesparse
 		SYSTEM
     GIT_REPOSITORY https://github.com/DrTimothyAldenDavis/SuiteSparse.git
 		GIT_TAG stable
 )
-FetchContent_MakeAvailable(suitesparse)
+FetchContent_MakeAvailable(suitesparse)

scripts/__pycache__/script_util.cpython-39.pyc

@@ -23,6 +23,8 @@ def add_mesh_arguments(parser):
     # Input and output directories
     parser.add_argument("-f", "--fname",         help="filenames of the obj file", 
                                                      nargs='+')
+    parser.add_argument("--skip_fname",         help="filenames of the obj files to skip", 
+                                                     nargs='+')
     parser.add_argument("-i", "--input_dir",     help="input folder that stores obj files and Th_hat")
     parser.add_argument("--lambdas_dir",         help="directory for initial lambdas")
     parser.add_argument("--free_bd_angles",      help="let boundary angles be free",
@@ -348,7 +350,7 @@ def count_running_processes(p_list):
     return count
 
 def run_many(method, args):
-    args_list = [(args, fname) for fname in args['fname']]
+    args_list = [(args, fname) for fname in args['fname'] if fname not in args['skip_fname']]
     with multiprocessing.Pool(processes=args['num_processes']) as pool:
         pool.starmap(method, args_list, chunksize=1)
 

scripts/script_util.py

@@ -563,6 +563,60 @@ namespace CurvatureMetric
     return energy;
   }
 
+Scalar
+root_mean_square_error(const VectorX &x, const VectorX &x0)
+{
+    assert(x.size() == x0.size());
+    int num_var = x.size();
+    if (num_var == 0) return 0.0;
+
+    Scalar error = 0.0;
+    for (int i = 0; i < num_var; ++i)
+    {
+        Scalar diff = x[i] - x0[i];
+        error += diff * diff; 
+    }
+    
+    return std::sqrt(error / num_var);
+}
+
+Scalar
+relative_root_mean_square_error(const VectorX &x, const VectorX &x0)
+{
+    assert(x.size() == x0.size());
+    int num_var = x.size();
+    if (num_var == 0) return 0.0;
+
+    Scalar error = 0.0;
+    Scalar denom = 0.0;
+    for (int i = 0; i < num_var; ++i)
+    {
+        Scalar diff = x[i] - x0[i];
+        error += (diff * diff);
+        denom += (x0[i] * x0[i]);
+    }
+    
+    return std::sqrt(error / (num_var * denom));
+}
+
+Scalar
+root_mean_square_relative_error(const VectorX &x, const VectorX &x0)
+{
+    assert(x.size() == x0.size());
+    int num_var = x.size();
+    if (num_var == 0) return 0.0;
+
+    Scalar error = 0.0;
+    for (int i = 0; i < num_var; ++i)
+    {
+        Scalar rel_err = (!float_equal(x0[i], 0.0)) ? (x[i] - x0[i]) / x0[i] : 1e10;
+        error += (rel_err * rel_err);
+    }
+    
+    return std::sqrt(error / num_var);
+    
+}
+
 // FIXME Rename these variables
 // FIXME Ensure all pybind functions for the entire interface are in place
 #ifdef PYBIND

src/optimization/energies.cpp

@@ -225,6 +225,15 @@ symmetric_dirichlet_energy(
   VectorX& f2energy
 );
 
+Scalar
+root_mean_square_error(const VectorX &x, const VectorX &x0);
+
+Scalar
+relative_root_mean_square_error(const VectorX &x, const VectorX &x0);
+
+Scalar
+root_mean_square_relative_error(const VectorX &x, const VectorX &x0);
+
 #ifdef PYBIND
 
 VectorX

src/optimization/energies.hh

@@ -44,6 +44,9 @@ void initialize_data_log(const std::filesystem::path& data_log_path)
     output_file << "step_size,";
     output_file << "energy,";
     output_file << "max_error,";
+    output_file << "rmse,";
+    output_file << "rrmse,";
+    output_file << "rmsre,";
     output_file << "convergence_ratio,";
     output_file << "max_change_in_metric_coords,";
     output_file << "actual_to_unconstrained_direction_ratio,";
@@ -65,6 +68,19 @@ void update_data_log(
     const VectorX& descent_direction,
     Scalar convergence_ratio)
 {
+    // Get initial metric coordinates
+    VectorX metric_init = initial_cone_metric.get_reduced_metric_coordinates();
+
+    // Get edge lengths
+    int num_reduced_edges = metric_init.size();
+    VectorX l_init(num_reduced_edges);
+    VectorX l(num_reduced_edges);
+    for (int E = 0; E < num_reduced_edges; ++E)
+    {
+        l_init[E] = exp(metric_init[E] / 2.0);
+        l[E] = exp(reduced_metric_coords[E] / 2.0);
+    }
+
     // Copy metric and make discrete
     std::unique_ptr<DifferentiableConeMetric> cone_metric =
         initial_cone_metric.set_metric_coordinates(reduced_metric_coords);
@@ -84,6 +100,9 @@ void update_data_log(
     // Compute numerics
     log.energy = opt_energy.energy(initial_cone_metric);
     log.error = sup_norm(constraint);
+    log.rmse = root_mean_square_error(l, l_init);
+    log.rrmse = relative_root_mean_square_error(l, l_init);
+    log.rmsre = root_mean_square_relative_error(l, l_init);
     log.num_flips = cone_metric->num_flips();
     log.convergence_ratio = convergence_ratio;
     log.max_change_in_metric_coords = sup_norm(change_in_metric_coords);
@@ -175,6 +194,9 @@ void write_data_log_entry(const OptimizationLog& log, const std::filesystem::pat
     output_file << std::fixed << std::setprecision(17) << log.beta << ",";
     output_file << std::fixed << std::setprecision(17) << log.energy << ",";
     output_file << std::fixed << std::setprecision(17) << log.error << ",";
+    output_file << std::fixed << std::setprecision(17) << log.rmse << ",";
+    output_file << std::fixed << std::setprecision(17) << log.rrmse << ",";
+    output_file << std::fixed << std::setprecision(17) << log.rmsre << ",";
     output_file << std::fixed << std::setprecision(17) << log.convergence_ratio << ",";
     output_file << std::fixed << std::setprecision(17) << log.max_change_in_metric_coords << ",";
     output_file << std::fixed << std::setprecision(17)

src/optimization/implicit_optimization.cpp

@@ -20,6 +20,9 @@ namespace CurvatureMetric
     Scalar beta;                                    // Step size
     Scalar energy;                                  // Energy
     Scalar error;                                   // Max angle constraint error
+    Scalar rmse;                                    // Root mean squared edge lengths
+    Scalar rrmse;                                   // Relative root mean squared edge lengths
+    Scalar rmsre;                                   // Root mean squared relative edge lengths
     int num_flips;                                  // Number of flips for MakeDelaunay
     Scalar convergence_ratio;                       // Convergence ratio
     Scalar max_change_in_metric_coords;             // Maximum change per edge in the metric coordinates

src/optimization/implicit_optimization.hh

@@ -73,4 +73,23 @@ TEST_CASE( "The energies of a triangle can be computed", "[energies]" )
     REQUIRE( float_equal(J_f2J2.coeff(0, 1), A_sq / A0_sq) );
     REQUIRE( float_equal(J_f2J2.coeff(0, 2), 0.0) );
   }
-}
+}
+
+TEST_CASE( "The errors of a triangle can be computed", "[energies]" ) {
+  SECTION ( "Errors" )
+  {
+    VectorX x(3);
+    VectorX x0(3);
+    x << 1.0, 2.0, 0.5;
+    x0 << 2.0, 2.0, 1.0;
+    // err is 1.0, 0.0, 0.5
+    // Sq err is 1.0, 0.0, 0.25
+    // Rel err is 0.5, 0.0, 0.5
+    // Rel sq err is 0.25, 0.0, 0.25
+    // total Sq err is 1.25
+    // rel factor is 9
+    CHECK( float_equal(root_mean_square_error(x, x0), std::sqrt(1.25 / 3.0)) );
+    CHECK( float_equal(relative_root_mean_square_error(x, x0), std::sqrt(1.25 / (3.0 * 9.0))) );
+    CHECK( float_equal(root_mean_square_relative_error(x, x0), std::sqrt(0.5 / 3.0)) );
+  }
+}

src/tests/test_energies.cpp

@@ -7,4 +7,4 @@
 #include "test_optimize.cpp"
 #include "test_shear.cpp"
 #include "test_refinement.cpp"
-#include "test_regression.cpp"
+#include "test_regression.cpp"