Added Gaussian kernels energy

Changelog.md

@@ -2,11 +2,12 @@
 
 ### 2.0.1
 
+- 2024 May (Bastien DOIGNIES):
+    - Fixed spectrum multipointset calculation
+    - Added gaussiankernels energy
+
 - 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):
     - Binary file format support
     - Updated documentation

doc-sources/changelog.md

@@ -2,6 +2,10 @@
 
 ### 2.0.1
 
+- 2024 May (Bastien DOIGNIES):
+    - Fixed spectrum multipointset calculation
+    - Added gaussiankernels energy
+
 - 2024 March (Bastien DOIGNIES):
     - Added semidiscrete optimal transport 2D 
 - 2024 Fev (Bastien DOIGNIES):

doc-sources/metrics/GaussianKernels.md

@@ -0,0 +1,81 @@
+# GaussianKernels
+
+## Description
+
+Compute a weighted average of distances among points. Weights are given by : $w_ij = e^{-\frac{d(x_i, x_j)^2}{2\sigma^2}}$ and correspond to a gaussian centered of variance $\sigma$ around each points.
+
+## Files
+
+```
+src/metrics/GaussianKernels.cpp  
+include/utk/metrics/GaussianKernels.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
+GaussianKernels calculator
+Usage: ./GaussianKernels [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
+  --euclidean [0]             Use euclidean distance (default is wrap around)
+  --sigma [1]                 Variance of gaussians (relative to nominal grid N^(1/D))
+
+```
+
+</div>
+
+<div class="cpp tabcontent">
+
+```  cpp
+#include <utk/utils/PointsetIO.hpp>
+#include <utk/utils/Pointset.hpp>
+#include <utk/samplers/SamplerWhitenoise.hpp>
+#include <utk/metrics/MinDist.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::MinDist(true /* Toroidal / euclidean distance */, 1.0 /* sigma */).compute(pts);
+    }
+}
+```  
+
+</div>
+
+<div class="py tabcontent">
+
+``` python
+import pyutk
+import numpy as np
+
+d = pyutk.GaussianKernels(
+  True, # Toroidal / Euclidean Distance, 
+  1.0   # Sigma (relative to nominal grid !)
+).compute(np.random.uniform(0, 1, (128, 2)))
+```  
+
+</div>
+
+## License
+
+BSD, see `MinDist.hpp`

include/utk/metrics/GaussianKernels.hpp

