Merge pull request #99 from ecmwf/feature/qhull

Feature/qhull
ecmwf · Jan 31, 2024 · 320627f · 320627f
2 parents 547577f + 19701fe
commit 320627f
Show file tree

Hide file tree

Showing 9 changed files with 755 additions and 21 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -143,14 +143,32 @@ ecbuild_add_option( FEATURE AEC
 ecbuild_add_option( FEATURE XXHASH
                     DESCRIPTION "xxHash support for hashing" )
 
-### Codec options
+### eckit::codec
 
 ecbuild_add_option( FEATURE ECKIT_CODEC
                     DEFAULT OFF
-                    DESCRIPTION "Encoding/Decoding library ported from atlas_io" )
+                    DESCRIPTION "eckit::codec encoding/decoding library" )
 
 set( eckit_CODEC_STATIC_ASSERT ON CACHE BOOL "eckit::codec static assertions" )
 
+### eckit::maths
+
+ecbuild_add_option( FEATURE CONVEX_HULL
+                    DEFAULT OFF
+                    DESCRIPTION "eckit::maths library convex hull/Delaunay triangulation" )
+                    # REQUIRED_PACKAGES "Qhull COMPONENTS C CXX"
+
+if( eckit_HAVE_CONVEX_HULL )
+    find_package( Qhull REQUIRED CONFIG )
+
+    if( NOT TARGET Qhull::qhullcpp OR NOT TARGET Qhull::qhullstatic_r )
+        message( FATAL_ERROR "eckit::maths ENABLE_CONVEX_HULL requires Qhull C/C++ libraries" )
+    endif()
+
+    add_library(Qhull::Qhull INTERFACE IMPORTED)
+    target_link_libraries(Qhull::Qhull INTERFACE Qhull::qhullcpp Qhull::qhullstatic_r )
+endif()
+
 ### LAPACK
 
 if( eckit_HAVE_MKL )

diff --git a/src/eckit/maths/CMakeLists.txt b/src/eckit/maths/CMakeLists.txt
@@ -1,16 +1,24 @@
-list( APPEND eckit_maths_lib_srcs
-Eigen.h
-Lapack.h
-Lapack.cc
-Matrix.h
-MatrixEigen.h
-MatrixLapack.h
-)
+list(APPEND eckit_maths_private_libs "${LAPACK_LIBRARIES}" "${BLAS_LIBRARIES}")
+list(APPEND eckit_maths_sources
+    Eigen.h
+    Lapack.h
+    Lapack.cc
+    Matrix.h
+    MatrixEigen.h
+    MatrixLapack.h)
+
+if(eckit_HAVE_CONVEX_HULL)
+    list(APPEND eckit_maths_sources ConvexHull.h ConvexHullN.h Qhull.cc Qhull.h)
+    list(APPEND eckit_maths_private_libs Qhull::Qhull)
+endif()
+
+ecbuild_add_library(
+    TARGET              eckit_maths
+    TYPE                SHARED
+    INSTALL_HEADERS     ALL
+    HEADER_DESTINATION  ${INSTALL_INCLUDE_DIR}/eckit/maths
+    SOURCES             ${eckit_maths_sources}
+    PRIVATE_LIBS        ${eckit_maths_private_libs}
+    PUBLIC_INCLUDES     "${EIGEN3_INCLUDE_DIRS}"
+    PUBLIC_LIBS         eckit)
 
