eckit::geo::range

src/eckit/geo/Range.h

@@ -58,7 +58,8 @@ class Range : protected std::vector<double> {
 
     static std::string className() { return "range"; }
 
-    bool empty() const { return vector::empty(); }
+    using vector::empty;
+
     size_t size() const { return empty() ? n_ : vector::size(); }
 
     virtual const std::vector<double>& values() const = 0;

src/eckit/geo/range/Gaussian.cc

@@ -22,59 +22,35 @@
 namespace eckit::geo::range {
 
 
-Gaussian::Gaussian(size_t N) :
-    Gaussian(N, 90., -90) {
-}
-
-
 Gaussian::Gaussian(size_t N, double a, double b, double precision) :
     Range(2 * N),
     N_(N),
     a_(a),
     b_(b),
-    precision_(precision) {
+    eps_(precision) {
     ASSERT(N > 0);
-    ASSERT(precision_ >= 0.);
+    ASSERT(eps_ >= 0.);
+
+    // pre-calculate on cropping
+    auto [min, max] = std::minmax(a_, b_);
+    if (!types::is_approximately_equal(min, -90., eps_) || !types::is_approximately_equal(max, 90., eps_)) {
+        auto [from, to] = util::monotonic_crop(values(), min, max, precision);
+        if (to != end()) {
+            erase(to);
+        }
+        if (from != begin()) {
+            erase(begin(), from);
+        }
+        ASSERT(!empty());
+    }
 }
 
 
 const std::vector<double>& Gaussian::values() const {
     if (empty()) {
-        auto& v = const_cast<std::vector<double>&>(valuesVector());
-        v       = util::gaussian_latitudes(N_, a_ < b_);
-        ASSERT(v.size() == 2 * N_);
-
-        const bool same(a_ == b_);
-
-        auto a = a_;
-        auto b = b_;
-
-        if (a < v.back()) {
-            a = v.back();
-        }
-        else {
-            auto best = std::lower_bound(v.begin(), v.end(), a, [&](double l1, double l2) {
-                return !types::is_approximately_equal(l1, l2, precision_) && !(l1 < l2);
-            });
-            ASSERT(best != v.end());
-            a = *best;
-        }
-
-        if (same) {
-            b = a;
-        }
-        else if (b > v.front()) {
-            b = v.front();
-        }
-        else {
-            auto best = std::lower_bound(v.rbegin(), v.rend(), b, [&](double l1, double l2) {
-                return !types::is_approximately_equal(l1, l2, precision_) && !(l1 > l2);
-            });
-            ASSERT(best != v.rend());
-            b = *best;
-        }
-
-        ASSERT((a_ <= b_) == (a <= b));
+        const_cast<std::vector<double>&>(valuesVector()) = util::gaussian_latitudes(N_, a_ < b_);
+        ASSERT(!empty());
+        ASSERT(size() == 2 * N_);
     }
 
     return *this;

src/eckit/geo/range/Gaussian.h

@@ -28,7 +28,9 @@ class Gaussian final : public Range {
 
     // -- Constructors
 
-    explicit Gaussian(size_t N);
+    explicit Gaussian(size_t N, double precision = 0.) :
+        Gaussian(N, 90., -90., precision) {}
+
     Gaussian(size_t N, double a, double b, double precision = 0.);
 
     // -- Destructor
@@ -60,7 +62,7 @@ class Gaussian final : public Range {
     const size_t N_;
     const double a_;
     const double b_;
-    const double precision_;
+    const double eps_;
 
     // -- Methods
     // None

tests/geo/CMakeLists.txt

@@ -9,6 +9,7 @@ foreach( _test
         points
         polygon
         projection
+        range
         search
         sphere
         types

tests/geo/test_range.cc

@@ -0,0 +1,64 @@
+/*
+ * (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 <iostream>
+#include <vector>
+
+#include "eckit/geo/range/Gaussian.h"
+#include "eckit/geo/range/GlobalRegular.h"
+#include "eckit/geo/range/LocalRegular.h"
+#include "eckit/testing/Test.h"
+#include "eckit/types/FloatCompare.h"
+
+
+#define EXPECT_APPROX(a, b) EXPECT(::eckit::types::is_approximately_equal((a), (b), 1e-12))
+
+
+int main(int argc, const char* argv[]) {
+    using namespace eckit::geo::range;
+
+
+    {
+        std::vector<double> ref{59.44440828916676,
+                                19.87571914744090,
+                                -19.87571914744090,
+                                -59.44440828916676};
+
+        auto global = Gaussian(2);
+        EXPECT(global.size() == ref.size());
+
+        size_t i = 0;
+        for (const auto& test : global.values()) {
+            EXPECT_APPROX(test, ref[i++]);
+        }
+
+        auto cropped = Gaussian(2, 50., -50., 1e-3);
+        EXPECT(cropped.size() == ref.size() - 2);
+
+        i = 1;
+        for (const auto& test : cropped.values()) {
+            EXPECT_APPROX(test, ref[i++]);
+        }
+
+        EXPECT(Gaussian(2, 59.444, -59.444, 1e-3).size() == 4);
+        EXPECT(Gaussian(2, 59.444, -59.444, 1e-6).size() == 2);
+        EXPECT(Gaussian(2, 59.444, -59.445, 1e-6).size() == 3);
+
+        auto single = Gaussian(2, -59.444, -59.444, 1e-3);
+        EXPECT(single.size() == 1);
+
+        EXPECT_APPROX(single.values().front(), ref.back());
+    }
+
+
+    return 0;
+}

tests/geo/test_util.cc

@@ -27,46 +27,6 @@ std::ostream& operator<<(std::ostream& out, const std::vector<T>& v) {
 }
 
 
-template <typename X>
-struct iterator_t {
-    X& cnt;
-    size_t pos;
-
-    bool operator!=(const iterator_t& other) const { return &cnt != &other.cnt || pos != other.pos; }
-
-    iterator_t& operator++() {
-        pos++;
-        return *this;
-    }
-
-    iterator_t operator++(int) {
-        auto old = *this;
-        operator++();
-        return old;
-    }
-
-    typename X::value_type& operator*() { return cnt.at(pos); }
-
-    const typename X::value_type& operator*() const { return cnt.at(pos); }
-};
-
-
-class iterable_t : public std::vector<double> {
-    using iterator       = iterator_t<iterable_t>;
-    using const_iterator = iterator_t<const iterable_t>;
-
-public:
-    using vector::vector;
-
-    iterator begin() { return {*this, 0}; }
-    iterator end() { return {*this, this->size()}; }
-    const_iterator cbegin() const { return {*this, 0}; }
-    const_iterator cend() const { return {*this, this->size()}; }
-    const_iterator begin() const { return cbegin(); }
-    const_iterator end() const { return cend(); }
-};
-
-
 int main(int argc, char* argv[]) {
     using namespace eckit::geo::util;