Commit current state of repo

Author: ManifoldFR

CMakeLists.txt

@@ -29,7 +29,7 @@ set(DOXYGEN_USE_MATHJAX YES)
 
 if(UNIX)
 
-    add_compile_options(-Wall -march=native)
+    add_compile_options(-Wall -Wpedantic -Wextra -march=native)
 
     add_library(${PROJECT_NAME} SHARED ${EntropicMFG_SOURCES} ${EntropicMFG_HEADERS})
     set_target_properties(${PROJECT_NAME} PROPERTIES LINKER_LANGUAGE CXX)

include/Kernels.h

@@ -21,7 +21,7 @@ class KernelBase {
 using KernelPtr = std::shared_ptr<KernelBase<double>>;
 
 
-/// Heat kernel on Euclidean space.
+/// Gaussian (heat) kernel on Euclidean space.
 template <size_t Dim, typename S = double>
 class EuclideanHeatKernelTpl : public KernelBase<S> {
     typedef S Scalar;
@@ -40,10 +40,11 @@ class EuclideanHeatKernelTpl : public KernelBase<S> {
 };
 
 
-
-
+/// Separable Gaussian kernel on Euclidean space.
+/// Assumes a kernel \f$ K = K_1\otimes K_2\f$.
 class SeparableEuclideanKernel2D : public EuclideanHeatKernelTpl<2, double> {
 public:
+    /// Individual sub-kernels.
     MatrixXd K1_, K2_;
     double dx, dy;
     SeparableEuclideanKernel2D(

include/RectangularGridLaplacian.h

@@ -0,0 +1,13 @@
+#pragma once
+
+#include "Kernels.h"
+#include "RectangularGridLaplacian.h"
+
+
+namespace kernels
+{
+
+
+
+
+}

include/kernels/Geometric.h

@@ -0,0 +1,94 @@
+#pragma once
+
+#include <Eigen/Core>
+#include <Eigen/SparseCore>
+#include <Eigen/SparseCholesky>
+
+
+using namespace Eigen;
+
+namespace kernels
+{
+    namespace laplacian
+    {
+
+    /// Define a matrix for the 1D Laplacian using the central
+    /// differences stencil.
+    auto make_laplacian_matrix1d(size_t nx, double h)
+    {
+        double h_2 = h * h;
+        SparseMatrix<double> out(nx, nx);
+
+        for (size_t i = 0; i < nx; i++) {
+            out.coeffRef(i, i) = 2. / h_2;
+            if (i > 0) {
+                out.coeffRef(i, i-1) = 1. / h_2;
+            }
+            if (i < nx) {
+                out.coeffRef(i, i+1) = 1. / h_2;
+            }
+        }
+        out.makeCompressed();
+        return out;
+    }
+
+    /// Given a mask of the discretization grid for the domain, define the graph obtained by
+    /// removing the masked nodes from the classical grid graph.
+    ///
+    /// Our convention is that a 1 in the mask indicates an obstacle.
+    /// The \a x coord index is \c i.
+    struct MaskedStencilGraph
+    {
+        /// Array type for grid mask.
+        typedef Array<bool, -1, -1> Mask_t;
+        Mask_t mask_;
+        size_t n_nodes;
+        std::vector<std::vector<size_t>> neighbors_;
+    
+        MaskedStencilGraph(const Mask_t& mask
+        ) : mask_(mask) {
+            size_t nx = mask.cols();
+            size_t ny = mask.rows();
+            n_nodes = 0;
+            int v_id;
+            for (size_t i = 0; i < nx; i++)
+            {
+                for (size_t j = 0; j < ny; j++)
+                {
+                    if (mask_(i, j) == 0) {
+                        std::vector<size_t> nlist { n_nodes };
+                        n_nodes++;
+                        neighbors_.push_back(nlist);
+                    }
+                }
+            }
+
+            // 
+            for (size_t i = 0; i < nx; i++)
+            {
+                for (size_t j = 0; j < ny; j++)
+                {
+                    if (mask_(i, j) == 0)
+                    {
+                        std::vector<size_t> nlist{n_nodes};
+                        n_nodes++;
+                        neighbors_.push_back(nlist);
+                    }
+                }
+            }
+        };
+
+        /// @param i row index
+        /// @param j columns index
+        size_t grid_idx_to_node_id(size_t i, size_t j) {
+            size_t nx = mask_.cols();
+            size_t ny = mask_.rows();
+            int v = j * nx + i;
+            return v;
+        }
+
+    };
+
+    }
+
+}

include/kernels/RectangularGridLaplacian.h

@@ -49,9 +49,9 @@ bool MultimarginalSinkhorn::solve(size_t num_iterations, std::vector<MatrixXd> &
         std::cout << BLUE << "> MultiSinkhorn solver" << RESET << std::endl;
 
     nsteps_ = potentials.size();
-    metric_vals_.resize(0);
-    _old_potentials.resize(nsteps_);
-    potentials_.resize(nsteps_);
+    metric_vals_.reserve(nsteps_);
+    _old_potentials.reserve(nsteps_);
+    potentials_.reserve(nsteps_);
 
     for (size_t t = 0; t < nsteps_; t++) {
         potentials_[t].noalias() = potentials[t];