Semidiscrete ot

Changelog.md

@@ -2,6 +2,9 @@
 
 ### 2.0.1
 
+- 2024 March (Bastien DOIGNIES):
+    - Added semidiscrete optimal transport 2D 
+
 - 2024 Feb (David Coeurjolly)
     - cmake CPM instead of `FETCH_CONTENT` for dependencies 
 - 2024 Fev (Bastien DOIGNIES):

cmake/externals.cmake

@@ -1 +1,2 @@
-include(cmake/externals/sot.cmake)
+include(cmake/externals/sot.cmake)
+include(cmake/externals/semidiscrete_ot_2d.cmake)

cmake/externals/semidiscrete_ot_2d.cmake

@@ -0,0 +1,15 @@
+IF (UTK_USE_OPENMP)
+    include(CPM)
+
+    set(semidiscrete_ot_2d_patch git apply ${CMAKE_CURRENT_SOURCE_DIR}/cmake/externals/semidiscrete_ot_2d.patch)
+    CPMAddPackage(
+        NAME SEMIDISCRETE_OT_2D
+        GITHUB_REPOSITORY nbonneel/semi_discrete_ot_2d
+        GIT_TAG af331f8f1058291e359e51830fc3fe5c7009b2f4
+        UPDATE_DISCONNECTED 1 
+        PATCH_COMMAND ${semidiscrete_ot_2d_patch}
+    )
+
+    SET(UTK_EXTERNALS_SEMIDISCRETE_OT semidiscrete_ot_2d)
+    add_definitions(-DUTK_EXTERNALS_SEMIDISCRETE_OT)
+ENDIF()

cmake/externals/semidiscrete_ot_2d.patch

@@ -0,0 +1,62 @@
+From 645b95cccd1d341a45a9f1ba51aa4f2010c7e048 Mon Sep 17 00:00:00 2001
+From: Bastien DOIGNIES <[email protected]>
+Date: Fri, 8 Mar 2024 14:08:09 +0100
+Subject: [PATCH] Added CMakeLists.txt
+
+---
+ CMakeLists.txt | 41 +++++++++++++++++++++++++++++++++++++++++
+ 1 file changed, 41 insertions(+)
+ create mode 100644 CMakeLists.txt
+
+diff --git a/CMakeLists.txt b/CMakeLists.txt
+new file mode 100644
+index 0000000..11a90c3
+--- /dev/null
++++ b/CMakeLists.txt
+@@ -0,0 +1,41 @@
++cmake_minimum_required(VERSION 3.12)
++
++project(semidiscrete_ot_2d)
++
++set(CMAKE_CXX_STANDARD 17)
++set(CMAKE_CXX_STANDARD_REQUIRED ON)
++set(CMAKE_CXX_EXTENSIONS ON)
++
++# OpenMP management
++if(APPLE)
++    message("-- Setting OpenMP flags on MacOs. Assuming `brew install libomp`")
++    if(CMAKE_C_COMPILER_ID MATCHES "Clang\$")
++    IF(EXISTS "/opt/homebrew/")
++        message(STATUS "Using M1-M2/Homebrew C paths")
++        set(OpenMP_C_FLAGS "-Xpreprocessor -fopenmp  -I/opt/homebrew/opt/libomp/include" CACHE INTERNAL "OpenMP flags for #Xcode toolchain.")
++        set(OpenMP_C_LIB_NAMES "omp" CACHE INTERNAL "OpenMP lib name for Xcode toolchain.")
++        set(OpenMP_omp_LIBRARY "/opt/homebrew/opt/libomp/lib/libomp.dylib" CACHE INTERNAL "OpenMP lib name for Xcode toolchain.")
++      else()
++        set(OpenMP_C_FLAGS "-Xpreprocessor -fopenmp -I/usr/local/opt/libomp/include" CACHE INTERNAL "OpenMP flags for Xcode toolchain.")
++        set(OpenMP_C_LIB_NAMES "omp" CACHE INTERNAL "OpenMP lib name for Xcode toolchain.")
++        set(OpenMP_omp_LIBRARY "/usr/local/opt/libomp/lib/libomp.dylib" CACHE INTERNAL "OpenMP lib name for Xcode toolchain.")
++      endif()
++    endif()
++   if(CMAKE_CXX_COMPILER_ID MATCHES "Clang\$")
++      IF(EXISTS "/opt/homebrew/")
++        message(STATUS "Using M1-M2/Homebrew C++ paths")
++        set(OpenMP_CXX_FLAGS "-Xpreprocessor -fopenmp  -I/opt/homebrew/opt/libomp/include"  CACHE INTERNAL "OpenMP flags for Xcode toolchain.")
++        set(OpenMP_CXX_LIB_NAMES "omp" CACHE INTERNAL "OpenMP lib name for Xcode toolchain.")
++        set(OpenMP_omp_LIBRARY "/opt/homebrew/opt/libomp/lib/libomp.dylib" CACHE INTERNAL "OpenMP lib name for Xcode toolchain.")
++      else()
++        set(OpenMP_CXX_FLAGS "-Xpreprocessor -fopenmp -I/usr/local/opt/libomp/include" CACHE INTERNAL "OpenMP flags for Xcode toolchain.")
++        set(OpenMP_CXX_LIB_NAMES "omp" CACHE INTERNAL "OpenMP lib name for Xcode toolchain.")
++        set(OpenMP_omp_LIBRARY "/usr/local/opt/libomp/lib/libomp.dylib" CACHE INTERNAL "OpenMP lib name for Xcode toolchain.")
++      endif()
++    endif()
++endif()
++
++find_package(OpenMP REQUIRED)
++
++add_library(semidiscrete_ot_2d transport.cpp transport.h)
++target_link_libraries(semidiscrete_ot_2d PRIVATE OpenMP::OpenMP_CXX)
++target_include_directories(semidiscrete_ot_2d PUBLIC ${CMAKE_CURRENT_LIST_DIR})
+\ No newline at end of file
+-- 
+2.39.2
+

