From dfd12e4d0de8b9f999fa7e48a67a6aada4edcce5 Mon Sep 17 00:00:00 2001 From: Oscar Torreno Date: Tue, 12 Mar 2019 16:09:23 +0100 Subject: [PATCH] Fix kShape for double array (#91) --- src/khiva/clustering.cpp | 77 +++++++++++++++++++++++++--------------- test/clusteringTest.cpp | 34 ++++++++++++++++-- 2 files changed, 80 insertions(+), 31 deletions(-) diff --git a/src/khiva/clustering.cpp b/src/khiva/clustering.cpp index 7789fd4..4ff0454 100644 --- a/src/khiva/clustering.cpp +++ b/src/khiva/clustering.cpp @@ -18,7 +18,7 @@ * @param k The number of centroids. * @return The new centroids. */ -af::array calculateInitialMeans(af::array tss, int k) { return af::constant(0.0f, tss.dims(0), k); } +af::array calculateInitialMeans(af::array tss, int k) { return af::constant(0, tss.dims(0), k, tss.type()); } /** * Computes The euclidean distance for a tiled time series agains k-means. @@ -40,13 +40,13 @@ af::array kEuclideanDistance(af::array tts, af::array means) { * @param labels The ids of the closes mean for all time series. */ void euclideanDistance(af::array tss, af::array means, af::array &minDistance, af::array &idxs) { - int nseries = tss.dims(1); - af::array kDistances = af::constant(0.0, means.dims(1), nseries); + int nSeries = tss.dims(1); + af::array kDistances = af::constant(0.0, means.dims(1), nSeries, tss.type()); // This for loop could be parallel, not parallelized to keep memory footprint low - for (int i = 0; i < nseries; i++) { - af::array tiledSerie = af::tile(tss.col(i), 1, means.dims(1)); - kDistances(af::span, i) = kEuclideanDistance(tiledSerie, means); + for (int i = 0; i < nSeries; i++) { + af::array tiledSeries = af::tile(tss.col(i), 1, means.dims(1)); + kDistances(af::span, i) = kEuclideanDistance(tiledSeries, means); } af::min(minDistance, idxs, kDistances, 0); } @@ -61,7 +61,7 @@ void euclideanDistance(af::array tss, af::array means, af::array &minDistance, a */ af::array computeNewMeans(af::array tss, af::array labels, int k) { af::array labelsTiled = af::tile(labels, tss.dims(0)); - af::array newMeans = af::constant(0.0, tss.dims(0), k); + af::array newMeans = af::constant(0.0, tss.dims(0), k, tss.type()); gfor(af::seq ii, k) { newMeans(af::span, ii) = af::sum(tss * (labelsTiled == ii), 1) / (af::sum(af::sum((labels == ii), 3), 1)); @@ -133,7 +133,7 @@ void khiva::clustering::kMeans(af::array tss, int k, af::array ¢roids, af::a labels = generateRandomLabels(tss.dims(1), k); } - af::array distances = af::constant(0.0f, tss.dims(1)); + af::array distances = af::constant(0, tss.dims(1), tss.type()); af::array newMeans; int iter = 0; @@ -204,18 +204,36 @@ af::array eigenValues(af::array matrix) { * @return The first Eigen vector. */ af::array getFirstEigenVector(af::array m) { - float *matHost = m.host(); - Eigen::MatrixXf mat = Eigen::Map(matHost, m.dims(0), m.dims(1)); - - // Compute Eigen Values. - Eigen::VectorXcf eivals = mat.eigenvalues(); - Eigen::VectorXf reEIVals = eivals.real(); - af::array eigenValues = af::array(m.dims(0), reEIVals.data()); - - // Compute Eigen Vectors. - Eigen::EigenSolver solution(mat); - Eigen::MatrixXf reEIVectors = solution.eigenvectors().real(); - af::array eigenVectors = af::array(m.dims(0), m.dims(1), reEIVectors.data()); + af::array eigenValues; + af::array eigenVectors; + + if (m.type() == af::dtype::f64) { + double *matHost = m.host(); + Eigen::MatrixXd mat = Eigen::Map(matHost, m.dims(0), m.dims(1)); + + // Compute Eigen Values. + Eigen::VectorXcd eivals = mat.eigenvalues(); + Eigen::VectorXd reEIVals = eivals.real(); + eigenValues = af::array(m.dims(0), reEIVals.data()); + + // Compute Eigen Vectors. + Eigen::EigenSolver solution(mat); + Eigen::MatrixXd reEIVectors = solution.eigenvectors().real(); + eigenVectors = af::array(m.dims(0), m.dims(1), reEIVectors.data()); + } else if (m.type() == af::dtype::f32) { + float *matHost = m.host(); + Eigen::MatrixXf mat = Eigen::Map(matHost, m.dims(0), m.dims(1)); + + // Compute Eigen Values. + Eigen::VectorXcf eivals = mat.eigenvalues(); + Eigen::VectorXf reEIVals = eivals.real(); + eigenValues = af::array(m.dims(0), reEIVals.data()); + + // Compute Eigen Vectors. + Eigen::EigenSolver solution(mat); + Eigen::MatrixXf reEIVectors = solution.eigenvectors().real(); + eigenVectors = af::array(m.dims(0), m.dims(1), reEIVectors.data()); + } // Get maximum Eigen Value af::array maxEigenValue; @@ -256,7 +274,7 @@ af::array ncc3Dim(af::array tss, af::array centroids) { int distanceSize = static_cast(centroids.dims(0)) * 2 - 1; af::array cc = af::constant(0, static_cast(centroids.dims(1)), static_cast(tss.dims(1)), - distanceSize); + distanceSize, tss.type()); for (unsigned int i = 0; i < static_cast(centroids.dims(1)); i++) { for (unsigned int j = 0; j < static_cast(tss.dims(1)); j++) { cc(i, j, af::span) = af::convolve(tss.col(j), af::flip(centroids.col(i), 0), AF_CONV_EXPAND); @@ -286,9 +304,9 @@ af::array SBDShifted(af::array ts, af::array centroid) { int shift = index - tsLength + 1; if (shift >= 0) { - shiftedTS = af::join(0, af::constant(0, shift), ts(af::range(tsLength - shift), 0)); + shiftedTS = af::join(0, af::constant(0, shift, ts.type()), ts(af::range(tsLength - shift), 0)); } else { - shiftedTS = af::join(0, ts(af::range(tsLength + shift) - shift, 0), af::constant(0, -shift)); + shiftedTS = af::join(0, ts(af::range(tsLength + shift) - shift, 0), af::constant(0, -shift, ts.type())); } return shiftedTS; @@ -304,8 +322,8 @@ af::array SBDShifted(af::array ts, af::array centroid) { af::array shapeExtraction(af::array tss, af::array centroid) { int ntss = tss.dims(1); int nelements = tss.dims(0); - af::array shiftedTSS = af::constant(0.0f, tss.dims(0), tss.dims(1)); - af::array shiftedTSi = af::constant(0.0f, tss.dims(0), 1); + af::array shiftedTSS = af::constant(0, tss.dims(0), tss.dims(1), tss.type()); + af::array shiftedTSi = af::constant(0, tss.dims(0), 1, tss.type()); af::array dist; for (int i = 0; i < ntss; i++) { @@ -316,7 +334,8 @@ af::array shapeExtraction(af::array tss, af::array centroid) { shiftedTSS = khiva::normalization::znorm(shiftedTSS); af::array s = af::matmul(shiftedTSS, shiftedTSS.T()); - af::array q = af::identity(nelements, nelements) - (af::constant(1.0, nelements, nelements) / nelements); + af::array q = + af::identity(nelements, nelements) - (af::constant(1.0, nelements, nelements, tss.type()) / nelements); af::array m = af::matmul(q.T(), s, q); af::array first = getFirstEigenVector(m); @@ -365,7 +384,7 @@ af::array refinementStep(af::array tss, af::array centroids, af::array labels) { * @return The new set of labels. */ af::array assignmentStep(af::array tss, af::array centroids, af::array labels) { - af::array min = af::constant(std::numeric_limits::max(), tss.dims(1)); + af::array min = af::constant(std::numeric_limits::max(), tss.dims(1), tss.type()); af::array distances = 1 - af::max(ncc3Dim(tss, centroids), 2); af::min(min, labels, distances, 0); labels = labels.T(); @@ -379,7 +398,7 @@ void khiva::clustering::kShape(af::array tss, int k, af::array ¢roids, af::a unsigned int nElements = static_cast(tss.dims(0)); if (centroids.isempty()) { - centroids = af::constant(0.0f, nElements, k); + centroids = af::constant(0, nElements, k, tss.type()); } if (labels.isempty()) { @@ -387,7 +406,7 @@ void khiva::clustering::kShape(af::array tss, int k, af::array ¢roids, af::a } af::array normTSS = khiva::normalization::znorm(tss); - af::array distances = af::constant(0.0f, nTimeseries, k); + af::array distances = af::constant(0, nTimeseries, k, tss.type()); af::array newCentroids; float error = std::numeric_limits::max(); diff --git a/test/clusteringTest.cpp b/test/clusteringTest.cpp index c2858b8..9a200af 100644 --- a/test/clusteringTest.cpp +++ b/test/clusteringTest.cpp @@ -54,7 +54,7 @@ void kmeans2() { } } -void kShape() { +void kShapeFloat() { float tolerance = 1e-10; int maxIter = 100; float a[35] = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 0.0f, 10.0f, 4.0f, 5.0f, 7.0f, @@ -84,6 +84,36 @@ void kShape() { } } +void kShapeDouble() { + double tolerance = 1e-10; + int maxIter = 100; + double a[35] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 0.0, 10.0, 4.0, 5.0, 7.0, -3.0, 0.0, -1.0, 15.0, -12.0, 8.0, + 9.0, 4.0, 5.0, 2.0, 8.0, 7.0, -6.0, -1.0, 2.0, 9.0, -5.0, -5.0, -6.0, 7.0, 9.0, 9.0, 0.0}; + + unsigned int idxh[5] = {0, 1, 2, 0, 1}; + + double expected_c[21] = {-0.5234, 0.1560, -0.3627, -1.2764, -0.7781, 0.9135, 1.8711, + -0.7825, 1.5990, 0.1701, 0.4082, 0.8845, -1.4969, -0.7825, + -0.6278, 1.3812, -2.0090, 0.5022, 0.6278, -0.0000, 0.1256}; + int k = 3; + int nElements = 7; + int ntss = 5; + + af::array data = af::array(nElements, ntss, a); + af::array idx = af::array(ntss, 1, idxh); + af::array centroids = af::constant(0, nElements, k, data.type()); + + khiva::clustering::kShape(data, k, centroids, idx, tolerance, maxIter); + + double *calculated_c = centroids.host(); + unsigned int *calculated_l = idx.host(); + + for (size_t i = 0; i < 21; i++) { + ASSERT_NEAR(calculated_c[i], expected_c[i], 1e-3); + } +} + KHIVA_TEST(ClusteringTests, KMeans, kmeans) KHIVA_TEST(ClusteringTests, KMeans2, kmeans2) -KHIVA_TEST(ClusteringTests, KShape, kShape) +KHIVA_TEST(ClusteringTests, KShapeFloat, kShapeFloat) +KHIVA_TEST(ClusteringTests, KShapeDouble, kShapeDouble)