Minor support improvements and fix non-determinism introduced in comp…

…iler fixes. (#32)

* Add support for `debug` logger level.

* Limit precision of progress log.

* Add a test for deterministic behavior.

* Improve determinism.

public/coacd.cpp

@@ -90,6 +90,8 @@ std::vector<Mesh> CoACD(Mesh const &input, double threshold,
 void set_log_level(std::string_view level) {
   if (level == "off") {
     logger::get()->set_level(spdlog::level::off);
+  } else if (level == "debug") {
+    logger::get()->set_level(spdlog::level::debug);
   } else if (level == "info") {
     logger::get()->set_level(spdlog::level::info);
   } else if (level == "warn" || level == "warning") {

run_tests.py

@@ -2,6 +2,7 @@
 import os
 import glob
 import coacd
+import numpy
 import unittest
 
 
@@ -17,6 +18,23 @@ def _test(self: 'TestExamples'):
 
         return _test
 
+    @unittest.skip
+    def test_deterministic(self):
+        mesh = trimesh.load('examples/SnowFlake.obj', force="mesh")
+        coacd.set_log_level("debug")
+        mesh = coacd.Mesh(mesh.vertices, mesh.faces)
+        acd1 = coacd.run_coacd(mesh, merge=False, seed=1337)
+        acd2 = coacd.run_coacd(mesh, merge=False, seed=1337)
+        self.assertEqual(len(acd1), len(acd2))
+        acd1 = sorted(acd1, key=lambda x: (len(x[0]), len(x[1]), x[0].sum()))
+        acd2 = sorted(acd2, key=lambda x: (len(x[0]), len(x[1]), x[0].sum()))
+        for (v1, t1), (v2, t2) in zip(acd1, acd2):
+            self.assertEqual(len(v1), len(v2))
+            self.assertEqual(len(t1), len(t2))
+        for (v1, t1), (v2, t2) in zip(acd1, acd2):
+            self.assertTrue(numpy.allclose(numpy.sort(v1), numpy.sort(v2)))
+            self.assertTrue(numpy.allclose(numpy.sort(t1), numpy.sort(t2)))
+
 
 for f in glob.glob("examples/*.obj"):
     setattr(TestExamples, f'test_{os.path.splitext(os.path.basename(f))[0].lower()}', TestExamples.single(f))

src/clip.cpp

@@ -1,4 +1,5 @@
 #include "clip.h"
+#include "process.h"
 #include "include/CDTUtils.h"
 #include "include/CDT.h"
 

src/mcts.cpp

@@ -1,6 +1,6 @@
 #include <stdint.h>
-#include <random>
 #include "mcts.h"
+#include "process.h"
 
 namespace coacd
 {
@@ -629,9 +629,7 @@ namespace coacd
 
         if (shuffle)
         {
-            std::random_device rd;
-            std::mt19937 g(rd());
-            std::shuffle(planes.begin(), planes.end(), g);
+            std::shuffle(planes.begin(), planes.end(), coacd::random_engine);
         }
     }
 
@@ -854,7 +852,6 @@ namespace coacd
         initial_mesh.ComputeCH(initial_ch);
         double cost = ComputeRv(initial_mesh, initial_ch, params.rv_k) / params.mcts_max_depth;
         vector<Plane> current_path;
-        srand(params.seed);
 
         for (int i = 0; i < computation_budget; i++)
         {

src/model_obj.cpp

@@ -1,5 +1,6 @@
 #include <stdint.h>
 #include "model_obj.h"
+#include "process.h"
 #include "quickhull/QuickHull.hpp"
 #include "btConvexHull/btConvexHullComputer.h"
 #include "nanoflann.hpp"
@@ -315,7 +316,6 @@ namespace coacd
         if (base != 0)
             resolution = size_t(max(1000, int(resolution * (aObj / base))));
 
-        srand(seed);
         for (int i = 0; i < (int)triangles.size(); i++)
         {
             if (flag && plane.Side(points[triangles[i][0]], 1e-3) == 0 &&
@@ -331,16 +331,18 @@ namespace coacd
             else
                 N = max(int(i % 2 == 0), int(resolution / aObj * area));
 
-            int seed = rand() % 1000;
+            std::uniform_int_distribution<int> seeder(0, 1000);
+            int seed = seeder(coacd::random_engine);
             float r[2];
             for (int k = 0; k < N; k++)
             {
                 double a, b;
                 if (k % 3 == 0)
                 {
+                    std::uniform_real_distribution<double> uniform(0.0, 1.0);
                     //// random sample
-                    a = rand() / double(RAND_MAX);
-                    b = rand() / double(RAND_MAX);
+                    a = uniform(coacd::random_engine);
+                    b = uniform(coacd::random_engine);
                 }
                 else
                 {

src/process.cpp

@@ -8,6 +8,7 @@
 
 namespace coacd
 {
+    thread_local std::mt19937 random_engine;
 
     void ManifoldPreprocess(Params &params, Model &m)
     {
@@ -220,8 +221,9 @@ namespace coacd
 #endif
             for (int p = 0; p < (int)InputParts.size(); p++)
             {
+                random_engine.seed(params.seed);
                 if (p % ((int)InputParts.size() / 10 + 1) == 0)
-                    logger::info("Processing [{}%]", p * 100.0 / (int)InputParts.size());
+                    logger::info("Processing [{:.1f}%]", p * 100.0 / (int)InputParts.size());
 
                 Model pmesh = InputParts[p], pCH;
                 Plane bestplane;
@@ -286,7 +288,7 @@ namespace coacd
 #endif
                 }
             }
-            logger::info("Processing [100%]");
+            logger::info("Processing [100.0%]");
             InputParts.clear();
             InputParts = tmp;
             tmp.clear();

src/process.h

@@ -2,6 +2,7 @@
 
 #include <iostream>
 #include <string>
+#include <random>
 #include <fstream>
 #include <vector>
 #include <math.h>
@@ -19,6 +20,7 @@
 
 namespace coacd
 {
+  extern thread_local std::mt19937 random_engine;
 
   void ManifoldPreprocess(Params &params, Model &m);
   void MergeCH(Model &ch1, Model &ch2, Model &ch);