Add correlation distance algorithm (#3131)

* Adds correlation distance algorithm, examples, docs, and unit tests

Signed-off-by: North Iii <[email protected]>
Co-authored-by: Victoriya Fedotova <[email protected]>

Author: richardnorth3

cpp/daal/src/algorithms/cordistance/cordistance_batch_impl.i

@@ -56,13 +56,25 @@ template <typename algorithmFPType, Method method, CpuType cpu>
 services::Status DistanceKernel<algorithmFPType, method, cpu>::compute(const size_t na, const NumericTable * const * a, const size_t nr,
                                                                        NumericTable * r[], const daal::algorithms::Parameter * par)
 {
-    NumericTable * xTable                          = const_cast<NumericTable *>(a[0]); /* Input data */
+    NumericTable * xTable                          = const_cast<NumericTable *>(a[0]); /* x Input data */
     NumericTable * rTable                          = const_cast<NumericTable *>(r[0]); /* Result */
     const NumericTableIface::StorageLayout rLayout = r[0]->getDataLayout();
 
     if (isFull<algorithmFPType, cpu>(rLayout))
     {
-        return corDistanceFull<algorithmFPType, cpu>(xTable, rTable);
+        if (na == 1)
+        {
+            return corDistanceFull<algorithmFPType, cpu>(xTable, rTable);
+        }
+        else if (na == 2)
+        {
+            NumericTable * yTable = const_cast<NumericTable *>(a[1]); /* y Input data */
+            return corDistanceFull<algorithmFPType, cpu>(xTable, yTable, rTable);
+        }
+        else
+        {
+            return services::Status(services::ErrorIncorrectNumberOfInputNumericTables);
+        }
     }
     else
     {

cpp/daal/src/algorithms/cordistance/cordistance_dense_default_batch_fpt_cpu.cpp

@@ -38,7 +38,7 @@ template class BatchContainer<DAAL_FPTYPE, defaultDense, DAAL_CPU>;
 }
 namespace internal
 {
-template class DistanceKernel<DAAL_FPTYPE, defaultDense, DAAL_CPU>;
+template class DAAL_EXPORT DistanceKernel<DAAL_FPTYPE, defaultDense, DAAL_CPU>;
 
 } // namespace internal
 

cpp/daal/src/algorithms/cordistance/cordistance_full_impl.i

@@ -22,6 +22,7 @@
 */
 #include "src/services/service_defines.h"
 using namespace daal::internal;
+using namespace daal::services;
 
 namespace daal
 {
@@ -306,6 +307,139 @@ services::Status corDistanceFull(const NumericTable * xTable, NumericTable * rTa
     return safeStat.detach();
 }
 
+template <typename algorithmFPType, CpuType cpu>
+services::Status corDistanceFull(const NumericTable * xTable, const NumericTable * yTable, NumericTable * rTable)
+{
+    size_t p         = xTable->getNumberOfColumns(); /* Dimension of input feature vector */
+    size_t nVectors1 = xTable->getNumberOfRows();    /* Number of input vectors in X */
+    size_t nVectors2 = yTable->getNumberOfRows();    /* Number of input vectors in Y */
+
+    size_t nBlocks1 = nVectors1 / blockSizeDefault;
+    nBlocks1 += (nBlocks1 * blockSizeDefault != nVectors1);
+
+    size_t nBlocks2 = nVectors2 / blockSizeDefault;
+    nBlocks2 += (nBlocks2 * blockSizeDefault != nVectors2);
+
+    SafeStatus safeStat;
+
+    /* Allocate yMean for all Y vectors before the loop */
+    TArray<algorithmFPType, cpu> yMeanArr(nVectors2);
+    algorithmFPType * yMean = yMeanArr.get();
+    DAAL_CHECK(yMean, ErrorMemoryAllocationFailed);
+
+    /* Compute means for all Y vectors before the loop */
+    for (size_t k2 = 0; k2 < nBlocks2; k2++)
+    {
+        DAAL_INT blockSize2 = blockSizeDefault;
+        if (k2 == nBlocks2 - 1)
+        {
+            blockSize2 = nVectors2 - k2 * blockSizeDefault;
+        }
+
+        size_t shift2 = k2 * blockSizeDefault;
+
+        /* read access to blockSize2 rows in input dataset Y */
+        ReadRows<algorithmFPType, cpu> yBlock(*const_cast<NumericTable *>(yTable), shift2, blockSize2);
+        DAAL_CHECK_BLOCK_STATUS(yBlock);
+        const algorithmFPType * y = yBlock.get();
+
+        for (size_t j = 0; j < blockSize2; j++)
+        {
+            yMean[shift2 + j] = 0.0;
+            for (size_t k = 0; k < p; k++)
+            {
+                yMean[shift2 + j] += y[j * p + k];
+            }
+            yMean[shift2 + j] /= p;
+        }
+    }
+
+    /* compute results for blocks of the distance matrix */
+    daal::threader_for(nBlocks1, nBlocks1, [=, &safeStat](size_t k1) {
+        DAAL_INT blockSize1 = blockSizeDefault;
+        if (k1 == nBlocks1 - 1)
+        {
+            blockSize1 = nVectors1 - k1 * blockSizeDefault;
+        }
+
+        /* read access to blockSize1 rows in input dataset X at k1*blockSizeDefault*p row */
+        ReadRows<algorithmFPType, cpu> xBlock(*const_cast<NumericTable *>(xTable), k1 * blockSizeDefault, blockSize1);
+        DAAL_CHECK_BLOCK_STATUS_THR(xBlock);
+        const algorithmFPType * x = xBlock.get();
+
+        /* write access to blockSize1 rows in output dataset */
+        WriteOnlyRows<algorithmFPType, cpu> rBlock(rTable, k1 * blockSizeDefault, blockSize1);
+        DAAL_CHECK_BLOCK_STATUS_THR(rBlock);
+        algorithmFPType * r = rBlock.get();
+
+        /* Compute means for rows in X block */
+        TArrayScalable<algorithmFPType, cpu> xMeanArr(blockSize1);
+        algorithmFPType * xMean = xMeanArr.get();
+        DAAL_CHECK_THR(xMean, ErrorMemoryAllocationFailed);
+
+        for (size_t i = 0; i < blockSize1; i++)
+        {
+            xMean[i] = 0.0;
+            for (size_t j = 0; j < p; j++)
+            {
+                xMean[i] += x[i * p + j];
+            }
+            xMean[i] /= p;
+        }
+
+        for (size_t k2 = 0; k2 < nBlocks2; k2++)
+        {
+            DAAL_INT blockSize2 = blockSizeDefault;
+            if (k2 == nBlocks2 - 1)
+            {
+                blockSize2 = nVectors2 - k2 * blockSizeDefault;
+            }
+
+            size_t shift2 = k2 * blockSizeDefault;
+
+            /* read access to blockSize2 rows in input dataset Y */
+            ReadRows<algorithmFPType, cpu> yBlock(*const_cast<NumericTable *>(yTable), shift2, blockSize2);
+            DAAL_CHECK_BLOCK_STATUS_THR(yBlock);
+            const algorithmFPType * y = yBlock.get();
+
+            for (size_t i = 0; i < blockSize1; i++)
+            {
+                for (size_t j = 0; j < blockSize2; j++)
+                {
+                    algorithmFPType numerator = 0.0;
+                    algorithmFPType xNorm     = 0.0;
+                    algorithmFPType yNorm     = 0.0;
+
+                    for (size_t k = 0; k < p; k++)
+                    {
+                        algorithmFPType x_centered = x[i * p + k] - xMean[i];
+                        algorithmFPType y_centered = y[j * p + k] - yMean[shift2 + j];
+
+                        numerator += x_centered * y_centered;
+                        xNorm += x_centered * x_centered;
+                        yNorm += y_centered * y_centered;
+                    }
+
+                    algorithmFPType denominator = xNorm * yNorm;
+                    if (denominator > 0.0)
+                    {
+                        r[i * nVectors2 + shift2 + j] = 1.0
+                                                        - numerator
+                                                              / (daal::internal::MathInst<algorithmFPType, cpu>::sSqrt(xNorm)
+                                                                 * daal::internal::MathInst<algorithmFPType, cpu>::sSqrt(yNorm));
+                    }
+                    else
+                    {
+                        r[i * nVectors2 + shift2 + j] = 1.0; // Maximum distance when no variance
+                    }
+                }
+            }
+        }
+    });
+
+    return safeStat.detach();
+}
+
 } // namespace internal
 
 } // namespace correlation_distance

cpp/oneapi/dal/algo/BUILD

@@ -16,6 +16,7 @@ ALGOS = [
     "basic_statistics",
     "chebyshev_distance",
     "connected_components",
+    "correlation_distance",
     "cosine_distance",
     "dbscan",
     "decision_tree",

cpp/oneapi/dal/algo/correlation_distance.hpp

@@ -0,0 +1,19 @@
+/*******************************************************************************
+* Copyright contributors to the oneDAL project
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*******************************************************************************/
+
+#pragma once
+
+#include "oneapi/dal/algo/correlation_distance/compute.hpp"

cpp/oneapi/dal/algo/correlation_distance/BUILD

@@ -0,0 +1,38 @@
+package(default_visibility = ["//visibility:public"])
+load("@onedal//dev/bazel:dal.bzl",
+    "dal_module",
+    "dal_test_suite",
+)
+
+dal_module(
+    name = "correlation_distance",
+    auto = True,
+    dal_deps = [
+        "@onedal//cpp/oneapi/dal:core",
+	    "@onedal//cpp/oneapi/dal/backend/primitives:blas",
+	    "@onedal//cpp/oneapi/dal/backend/primitives:reduction",
+	    "@onedal//cpp/oneapi/dal/backend/primitives:distance",
+    ],
+    extra_deps = [
+	    "@onedal//cpp/daal/src/algorithms/cordistance:kernel",
+	    "@onedal//cpp/daal:data_management",
+    ]
+)
+
+dal_test_suite(
+    name = "interface_tests",
+    framework = "catch2",
+    srcs = glob([
+        "test/*.cpp",
+    ]),
+    dal_deps = [
+        ":correlation_distance",
+    ],
+)
+
+dal_test_suite(
+    name = "tests",
+    tests = [
+        ":interface_tests",
+    ],
+)

cpp/oneapi/dal/algo/correlation_distance/backend/cpu/compute_kernel.hpp

@@ -0,0 +1,40 @@
+/*******************************************************************************
+* Copyright contributors to the oneDAL project
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*******************************************************************************/
+
+#pragma once
+
+#include "oneapi/dal/algo/correlation_distance/compute_types.hpp"
+#include "oneapi/dal/backend/dispatcher.hpp"
+#include "oneapi/dal/table/homogen.hpp"
+
+namespace oneapi::dal::correlation_distance::backend {
+
+template <typename Float, typename Method, typename Task>
+struct compute_kernel_cpu {
+    compute_result<Task> operator()(const dal::backend::context_cpu& ctx,
+                                    const detail::descriptor_base<Task>& params,
+                                    const compute_input<Task>& input) const;
+
+#ifdef ONEDAL_DATA_PARALLEL
+    void operator()(const dal::backend::context_cpu& ctx,
+                    const detail::descriptor_base<Task>& params,
+                    const table& x,
+                    const table& y,
+                    homogen_table& res) const;
+#endif
+};
+
+} // namespace oneapi::dal::correlation_distance::backend

cpp/oneapi/dal/algo/correlation_distance/backend/cpu/compute_kernel_dense.cpp

@@ -0,0 +1,101 @@
+/*******************************************************************************
+* Copyright contributors to the oneDAL project
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*******************************************************************************/
+
+#include <daal/src/algorithms/cordistance/cordistance_kernel.h>
+
+#include "oneapi/dal/algo/correlation_distance/backend/cpu/compute_kernel.hpp"
+#include "oneapi/dal/backend/interop/common.hpp"
+#include "oneapi/dal/backend/interop/error_converter.hpp"
+#include "oneapi/dal/backend/interop/table_conversion.hpp"
+#include "oneapi/dal/exceptions.hpp"
+
+#include "oneapi/dal/table/row_accessor.hpp"
+
+namespace oneapi::dal::correlation_distance::backend {
+
+using dal::backend::context_cpu;
+using input_t = compute_input<task::compute>;
+using result_t = compute_result<task::compute>;
+using descriptor_t = detail::descriptor_base<task::compute>;
+
+namespace daal_correlation = daal::algorithms::correlation_distance;
+namespace interop = dal::backend::interop;
+
+template <typename Float, daal::CpuType Cpu>
+using daal_correlation_t =
+    daal_correlation::internal::DistanceKernel<Float, daal_correlation::defaultDense, Cpu>;
+
+template <typename Float>
+static result_t call_daal_kernel(const context_cpu& ctx,
+                                 const descriptor_t& desc,
+                                 const table& x,
+                                 const table& y) {
+    const std::int64_t row_count_x = x.get_row_count();
+    const std::int64_t row_count_y = y.get_row_count();
+
+    dal::detail::check_mul_overflow(row_count_x, row_count_y);
+    auto arr_values = array<Float>::empty(row_count_x * row_count_y);
+
+    const auto daal_x = interop::convert_to_daal_table<Float>(x);
+    const auto daal_y = interop::convert_to_daal_table<Float>(y);
+    const auto daal_values =
+        interop::convert_to_daal_homogen_table(arr_values, row_count_x, row_count_y);
+
+    daal::algorithms::Parameter param;
+    const daal::data_management::NumericTable* daal_input_tables[2] = { daal_x.get(),
+                                                                        daal_y.get() };
+    daal::data_management::NumericTable* daal_result_table[1] = { daal_values.get() };
+
+    interop::status_to_exception(
+        interop::call_daal_kernel<Float, daal_correlation_t>(ctx,
+                                                             2,
+                                                             daal_input_tables,
+                                                             1,
+                                                             daal_result_table,
+                                                             &param));
+
+    return result_t().set_values(
+        dal::detail::homogen_table_builder{}.reset(arr_values, row_count_x, row_count_y).build());
+}
+
+template <typename Float>
+static result_t compute(const context_cpu& ctx, const descriptor_t& desc, const input_t& input) {
+    return call_daal_kernel<Float>(ctx, desc, input.get_x(), input.get_y());
+}
+
+template <typename Float>
+struct compute_kernel_cpu<Float, method::dense, task::compute> {
+    result_t operator()(const context_cpu& ctx,
+                        const descriptor_t& desc,
+                        const input_t& input) const {
+        return compute<Float>(ctx, desc, input);
+    }
+
+#ifdef ONEDAL_DATA_PARALLEL
+    void operator()(const context_cpu& ctx,
+                    const descriptor_t& desc,
+                    const table& x,
+                    const table& y,
+                    homogen_table& res) const {
+        throw unimplemented(dal::detail::error_messages::method_not_implemented());
+    }
+#endif
+};
+
+template struct compute_kernel_cpu<float, method::dense, task::compute>;
+template struct compute_kernel_cpu<double, method::dense, task::compute>;
+
+} // namespace oneapi::dal::correlation_distance::backend