Skip to content

Visvalingam runs on both 64 and 32 floats. #169

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Visvalingam runs on both 64 and 32 floats. #169

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3acacec
Visvalingam runs on both 64 and 32 floats.
justoruiz Sep 2, 2021
e3c2f9d
try to resolve problem with narrowing types on test.
justoruiz Sep 3, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 6 additions & 2 deletions include/khiva/dimensionality.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,9 @@ namespace khiva {

namespace dimensionality {

template <typename T>
using TPoint = std::pair<T, T>;

using Point = std::pair<float, float>;

using Segment = std::pair<int, int>;
Expand Down Expand Up @@ -191,8 +194,9 @@ KHIVAAPI af::array SAX(const af::array &a, int alphabetSize);
*
* @return std:vector<khiva::dimensionality::Point> where the number of points has been reduced to numPoints.
*/
KHIVAAPI std::vector<Point> visvalingam(const std::vector<Point> &pointList, int64_t numPoints,
int64_t scale = 1000000000);
template <typename T>
KHIVAAPI std::vector<TPoint<T>> visvalingam(const std::vector<TPoint<T>> &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
Expand Down
80 changes: 49 additions & 31 deletions src/khiva/dimensionality.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,10 @@ using namespace khiva::dimensionality;

namespace {

template <typename T>
struct VisvalingamSummaryPoint {
float x;
float y;
T x;
T y;
int64_t area;
};

Expand Down Expand Up @@ -187,12 +188,13 @@ float calculateError(const std::vector<Point> &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<int64_t>(std::abs((static_cast<double>(f1) + f2 + f3) / 2.0f) * scale);
template <typename T>
int64_t computeTriangleArea(const VisvalingamSummaryPoint<T> &a, const VisvalingamSummaryPoint<T> &b,
const VisvalingamSummaryPoint<T> &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<int64_t>(std::abs((static_cast<T>(f1) + f2 + f3) / static_cast<T>(2.0)) * scale);
}

template <typename Iter, typename Distance>
Expand All @@ -208,17 +210,18 @@ class mapComparator {
}
};

void recomputeAreaNeighbor(std::map<int64_t, VisvalingamSummaryPoint>::iterator &iterator_point,
template <typename T>
void recomputeAreaNeighbor(typename std::map<int64_t, VisvalingamSummaryPoint<T>>::iterator &iterator_point,
std::set<std::pair<int64_t, int64_t>, mapComparator> &point_indexer,
std::map<int64_t, VisvalingamSummaryPoint> &points, const int64_t scale) {
std::map<int64_t, VisvalingamSummaryPoint<T>> &points, const int64_t scale) {
auto im1m1 = shiftIterator(iterator_point, -1);
auto im1p1 = shiftIterator(iterator_point, 1);
auto original_position_minus1 = iterator_point->first;

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<T>{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);
Expand All @@ -229,19 +232,19 @@ void recomputeAreaNeighbor(std::map<int64_t, VisvalingamSummaryPoint>::iterator
} // namespace

std::vector<Point> khiva::dimensionality::PAA(const std::vector<Point> &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;

std::vector<float> sum(bins, 0.0);
std::vector<int> counter(bins, 0);

// Iterating over the time series
for (const auto &p : points) {
for (const auto &p : points) {
auto pos = static_cast<size_t>(std::min(p.first * reduction, (float)(bins - 1)));
sum[pos] += p.second;
counter[pos] += 1;
}
}

std::vector<Point> result;
result.reserve(bins);
Expand Down Expand Up @@ -609,16 +612,17 @@ af::array khiva::dimensionality::SAX(const af::array &a, int alphabet_size) {
return result;
}

std::vector<Point> khiva::dimensionality::visvalingam(const std::vector<Point> &pointList, int64_t numPoints,
int64_t scale) {
std::map<int64_t, VisvalingamSummaryPoint> points;
template <typename T>
std::vector<TPoint<T>> khiva::dimensionality::visvalingam(const std::vector<TPoint<T>> &pointList, int64_t numPoints,
int64_t scale) {
std::map<int64_t, VisvalingamSummaryPoint<T>> points;
std::set<std::pair<int64_t, int64_t>, 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<int64_t>::max()});
counter++, VisvalingamSummaryPoint<T>{point.first, point.second, std::numeric_limits<int64_t>::max()});
});

auto points_to_be_deleted = pointList.size() - numPoints;
Expand Down Expand Up @@ -652,34 +656,34 @@ std::vector<Point> khiva::dimensionality::visvalingam(const std::vector<Point> &
}
}

std::vector<Point> out_vector;
std::vector<TPoint<T>> out_vector;
out_vector.reserve(numPoints);
std::transform(points.cbegin(), points.cend(), std::back_inserter(out_vector),
[](const std::pair<int64_t, VisvalingamSummaryPoint> &p) {
[](const std::pair<int64_t, VisvalingamSummaryPoint<T>> &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<Point> points;

template <typename T>
af::array visvalingamBridge(const af::array &pointList, int numPoints) {
std::vector<TPoint<T>> points;
points.reserve(pointList.dims(0));
auto x = khiva::utils::makeScopedHostPtr(pointList.col(0).host<float>());
auto y = khiva::utils::makeScopedHostPtr(pointList.col(1).host<float>());

auto x = khiva::utils::makeScopedHostPtr(pointList.col(0).host<T>());
auto y = khiva::utils::makeScopedHostPtr(pointList.col(1).host<T>());

for (int i = 0; i < pointList.dims(0); i++) {
points.emplace_back(x[i], y[i]);
}

std::vector<Point> rPoints = visvalingam(points, numPoints);
af::array out = af::constant(0, rPoints.size(), 2);
auto rPoints = visvalingam(points, numPoints);

std::vector<float> vx;
std::vector<T> vx;
vx.reserve(rPoints.size());
std::vector<float> vy;

std::vector<T> vy;
vy.reserve(rPoints.size());

for (const auto &rPoint : rPoints) {
Expand All @@ -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<double>(pointList, numPoints);
} else if (pointList.type() == af::dtype::f32) {
return visvalingamBridge<float>(pointList, numPoints);
} else {
throw std::invalid_argument("Visvalingam is only supported for f32 and f64 types");
}
}
64 changes: 40 additions & 24 deletions test/dimensionalityTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -68,7 +70,7 @@ void paaNorm() {
std::vector<khiva::dimensionality::Point> 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() {
Expand Down Expand Up @@ -253,34 +255,49 @@ void saxException() {
ASSERT_THROW(khiva::dimensionality::SAX(a, 6), std::invalid_argument);
}

void visvalingam() {
std::vector<khiva::dimensionality::Point> 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<khiva::dimensionality::Point> expected = {{0.0, 0.0}, {2.0, -0.1}, {3.0, 5.0}, {7.0, 9.0}, {9.0, 9.0}};
template <typename T>
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<khiva::dimensionality::TPoint<T>> 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<khiva::dimensionality::Point> 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<T>());
auto points_y = khiva::utils::makeScopedHostPtr(res.col(1).host<T>());

af::array res = khiva::dimensionality::visvalingam(points, 5);
auto points_x = khiva::utils::makeScopedHostPtr(res.col(0).host<float>());
auto points_y = khiva::utils::makeScopedHostPtr(res.col(1).host<float>());
auto expectedX = std::vector<T>{0.0, 2.0, 3.0, 7.0, 9.0};
auto expectedY = std::vector<T>{0.0, -0.1, 5.0, 9.0, 9.0};

auto expectedX = std::vector<float>{0.0f, 2.0f, 3.0f, 7.0f, 9.0f};
auto expectedY = std::vector<float>{0.0f, -0.1f, 5.0f, 9.0f, 9.0f};
auto pxVector = std::vector<T>(points_x.get(), points_x.get() + res.col(0).elements());
auto pyVector = std::vector<T>(points_y.get(), points_y.get() + res.col(1).elements());

auto pxVector = std::vector<float>(points_x.get(), points_x.get() + res.col(0).elements());
auto pyVector = std::vector<float>(points_y.get(), points_y.get() + res.col(1).elements());
std::vector<T> xsbs;
std::vector<T> 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<float>();
if (af::isDoubleAvailable(af::getDevice())) {
SCOPED_TRACE("f64 types...");
visvalingam_templated<double>();
}
}

void visvalingamException() {
Expand All @@ -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)