doc-sources/changelog.md

@@ -2,6 +2,8 @@
 
 ### 2.0.1
 
+- 2024 March (Bastien DOIGNIES):
+    - Added semidiscrete optimal transport 2D 
 - 2024 Fev (Bastien DOIGNIES):
     - Binary file format support
     - Updated documentation

doc-sources/metrics/SemidiscreteOT.md

@@ -0,0 +1,77 @@
+# Semidiscrete Optimal transport to uniform distribution
+
+## Description
+
+Compute the semidiscrete optimal transport to uniform distribution distance of a pointset. 
+
+The code originates from : [https://github.com/nbonneel/semi_discrete_ot_2d](https://github.com/nbonneel/semi_discrete_ot_2d). A patch is applied so that it is compatible with CMake. 
+
+This metrics is only available for 2D pointsets . 
+
+## Files
+
+```
+src/metrics/SemidiscreteOT.cpp  
+include/utk/metrics/SemidiscreteOT.hpp
+```
+
+## Usage
+
+<button class="tablink exebutton" onclick="openCode('exe', this)" markdown="1">Exe</button> 
+<button class="tablink cppbutton" onclick="openCode('cpp', this)" markdown="1">C++</button> 
+<button class="tablink pybutton" onclick="openCode('py', this)" markdown="1">Python</button> 
+<br/>
+
+
+<div class="exe tabcontent">
+
+```bash
+Semidiscrete Optimal Transport to uniform distribution calculator (2D)
+Usage: ./SemidiscreteOT [OPTIONS]
+
+Options:
+  -h,--help                   Print this help message and exit
+  -i,--input TEXT:FILE ... REQUIRED
+                              Input file(s)
+  -o,--output TEXT            Output file (empty is stdout)
+  --silent                    Silence UTK logs
+```
+
+</div>
+
+<div class="cpp tabcontent">
+
+```  cpp
+#include <utk/utils/PointsetIO.hpp>
+#include <utk/utils/Pointset.hpp>
+#include <utk/samplers/SamplerWhitenoise.hpp>
+#include <utk/metrics/SemidiscreteOT.hpp>
+
+int main()
+{
+    utk::Pointset<double> pts;
+
+    utk::SamplerWhitenoise wn(2 /* dimension */);
+    wn.setRandomSeed(/* empty means random, can also pass a number */);
+    // Check for no errors
+    if (wn.generateSamples(pts, 1024 /* Number of points */))
+    {
+        auto rslt = utk::SemidiscreteOT().compute(pts);
+    }
+}
+```  
+
+</div>
+
+<div class="py tabcontent">
+
+``` python
+import pyutk
+import numpy as np
+
+distance = pyutk.SemidiscreteOT().compute(np.random.uniform(0, 1, (128, 2)))
+```  
+
+</div>
+
+## License

include/utk/metrics/SemidiscreteOT.hpp

@@ -0,0 +1,89 @@
+/*
+ * Coded by Nicolas Bonneel [email protected] 
+ *  and Bastien Doignies [email protected] 
+ *
+ * Copyright (c) 2023 CNRS Université de Lyon
+ * All rights reserved.
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies,
+ * either expressed or implied, of the UTK project.
+ */
+#pragma once
+
+
+#include <limits>
+#include <utk/utils/Pointset.hpp>
+#include <utk/utils/Omp.hpp>
+#include <utk/utils/Log.hpp>
+#include <cmath>
+
+#include "transport.h"
+
+namespace utk
+{
+    class SemidiscreteOT
+    {
+    public:
+        SemidiscreteOT()
+        {}
+
+        void setMaxNewtonIter(uint32_t it) { maxNewtonIters = it; }
+        void setWorstAreaRelativeThreshold(double ar) { worstAreaRelativeThresholdPercent = ar; }
+
+        template<typename T>
+        T compute(const Pointset<T>& pts)
+        {
+            if (pts.Ndim() != 2)
+            {
+                UTK_ERROR("Only 2D pointsets are supported. Found {}-d pointsets", pts.Ndim());
+                return static_cast<T>(-1.0);
+            }
+
+            transport::OptimalTransport2D ot;
+
+            ot.V.vertices.resize(pts.Npts());
+            for (unsigned int i = 0; i < pts.Npts(); i++)
+                ot.V.vertices[i] = transport::Vector(pts[i][0], pts[i][1]);
+
+            return static_cast<T>(ot.optimize(maxNewtonIters, worstAreaRelativeThresholdPercent));
+        }
+
+        template<typename T>
+        std::vector<T> compute(const std::vector<Pointset<T>>& ptss)
+        {
+            if (ptss.size() == 0) return std::vector<T>();
+
+            std::vector<T> results(ptss.size());
+            for (uint32_t i = 0; i < results.size(); i++)
+            {
+                results[i] = compute(ptss[i]);
+            }
+
+            return results;
+        }
+    private:
+        uint32_t maxNewtonIters = 100;                  // Default value of original code !
+        double worstAreaRelativeThresholdPercent = 0.5; // Default value of original code !
+    };
+};

src/metrics/CMakeLists.txt

@@ -17,7 +17,8 @@ SET(Metrics
     PCF
 )
 
-SET(pyutk_metrics)
+SET(pytuk_metrics)
+SET(pyutk_metrics_includes)
 SET(pyutk_metrics_libs pybind11::module spdlog::spdlog)
 
 IF (OpenMP_FOUND)
@@ -45,6 +46,25 @@ FOREACH(FILE ${Metrics})
     ENDIF()
 ENDFOREACH(FILE)
 
+# Semidiscrete optimal transport
+IF (UTK_EXTERNALS_SEMIDISCRETE_OT)
+    IF (UTK_BUILD_EXE)
+        add_executable(${UTK_EXE_PREFIX}SemidiscreteOT SemidiscreteOT.cpp)
+        target_include_directories(${UTK_EXE_PREFIX}SemidiscreteOT PRIVATE "${PROJECT_SOURCE_DIR}/include/")
+        target_include_directories(${UTK_EXE_PREFIX}SemidiscreteOT PRIVATE "${semidiscrete_ot_2d_SOURCE_DIR}/")
+
+        target_link_libraries(${UTK_EXE_PREFIX}SemidiscreteOT PRIVATE semidiscrete_ot_2d spdlog::spdlog CLI11::CLI11 OpenMP::OpenMP_CXX)
+
+        install(TARGETS ${UTK_EXE_PREFIX}SemidiscreteOT
+            RUNTIME DESTINATION bin
+            LIBRARY DESTINATION lib
+            ARCHIVE DESTINATION lib)
+    ENDIF()
+
+    LIST(APPEND pyutk_metrics_libs semidiscrete_ot_2d OpenMP::OpenMP_CXX)
+    LIST(APPEND pyutk_metrics_includes "${semidiscrete_ot_2d_SOURCE_DIR}/")
+ENDIF()
+
 IF(UTK_BUILD_EXE)
     add_subdirectory(Integration)
 ENDIF()
@@ -53,6 +73,7 @@ IF (UTK_PYTHON)
     # TODO : OpenMP !!!
     SET(pyutk_metrics_bindings 
         "python/Metrics.cpp"
+        "python/ExternalsMetrics.cpp"
         "python/BaseMetrics.cpp"
     )
     add_library(pyutkMetrics ${pyutk_metrics} ${pyutk_metrics_bindings})
@@ -68,4 +89,5 @@ IF (UTK_PYTHON)
 
     target_include_directories(pyutkMetrics PRIVATE "${UTK_INCLUDE_DIR}")
     target_link_libraries(pyutkMetrics PRIVATE ${pyutk_metrics_libs})
+    target_include_directories(pyutkMetrics PRIVATE ${pyutk_metrics_includes})
 ENDIF()

src/metrics/SemidiscreteOT.cpp

@@ -0,0 +1,51 @@
+/*
+ * Coded by Nicolas Bonneel [email protected] 
+ *  and Bastien Doignies [email protected] 
+ *
+ * Copyright (c) 2023 CNRS Université de Lyon
+ * All rights reserved.
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies,
+ * either expressed or implied, of the UTK project.
+ */
+#include <utk/utils/MetricsArgumentParser.hpp>
+#include <utk/metrics/SemidiscreteOT.hpp>
+
+int main(int argc, char** argv)
+{
+    CLI::App app { "Semidiscrete Optimal Transport to uniform distribution calculator (2D)" };
+    auto* margs = utk::add_arguments(app);    
+    CLI11_PARSE(app, argc, argv);
+
+    auto ptss = margs->GetAllPointsets();
+    utk::CheckPointsets(ptss);
+
+    auto rslts = utk::SemidiscreteOT().compute(ptss);
+
+    auto& ostream = margs->GetOutputStream();
+    for (double rslt : rslts)
+        ostream << rslt << '\n';
+
+    delete margs;
+}

src/metrics/python/BaseMetrics.cpp

@@ -52,7 +52,6 @@
 #include <utk/metrics/Integrands/HeavisideIntegrands.hpp>
 #include <utk/metrics/Integrands/BlinnPhong.hpp>
 
-
 template<typename InputType, typename IntegrandType>
 auto GetBuildComputeFunction()
 {
@@ -84,7 +83,6 @@ auto GetBuildComputeFunction()
     };
 }
 
-
 void init_Metrics(py::module& m)
 {
     using namespace utk;