diff --git a/include/khiva/dimensionality.h b/include/khiva/dimensionality.h index 3054b91..5102822 100644 --- a/include/khiva/dimensionality.h +++ b/include/khiva/dimensionality.h @@ -16,6 +16,9 @@ namespace khiva { namespace dimensionality { +template +using TPoint = std::pair; + using Point = std::pair; using Segment = std::pair; @@ -191,8 +194,9 @@ KHIVAAPI af::array SAX(const af::array &a, int alphabetSize); * * @return std:vector where the number of points has been reduced to numPoints. */ -KHIVAAPI std::vector visvalingam(const std::vector &pointList, int64_t numPoints, - int64_t scale = 1000000000); +template +KHIVAAPI std::vector> visvalingam(const std::vector> &pointList, int64_t numPoints, + int64_t scale = 1000000000); /** * @brief Reduces a set of points by applying the Visvalingam method (minimum triangle area) until the number diff --git a/src/khiva/dimensionality.cpp b/src/khiva/dimensionality.cpp index 80e632a..da38285 100644 --- a/src/khiva/dimensionality.cpp +++ b/src/khiva/dimensionality.cpp @@ -23,9 +23,10 @@ using namespace khiva::dimensionality; namespace { +template struct VisvalingamSummaryPoint { - float x; - float y; + T x; + T y; int64_t area; }; @@ -187,12 +188,13 @@ float calculateError(const std::vector &ts, int start, int end) { Segment merge(Segment s1, Segment s2) { return {s1.first, s2.second}; } -int64_t computeTriangleArea(const VisvalingamSummaryPoint &a, const VisvalingamSummaryPoint &b, - const VisvalingamSummaryPoint &c, const long scale = 1e9) { - float f1 = a.x * (b.y - c.y); - float f2 = b.x * (c.y - a.y); - float f3 = c.x * (a.y - b.y); - return static_cast(std::abs((static_cast(f1) + f2 + f3) / 2.0f) * scale); +template +int64_t computeTriangleArea(const VisvalingamSummaryPoint &a, const VisvalingamSummaryPoint &b, + const VisvalingamSummaryPoint &c, const long scale = 1e9) { + auto f1 = a.x * (b.y - c.y); + auto f2 = b.x * (c.y - a.y); + auto f3 = c.x * (a.y - b.y); + return static_cast(std::abs((static_cast(f1) + f2 + f3) / static_cast(2.0)) * scale); } template @@ -208,9 +210,10 @@ class mapComparator { } }; -void recomputeAreaNeighbor(std::map::iterator &iterator_point, +template +void recomputeAreaNeighbor(typename std::map>::iterator &iterator_point, std::set, mapComparator> &point_indexer, - std::map &points, const int64_t scale) { + std::map> &points, const int64_t scale) { auto im1m1 = shiftIterator(iterator_point, -1); auto im1p1 = shiftIterator(iterator_point, 1); auto original_position_minus1 = iterator_point->first; @@ -218,7 +221,7 @@ void recomputeAreaNeighbor(std::map::iterator auto old_area_minus1 = iterator_point->second.area; auto new_area_minus1 = computeTriangleArea(im1m1->second, iterator_point->second, im1p1->second, scale); points[iterator_point->first] = - VisvalingamSummaryPoint{iterator_point->second.x, iterator_point->second.y, new_area_minus1}; + VisvalingamSummaryPoint{iterator_point->second.x, iterator_point->second.y, new_area_minus1}; auto it = point_indexer.find(std::make_pair(old_area_minus1, original_position_minus1)); point_indexer.erase(it); @@ -229,7 +232,7 @@ void recomputeAreaNeighbor(std::map::iterator } // namespace std::vector khiva::dimensionality::PAA(const std::vector &points, int bins) { - float xrange = points.back().first - points.front().first; + float xrange = points.back().first - points.front().first; float width_bin = xrange / bins; float reduction = bins / xrange; @@ -237,11 +240,11 @@ std::vector khiva::dimensionality::PAA(const std::vector &points, std::vector counter(bins, 0); // Iterating over the time series - for (const auto &p : points) { + for (const auto &p : points) { auto pos = static_cast(std::min(p.first * reduction, (float)(bins - 1))); sum[pos] += p.second; counter[pos] += 1; - } + } std::vector result; result.reserve(bins); @@ -609,16 +612,17 @@ af::array khiva::dimensionality::SAX(const af::array &a, int alphabet_size) { return result; } -std::vector khiva::dimensionality::visvalingam(const std::vector &pointList, int64_t numPoints, - int64_t scale) { - std::map points; +template +std::vector> khiva::dimensionality::visvalingam(const std::vector> &pointList, int64_t numPoints, + int64_t scale) { + std::map> points; std::set, mapComparator> point_indexer; int64_t counter = 0; std::transform( pointList.cbegin(), pointList.cend(), std::inserter(points, points.end()), [&counter](const Point &point) { return std::make_pair( - counter++, VisvalingamSummaryPoint{point.first, point.second, std::numeric_limits::max()}); + counter++, VisvalingamSummaryPoint{point.first, point.second, std::numeric_limits::max()}); }); auto points_to_be_deleted = pointList.size() - numPoints; @@ -652,34 +656,34 @@ std::vector khiva::dimensionality::visvalingam(const std::vector & } } - std::vector out_vector; + std::vector> out_vector; out_vector.reserve(numPoints); std::transform(points.cbegin(), points.cend(), std::back_inserter(out_vector), - [](const std::pair &p) { + [](const std::pair> &p) { return Point{p.second.x, p.second.y}; }); return out_vector; } -af::array khiva::dimensionality::visvalingam(const af::array &pointList, int numPoints) { - if (pointList.dims(1) != 2) { - throw std::invalid_argument("Invalid dims. Khiva array with two columns expected (x axis and y axis)."); - } - std::vector points; + +template +af::array visvalingamBridge(const af::array &pointList, int numPoints) { + std::vector> points; points.reserve(pointList.dims(0)); - auto x = khiva::utils::makeScopedHostPtr(pointList.col(0).host()); - auto y = khiva::utils::makeScopedHostPtr(pointList.col(1).host()); + + auto x = khiva::utils::makeScopedHostPtr(pointList.col(0).host()); + auto y = khiva::utils::makeScopedHostPtr(pointList.col(1).host()); for (int i = 0; i < pointList.dims(0); i++) { points.emplace_back(x[i], y[i]); } - std::vector rPoints = visvalingam(points, numPoints); - af::array out = af::constant(0, rPoints.size(), 2); + auto rPoints = visvalingam(points, numPoints); - std::vector vx; + std::vector vx; vx.reserve(rPoints.size()); - std::vector vy; + + std::vector vy; vy.reserve(rPoints.size()); for (const auto &rPoint : rPoints) { @@ -692,3 +696,17 @@ af::array khiva::dimensionality::visvalingam(const af::array &pointList, int num return af::join(1, ox, oy); } + +af::array khiva::dimensionality::visvalingam(const af::array &pointList, int numPoints) { + if (pointList.dims(1) != 2) { + throw std::invalid_argument("Invalid dims. Khiva array with two columns expected (x axis and y axis)."); + } + + if (pointList.type() == af::dtype::f64) { + return visvalingamBridge(pointList, numPoints); + } else if (pointList.type() == af::dtype::f32) { + return visvalingamBridge(pointList, numPoints); + } else { + throw std::invalid_argument("Visvalingam is only supported for f32 and f64 types"); + } +} diff --git a/test/dimensionalityTest.cpp b/test/dimensionalityTest.cpp index 60262ce..52cac17 100644 --- a/test/dimensionalityTest.cpp +++ b/test/dimensionalityTest.cpp @@ -42,6 +42,8 @@ void paaNonDivisibleFloat() { } void paaNonDivisibleDouble() { + if (!af::isDoubleAvailable(af::getDevice())) return; + double pointList[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 0.0, 0.1, -0.1, 5.0, 6.0, 7.0, 8.1, 9.0, 9.0, 9.0}; af::array a(10, 2, pointList); @@ -68,7 +70,7 @@ void paaNorm() { std::vector expected = {{0.75, 0.05}, {2.25, -0.1}, {3.75, 5.5}, {5.25, 7.0}, {6.75, 8.55}, {8.25, 9.0}}; - ASSERT_EQ(expected, out); + ASSERT_EQ(expected, out); } void pip() { @@ -253,34 +255,49 @@ void saxException() { ASSERT_THROW(khiva::dimensionality::SAX(a, 6), std::invalid_argument); } -void visvalingam() { - std::vector pointList = {{0.0, 0.0}, {1.0, 0.1}, {2.0, -0.1}, {3.0, 5.0}, {4.0, 6.0}, - {5.0, 7.0}, {6.0, 8.1}, {7.0, 9.0}, {8.0, 9.0}, {9.0, 9.0}}; - std::vector expected = {{0.0, 0.0}, {2.0, -0.1}, {3.0, 5.0}, {7.0, 9.0}, {9.0, 9.0}}; +template +void visvalingam_templated() { + T pointList[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, + 0.0, 0.1, -0.1, 5.0, 6.0, 7.0, 8.1, 9.0, 9.0, 9.0}; + af::array points(10, 2, pointList); - auto out = khiva::dimensionality::visvalingam(pointList, 5); + std::vector> expected = { + {0.0, 0.0}, {2.0, -0.1}, {3.0, 5.0}, {7.0, 9.0}, {9.0, 9.0}}; - ASSERT_EQ(expected, out); -} + af::array res = khiva::dimensionality::visvalingam(points, 5); -void visvalingam2() { - float pointList[] = {0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f, - 0.0f, 0.1f, -0.1f, 5.0f, 6.0f, 7.0f, 8.1f, 9.0f, 9.0f, 9.0f}; - af::array points(10, 2, pointList); - std::vector expected = {{0.0, 0.0}, {2.0, -0.1}, {3.0, 5.0}, {7.0, 9.0}, {9.0, 9.0}}; + auto points_x = khiva::utils::makeScopedHostPtr(res.col(0).host()); + auto points_y = khiva::utils::makeScopedHostPtr(res.col(1).host()); - af::array res = khiva::dimensionality::visvalingam(points, 5); - auto points_x = khiva::utils::makeScopedHostPtr(res.col(0).host()); - auto points_y = khiva::utils::makeScopedHostPtr(res.col(1).host()); + auto expectedX = std::vector{0.0f, 2.0f, 3.0f, 7.0f, 9.0f}; + auto expectedY = std::vector{0.0f, -0.1f, 5.0f, 9.0f, 9.0f}; - auto expectedX = std::vector{0.0f, 2.0f, 3.0f, 7.0f, 9.0f}; - auto expectedY = std::vector{0.0f, -0.1f, 5.0f, 9.0f, 9.0f}; + auto pxVector = std::vector(points_x.get(), points_x.get() + res.col(0).elements()); + auto pyVector = std::vector(points_y.get(), points_y.get() + res.col(1).elements()); - auto pxVector = std::vector(points_x.get(), points_x.get() + res.col(0).elements()); - auto pyVector = std::vector(points_y.get(), points_y.get() + res.col(1).elements()); + std::vector xsbs; + std::vector ysbs; + std::transform(expectedX.begin(), expectedX.end(), pxVector.begin(), std::back_inserter(xsbs), + [](const T &e, const T &a) { return std::abs(e - a); }); + std::transform(expectedY.begin(), expectedY.end(), pyVector.begin(), std::back_inserter(ysbs), + [](const T &e, const T &a) { return std::abs(e - a); }); - ASSERT_EQ(expectedX, pxVector); - ASSERT_EQ(expectedY, pyVector); + for (auto diff : xsbs) { + ASSERT_LT(diff, 1e-08); + } + + for (auto diff : ysbs) { + ASSERT_LT(diff, 1e-08); + } +} + +void visvalingam_af() { + SCOPED_TRACE("f32 types..."); + visvalingam_templated(); + if (af::isDoubleAvailable(af::getDevice())) { + SCOPED_TRACE("f64 types..."); + visvalingam_templated(); + } } void visvalingamException() { @@ -307,6 +324,5 @@ KHIVA_TEST(DimensionalityTests, RamerDouglasPeuckerException, ramerDouglasPeucke KHIVA_TEST(DimensionalityTests, SAX, sax) KHIVA_TEST(DimensionalityTests, SAX2, sax2) KHIVA_TEST(DimensionalityTests, SAXException, saxException) -KHIVA_TEST(DimensionalityTests, Visvalingam, visvalingam) -KHIVA_TEST(DimensionalityTests, Visvalingam2, visvalingam2) +KHIVA_TEST(DimensionalityTests, VisvalingamAF, visvalingam_af) KHIVA_TEST(DimensionalityTests, VisvalingamException, visvalingamException)