@@ -0,0 +1,129 @@
+/*
+ * Coded by Hélène Perrier helene.perrier@liris.cnrs.fr
+ * and Bastien Doignies bastien.doignies@liris.cnrs.fr 
+ * and David Coeurjolly David.coeurjolly@liris.cnrs.fr
+ *
+ * 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
+
+#ifdef UTK_USE_OPENMP
+    #include <omp.h>
+#endif
+
+#include <limits>
+#include <utk/utils/Pointset.hpp>
+#include <utk/utils/Omp.hpp>
+#include <cmath>
+
+namespace utk
+{
+    class GaussianKernels
+    {
+    public:
+        GaussianKernels(bool to = false, double _sigma = 1.0) :
+            toroidal(to), sigma(_sigma)
+        {}
+
+        void setToroidal(bool to) { toroidal = to; }
+        void setSigma(double _sigma) { sigma = _sigma; }
+
+        template<typename T>
+        T compute(const Pointset<T>& pts) const
+        {
+            const auto dist = toroidal ? &GaussianKernels::toricDistanceSquared<T> :
+                                         &GaussianKernels::distanceSquared<T>;
+
+
+            const double PI = 3.141592653589793238462;
+            const unsigned int N  = pts.Npts();
+            const unsigned int D  = pts.Ndim();
+            const double scaledSigma = sigma * std::pow(N, -1.0 / D);
+            const double invSigma = - 1.0 / (2.0 * scaledSigma * scaledSigma); 
+            const double norm     = PI * scaledSigma * scaledSigma / (2 * N);
+            // const double invN     = 1.0 / N;
+
+            double value = 0.;
+
+            #pragma omp parallel for reduction(+: value)
+            for (OPENMP_UINT k = 0; k < N; k++)
+            {
+                for (OPENMP_UINT l = 0; l < N; l++)
+                {
+                    const double weight = std::exp(invSigma * dist(pts.Ndim(), pts[k], pts[l]));
+                    value += weight * norm;
+                }   
+            }
+
+            return value;
+        }
+
+        template<typename T>
+        std::vector<T> compute(const std::vector<Pointset<T>>& ptss) const
+        {
+            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:
+        bool toroidal;
+        double sigma;
+
+        template <typename T>
+        static T toricDistanceSquared(uint32_t D, const T* a, const T* b)
+        {
+            T dist = 0;
+            for(uint32_t d = 0; d < D; d++)
+            {
+                T tmp = std::abs(a[d] - b[d]);
+                T diff = std::min(tmp, 1 - tmp);
+                dist += diff * diff;
+            }
+            return dist;
+        }
+
+        template <typename T>
+        static T distanceSquared(uint32_t D, const T* a, const T* b)
+        {
+            T dist = 0;
+            for(uint32_t d = 0; d < D; d++)
+            {
+                T diff = a[d] - b[d];
+                dist += diff * diff;
+            }
+            return dist;
+        }
+    };
+};

mkdocs.yml

@@ -104,6 +104,7 @@ nav:
     - CenteredL2: 'metrics/CenteredL2.md'
     - Diaphony: 'metrics/Diaphony.md'
     - MinDist: 'metrics/MinDist.md'
+    - GaussianKernels: 'metrics/GaussianKernels.md'
     - t-value: 'metrics/tvalue.md'
     - PCF: 'metrics/PCF.md'
     - Spectrum: 'metrics/Spectrum.md'

src/metrics/CMakeLists.txt

@@ -1,6 +1,7 @@
 SET(Metrics
 
     MinDist   
+    GaussianKernels
 
     TValue
 

src/metrics/GaussianKernels.cpp

@@ -0,0 +1,58 @@
+/*
+ * Coded by Hélène Perrier helene.perrier@liris.cnrs.fr
+ * and Bastien Doignies bastien.doignies@liris.cnrs.fr 
+ * and David Coeurjolly David.coeurjolly@liris.cnrs.fr
+ *
+ * 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/GaussianKernels.hpp>
+
+int main(int argc, char** argv)
+{
+    CLI::App app { "GaussianKernels calculator" };
+    auto* margs = utk::add_arguments(app);    
+
+    bool euclidean = false;
+    double sigma = 1.0;
+    app.add_flag("--euclidean", euclidean, "Use euclidean distance (default is wrap around)")->default_val(euclidean);
+    app.add_option("--sigma", sigma, "Variance of gaussians (relative to nominal grid N^(1/D))")->default_val(sigma);
+
+    CLI11_PARSE(app, argc, argv);
+
+    auto ptss = margs->GetAllPointsets();
+    utk::CheckPointsets(ptss);
+
+    auto rslts = utk::GaussianKernels(euclidean, sigma).compute(ptss);
+
+    auto& ostream = margs->GetOutputStream();
+    for (double rslt : rslts)
+        ostream << rslt << '\n';
+
+    delete margs;
+}

src/metrics/python/BaseMetrics.cpp

@@ -44,6 +44,7 @@
 #include <utk/metrics/RadialSpectrum.hpp>
 #include <utk/metrics/PCF.hpp>
 #include <utk/metrics/MinDist.hpp>
+#include <utk/metrics/GaussianKernels.hpp>
 
 #include <utk/metrics/TValue.hpp>
 
@@ -151,9 +152,14 @@ void init_Metrics(py::module& m)
         .def("compute", GetComputeFunction<BoundariesStarDiscrepancy, double, std::pair<double, double>>());
 
     py::class_<MinDist>(m, "MinDistance")
-        .def(py::init<double>(), py::arg("toroidal") = false)
+        .def(py::init<bool>(), py::arg("toroidal") = false)
         .def("compute", GetComputeFunction<MinDist, double>());
 
+    py::class_<GaussianKernels>(m, "GaussianKernels")
+        .def(py::init<bool, double>(), py::arg("toroidal") = false, py::arg("sigma") = 1.0)
+        .def("compute", GetComputeFunction<GaussianKernels, double>());
+
+
     py::class_<TValue>(m, "TValue")
         .def(py::init<std::uint8_t>(), py::arg("basis") = 2)
         .def("compute", GetComputeFunction<TValue, double, std::uint32_t>());