-ecbuild_add_library( TARGET              eckit_maths TYPE SHARED
-                     INSTALL_HEADERS     ALL
-                     HEADER_DESTINATION  ${INSTALL_INCLUDE_DIR}/eckit/maths
-                     SOURCES             ${eckit_maths_lib_srcs}
-                     PUBLIC_INCLUDES     "${EIGEN3_INCLUDE_DIRS}"
-                     PUBLIC_LIBS         eckit
-                     PRIVATE_LIBS        "${LAPACK_LIBRARIES}" "${BLAS_LIBRARIES}" )
diff --git a/src/eckit/maths/ConvexHull.h b/src/eckit/maths/ConvexHull.h
@@ -0,0 +1,86 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+
+#pragma once
+
+#include <vector>
+
+#include "eckit/exception/Exceptions.h"
+#include "eckit/maths/Qhull.h"
+
+
+namespace eckit::maths {
+
+
+class ConvexHull {
+public:
+    // -- Types
+
+    using coord_t    = typename Qhull::coord_t;
+    using facets_n_t = typename Qhull::facets_n_t;
+
+    struct Exception : eckit::Exception {
+        Exception(const std::string& what, int _errorCode, const std::string& _command) :
+            eckit::Exception(what), errorCode(_errorCode), command(_command) {}
+        const int errorCode;
+        const std::string& command;
+    };
+
+    struct DimensionError : Exception {
+        using Exception::Exception;
+    };
+
+    struct InputError : Exception {
+        using Exception::Exception;
+    };
+
+    struct OptionError : Exception {
+        using Exception::Exception;
+    };
+
+    struct PrecisionError : Exception {
+        using Exception::Exception;
+    };
+
+    struct TopologyError : Exception {
+        using Exception::Exception;
+    };
+
+    // -- Constructors
+
+    ConvexHull(const ConvexHull&) = delete;
+    ConvexHull(ConvexHull&&)      = delete;
+
+    // -- Destructor
+
+    virtual ~ConvexHull() = default;
+
+    // -- Operators
+
+    void operator=(const ConvexHull&) = delete;
+    void operator=(ConvexHull&&)      = delete;
+
+    // -- Methods
+
+    virtual std::vector<size_t> list_vertices() const            = 0;
+    virtual std::vector<std::vector<size_t>> list_facets() const = 0;
+
+    virtual facets_n_t facets_n() const                = 0;
+    virtual std::vector<size_t> facets(size_t n) const = 0;
+
+protected:
+    // -- Constructors
+
+    ConvexHull() /*noexcept*/ = default;
+};
+
+
+}  // namespace eckit::maths
diff --git a/src/eckit/maths/ConvexHullN.h b/src/eckit/maths/ConvexHullN.h
@@ -0,0 +1,66 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+
+#pragma once
+
+#include <algorithm>
+
+#include "eckit/maths/ConvexHull.h"
+
+
+namespace eckit::maths {
+
+
+template <size_t N>
+class ConvexHullN : public ConvexHull {
+public:
+    explicit ConvexHullN(const ConvexHull::coord_t& coord, const std::string& qhull_command = Qhull::COMMAND_DEFAULT) :
+        qhull_(N, coord, qhull_command) {}
+
+    explicit ConvexHullN(const std::vector<std::vector<double>>& coord_v,
+                         const std::string& qhull_command = Qhull::COMMAND_DEFAULT) :
+        ConvexHullN(convert_vector_v(coord_v), qhull_command) {}
+
+    explicit ConvexHullN(const std::vector<std::array<double, N>>& coord_a,
+                         const std::string& qhull_command = Qhull::COMMAND_DEFAULT) :
+        ConvexHullN(convert_vector_a(coord_a), qhull_command) {}
+
+    std::vector<size_t> list_vertices() const override { return qhull_.list_vertices(); }
+    std::vector<std::vector<size_t>> list_facets() const override { return qhull_.list_facets(); }
+
+    ConvexHull::facets_n_t facets_n() const override { return qhull_.facets_n(); }
+    std::vector<size_t> facets(size_t n) const override { return qhull_.facets(n); }
+
+private:
+    static coord_t convert_vector_v(const std::vector<std::vector<double>>& coord_v) {
+        coord_t coord;
+        coord.reserve(N * coord_v.size());
+        for (const auto& v : coord_v) {
+            ASSERT(N == v.size());
+            std::for_each_n(v.begin(), N, [&coord](auto x) { coord.emplace_back(x); });
+        }
+        return coord;
+    }
+
+    static coord_t convert_vector_a(const std::vector<std::array<double, N>>& coord_a) {
+        coord_t coord;
+        coord.reserve(N * coord_a.size());
+        for (const auto& a : coord_a) {
+            std::for_each_n(a.begin(), N, [&coord](auto x) { coord.emplace_back(x); });
+        }
+        return coord;
+    }
+
+    Qhull qhull_;
+};
+
+
+}  // namespace eckit::maths
diff --git a/src/eckit/maths/Qhull.cc b/src/eckit/maths/Qhull.cc
@@ -0,0 +1,144 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+
+#include "eckit/maths/Qhull.h"
+
+#include <set>
+#include <sstream>
+
+#include "eckit/log/Log.h"
+#include "eckit/maths/ConvexHull.h"
+
+#include "libqhullcpp/Qhull.h"
+#include "libqhullcpp/QhullFacetList.h"
+#include "libqhullcpp/QhullVertexSet.h"
+
+
+namespace eckit::maths {
+
+
+const char* Qhull::COMMAND_DEFAULT = "QJ";
+
+
+Qhull::Qhull(size_t N, const coord_t& coord, const std::string& command) {
+    ASSERT(0 < N && coord.size() % N == 0);
+
+    auto pointDimension = static_cast<int>(N);
+    auto pointCount     = static_cast<int>(coord.size() / N);
+
+    qh_ = new orgQhull::Qhull;
+    ASSERT(qh_ != nullptr);
+
+    std::ostringstream err;
+    qh_->setErrorStream(&err);
+    qh_->setOutputStream(&Log::info());
+    qh_->enableOutputStream();
+
+    try {
+        qh_->runQhull("", pointDimension, pointCount, coord.data(), command.c_str());
+        ASSERT(qh_->qhullStatus() == 0);
+    }
+    catch (const orgQhull::QhullError& e) {
+        static const std::set<int> DIMENSION_ERROR{6050};
+
+        static const std::set<int> INPUT_ERROR{6010, 6013, 6019, 6020, 6021, 6022, 6023, 6033, 6067, 6070,
+                                               6072, 6073, 6074, 6075, 6077, 6078, 6150, 6151, 6152, 6153,
+                                               6203, 6204, 6214, 6220, 6221, 6222, 6223, 6229, 6411, 6412};
+
+        static const std::set<int> OPTION_ERROR{6006, 6029, 6035, 6036, 6037, 6041, 6044, 6045, 6046,
+                                                6047, 6048, 6049, 6051, 6053, 6054, 6055, 6056, 6057,
+                                                6058, 6059, 6215, 6362, 6363, 6364, 6365, 6375};
+
+        static const std::set<int> PRECISION_ERROR{6012, 6109, 6110, 6111, 6112, 6113, 6114, 6115, 6116,
+                                                   6117, 6118, 6136, 6154, 6239, 6240, 6297, 6298, 6347,
+                                                   6348, 6354, 6379, 6380, 6417, 6418, 6422};
+
+        static const std::set<int> TOPOLOGY_ERROR{6001, 6107, 6155, 6168, 6170, 6208, 6227, 6260,
+                                                  6271, 6356, 6361, 6381, 6391, 6420, 6425};
+
+        auto is = [](int err, const auto& set) { return set.find(err) != set.end(); };
+
+        is(e.errorCode(), DIMENSION_ERROR)   ? throw ConvexHull::DimensionError(err.str(), e.errorCode(), command)
+        : is(e.errorCode(), INPUT_ERROR)     ? throw ConvexHull::InputError(err.str(), e.errorCode(), command)
+        : is(e.errorCode(), OPTION_ERROR)    ? throw ConvexHull::OptionError(err.str(), e.errorCode(), command)
+        : is(e.errorCode(), PRECISION_ERROR) ? throw ConvexHull::PrecisionError(err.str(), e.errorCode(), command)
+        : is(e.errorCode(), TOPOLOGY_ERROR)  ? throw ConvexHull::TopologyError(err.str(), e.errorCode(), command)
+                                             : throw ConvexHull::Exception(err.str(), e.errorCode(), command);
+    }
+}
+
+
+Qhull::~Qhull() {
+    delete qh_;
+}
+
+
+std::vector<size_t> Qhull::list_vertices() const {
+    std::vector<size_t> vertices;
+    vertices.reserve(qh_->vertexCount());
+
+    for (const auto& vertex : qh_->vertexList()) {
+        vertices.emplace_back(vertex.point().id());
+    }
+
+    return vertices;
+}
+
+
+std::vector<std::vector<size_t>> Qhull::list_facets() const {
+    std::vector<std::vector<size_t>> facets;
+    facets.reserve(qh_->facetCount());
+
+    for (const auto& facet : qh_->facetList()) {
+        const auto vertices = facet.vertices();
+
+        std::vector<size_t> f;
+        f.reserve(vertices.size());
+
+        for (const auto& vertex : vertices) {
+            f.emplace_back(vertex.point().id());
+        }
+
+        facets.emplace_back(std::move(f));
+    }
+
+    return facets;
+}
+
+
+Qhull::facets_n_t Qhull::facets_n() const {
+    facets_n_t counts;
+    for (const auto& facet : qh_->facetList()) {
+        counts[facet.vertices().size()]++;
+    }
+    return counts;
+}
+
+
+std::vector<size_t> Qhull::facets(size_t n) const {
+    ASSERT(0 < n);
+
+    std::vector<size_t> facets;
+    facets.reserve(qh_->facetCount() * n);
+
+    for (const auto& facet : qh_->facetList()) {
+        if (const auto vertices = facet.vertices(); vertices.size() == n) {
+            for (const auto& vertex : vertices) {
+                facets.emplace_back(vertex.point().id());
+            }
+        }
+    }
+
+    return facets;
+}
+
+
+}  // namespace eckit::maths