New OP: CDist (#1808)

Add a new op for scikit-learn converter. It's for scikit's cdist function:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html

Will add docs and shape-inference function later.
Will convert it to an ONNX function before pushing into ONNX.
This commit is contained in:
Changming Sun 2019-09-17 10:55:31 -07:00 committed by GitHub
parent d646dc73e4
commit dc03ce0278
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
9 changed files with 372 additions and 27 deletions

View file

@ -0,0 +1,16 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
#include "cdist.h"
namespace onnxruntime {
namespace contrib {
#define DEFINE_KERNEL(data_type) \
ONNX_OPERATOR_TYPED_KERNEL_EX(CDist, kMSDomain, 1, data_type, kCpuExecutionProvider, \
KernelDefBuilder().TypeConstraint("T", DataTypeImpl::GetTensorType<data_type>()), \
CDist<data_type>);
DEFINE_KERNEL(float);
DEFINE_KERNEL(double);
} // namespace contrib
} // namespace onnxruntime

View file

@ -0,0 +1,155 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
#pragma once
#include "core/common/common.h"
#include "core/util/distance.h"
#include "core/framework/op_kernel.h"
#include "core/framework/op_kernel_context_internal.h"
#include "core/util/math_cpuonly.h"
#include "assert.h"
#ifndef USE_OPENMP
#include "core/util/eigen_common_wrapper.h"
#endif
namespace onnxruntime {
namespace contrib {
// https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html
//\param a: matrix with shape of[ma,n]
//\param b: matrix with shape of[mb,n]
//\param dest: matrix with shape of [ma,mb]
template <typename T, typename ElemFunc>
void cdist_single_threaded(const T* a, const T* b, T* dest, size_t ma, size_t mb, size_t n) {
ElemFunc f;
for (size_t i = 0; i != ma; ++i) {
const T* a1 = a + n * i;
for (size_t j = 0; j != mb; ++j) {
const T* b1 = b + n * j;
*dest++ = f(a1, b1, n);
}
}
}
template <typename T, typename ElemFunc>
class CDistOneBlock {
public:
const T* a;
const T* b;
T* dest;
size_t mb;
size_t n;
CDistOneBlock(const T* a1, const T* b1, T* dest1, size_t mb1, size_t n1)
: a(a1), b(b1), dest(dest1), mb(mb1), n(n1) {}
void operator()(Eigen::Index start, Eigen::Index end) {
Eigen::Index mb_local = mb;
Eigen::Index i = start / mb_local;
Eigen::Index j = start - i * mb_local;
assert(i * mb_local + j == start);
Eigen::Index i_end = end / mb_local;
Eigen::Index j_end = end - i_end * mb_local;
assert(i_end * mb_local + j_end == end);
T* dest_local = dest + start;
ElemFunc f;
const T* a1 = a + n * i;
for (; i != i_end; ++i) {
a1 = a + n * i;
for (; j != mb_local; ++j) {
const T* b1 = b + n * j;
*dest_local++ = f(a1, b1, n);
}
j = 0;
}
a1 = a + n * i;
for (j = 0; j != j_end; ++j) {
const T* b1 = b + n * j;
*dest_local++ = f(a1, b1, n);
}
assert(dest_local == dest + end);
}
};
template <typename T, typename ElemFunc>
void cdist(const T* a, const T* b, T* dest, size_t ma, size_t mb, size_t n, concurrency::ThreadPool* tp) {
#ifndef USE_OPENMP
if (tp == nullptr) {
#else
(void)tp;
#endif
return cdist_single_threaded<T, ElemFunc>(a, b, dest, ma, mb, n);
#ifndef USE_OPENMP
}
Eigen::ThreadPoolDevice device(&tp->GetHandler(), tp->NumThreads());
device.parallelFor(ma * mb, Eigen::TensorOpCost(0, 0, static_cast<double>(3 * n)),
CDistOneBlock<T, ElemFunc>(a, b, dest, mb, n));
#endif
}
template <typename T>
class CDist final : public OpKernel {
private:
typedef void (*DistFunc)(const T* a, const T* b, T* dest, size_t ma, size_t mb, size_t n,
concurrency::ThreadPool* tp);
enum { EUCLIDEAN, SQEUCLIDEAN } mode_;
public:
CDist(const OpKernelInfo& info) : OpKernel(info) {
std::string metric;
ORT_ENFORCE(info.GetAttr<std::string>("metric", &metric).IsOK());
if (metric.compare("sqeuclidean") == 0)
mode_ = SQEUCLIDEAN;
else if (metric.compare("euclidean") == 0) {
mode_ = EUCLIDEAN;
} else
ORT_NOT_IMPLEMENTED();
}
common::Status Compute(OpKernelContext* context) const override {
auto ctx_internal = static_cast<OpKernelContextInternal*>(context);
concurrency::ThreadPool* tp = ctx_internal->GetOperatorThreadPool();
assert(context->InputCount() == 2);
const Tensor* A = context->Input<Tensor>(0);
const Tensor* B = context->Input<Tensor>(1);
const TensorShape& shape_a = A->Shape();
const TensorShape& shape_b = B->Shape();
if (shape_a.NumDimensions() != 2 || shape_a[1] <= 0) {
return ORT_MAKE_STATUS(ONNXRUNTIME, FAIL, "The first input of CDist kernel has wrong shape: ", shape_a);
}
if (shape_b.NumDimensions() != 2) {
return ORT_MAKE_STATUS(ONNXRUNTIME, FAIL, "The second input of CDist kernel has wrong shape: ", shape_b);
}
if (shape_a[1] != shape_b[1]) {
return ORT_MAKE_STATUS(ONNXRUNTIME, FAIL, "Input shape dimensions mismatch:", shape_a, " and ", shape_b);
}
TensorShape output_shape = {shape_a[0], shape_b[0]};
Tensor* C = context->Output(0, output_shape);
T* output = C->MutableData<T>();
switch (mode_) {
case EUCLIDEAN:
if (shape_a[1] >= 8)
cdist<T, EuclideanWithEigen<T> >(A->Data<T>(), B->Data<T>(), output, shape_a[0], shape_b[0], shape_a[1], tp);
else // for smaller vector size, a raw loop is better
cdist<T, Euclidean<T> >(A->Data<T>(), B->Data<T>(), output, shape_a[0], shape_b[0], shape_a[1], tp);
break;
case SQEUCLIDEAN:
if (shape_a[1] >= 8)
cdist<T, SqeuclideanWithEigen<T> >(A->Data<T>(), B->Data<T>(), output, shape_a[0], shape_b[0], shape_a[1],
tp);
else // for smaller vector size, a raw loop is better
cdist<T, Sqeuclidean<T> >(A->Data<T>(), B->Data<T>(), output, shape_a[0], shape_b[0], shape_a[1], tp);
break;
default:
return Status(ONNXRUNTIME, NOT_IMPLEMENTED);
}
return Status::OK();
}
};
} // namespace contrib
} // namespace onnxruntime

View file

@ -27,6 +27,8 @@ class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1,
class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, uint8_t, DequantizeLinear);
class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, int8_t, DequantizeLinear);
class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, uint8_t, QuantizeLinear);
class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, float, CDist);
class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, double, CDist);
// This section includes all opkernel declarations for former experimental ops which have now been removed from onnx.
// To maintain backward compatibility these are added as contrib ops.
@ -100,6 +102,8 @@ void RegisterCpuContribKernels(KernelRegistry& kernel_registry) {
BuildKernelCreateInfo<ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, uint8_t, DequantizeLinear)>,
BuildKernelCreateInfo<ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, int8_t, DequantizeLinear)>,
BuildKernelCreateInfo<ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, uint8_t, QuantizeLinear)>,
BuildKernelCreateInfo<ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, float, CDist)>,
BuildKernelCreateInfo<ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kMSDomain, 1, double, CDist)>,
// These ops were experimental ops in onnx domain which have been removed now. We add them here as
// contrib ops to main backward compatibility

View file

@ -1648,6 +1648,23 @@ Example 4:
output_counts = [1, 2, 2, 1]
)DOC");
//see:https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html
ONNX_CONTRIB_OPERATOR_SCHEMA(CDist)
.SetDomain(kMSDomain)
.SinceVersion(1)
.Attr("metric",
"The distance metric to use. If a string, the distance function can be \"braycurtis\", \"canberra\", "
"\"chebyshev\", \"cityblock\", \"correlation\", \"cosine\", \"dice\", \"euclidean\", \"hamming\", \"jaccard\", "
"\"jensenshannon\", \"kulsinski\", \"mahalanobis\", \"matching\", \"minkowski\", \"rogerstanimoto\", \"russellrao\", "
"\"seuclidean\", \"sokalmichener\", \"sokalsneath\", \"sqeuclidean\", \"wminkowski\", \"yule\".",
AttributeProto::STRING, std::string("sqeuclidean"))
.Input(0, "A", "2D matrix with shape (M,N)", "T")
.Input(1, "B", "2D matrix with shape (K,N)", "T")
.Output(0, "C",
"A 2D Matrix that represents the distance between each pair of the two collections of inputs.",
"T")
.TypeConstraint("T", {"tensor(float)", "tensor(double)"}, "Constrains input to only numeric types.");
ONNX_CONTRIB_OPERATOR_SCHEMA(CropAndResize)
.SetDomain(kMSDomain)
.SinceVersion(1)

View file

@ -0,0 +1,55 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
#pragma once
#include <cmath>
#include "math_cpuonly.h"
namespace onnxruntime {
// Computes the squared Euclidean distance between the vectors.
template <typename T>
class Sqeuclidean {
public:
T operator()(const T* a1, const T* b1, size_t n) const {
// if n is too small, Eigen is much slower than a plain loop
T sum = 0;
for (size_t k = 0; k != n; ++k) {
const T t = a1[k] - b1[k];
sum += t * t;
}
return sum;
}
};
// Computes the Euclidean distance between the vectors.
template <typename T>
class Euclidean {
public:
T operator()(const T* a1, const T* b1, size_t n) const {
// if n is too small, Eigen is much slower than a plain loop
T sum = 0;
for (size_t k = 0; k != n; ++k) {
const T t = a1[k] - b1[k];
sum += t * t;
}
return std::sqrt(sum);
}
};
template <typename T>
class SqeuclideanWithEigen {
public:
T operator()(const T* a1, const T* b1, size_t n) const {
return (ConstEigenVectorMap<T>(a1, n) - ConstEigenVectorMap<T>(b1, n)).array().square().sum();
}
};
template <typename T>
class EuclideanWithEigen {
public:
T operator()(const T* a1, const T* b1, size_t n) const {
return std::sqrt((ConstEigenVectorMap<T>(a1, n) - ConstEigenVectorMap<T>(b1, n)).array().square().sum());
}
};
} // namespace onnxruntime

View file

@ -1,18 +1,18 @@
/**
* Copyright (c) 2016-present, Facebook, Inc.
*
* 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.
*/
* Copyright (c) 2016-present, Facebook, Inc.
*
* 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.
*/
// Modifications Copyright (c) Microsoft.
#pragma once
@ -35,7 +35,7 @@
#endif
#include "Eigen/Dense"
#include "core/framework/tensor.h"
namespace onnxruntime {
// common Eigen types that we will often use
@ -56,16 +56,18 @@ using ConstEigenVectorMap = Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, 1>
template <typename T>
using ConstEigenVectorArrayMap = Eigen::Map<const Eigen::Array<T, Eigen::Dynamic, 1>>;
template <typename T>
using EigenMatrixMapRowMajor = Eigen::Map<
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>;
using EigenMatrixMapRowMajor = Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>;
template <typename T>
using ConstEigenMatrixMapRowMajor = Eigen::Map<
const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>;
using ConstEigenMatrixMapRowMajor = Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>;
template <typename T>
auto EigenMap(Tensor& t) { return EigenVectorMap<T>(t.template MutableData<T>(), t.Shape().Size()); }
auto EigenMap(Tensor& t) {
return EigenVectorMap<T>(t.template MutableData<T>(), t.Shape().Size());
}
template <typename T>
auto EigenMap(const Tensor& t) { return ConstEigenVectorMap<T>(t.template Data<T>(), t.Shape().Size()); }
auto EigenMap(const Tensor& t) {
return ConstEigenVectorMap<T>(t.template Data<T>(), t.Shape().Size());
}
class CPUMathUtil {
public:
@ -75,7 +77,7 @@ class CPUMathUtil {
static CPUMathUtil p;
return p;
}
//todo: the random generate interface.
// todo: the random generate interface.
private:
CPUMathUtil() = default;
};

View file

@ -0,0 +1,25 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
#include "core/util/distance.h"
#include <gtest/gtest.h>
using testing::Types;
namespace onnxruntime {
template <typename T>
class SqeuclideanTest : public ::testing::Test {
public:
T func;
};
using MyTypes = Types<Sqeuclidean<double>, SqeuclideanWithEigen<double> >;
TYPED_TEST_CASE(SqeuclideanTest, MyTypes);
TYPED_TEST(SqeuclideanTest, test1) {
double a = 10;
double b = 5;
ASSERT_DOUBLE_EQ(25, this->func(&a, &b, 1));
}
} // namespace onnxruntime

View file

@ -1,7 +1,6 @@
#!/usr/bin/env python3
# Copyright (c) Microsoft Corporation. All rights reserved.
# Copyright (c) Microsoft Corporation. All rights reserved.
# Licensed under the MIT License.
import onnx
import numpy as np
import os
@ -10,12 +9,27 @@ from onnx import numpy_helper
from onnx import helper
from onnx import utils
from onnx import AttributeProto, TensorProto, GraphProto
from scipy.spatial import distance
def parse_arguments():
parser = argparse.ArgumentParser()
parser.add_argument("--output_dir", required=True, help="Path to the build directory.")
return parser.parse_args()
def write_config(model_dir):
with open(os.path.join(model_dir,"config.txt"),"w") as f:
f.write("per_sample_tolerance:1e-6\n")
f.write("relative_per_sample_tolerance:1e-6\n")
def write_tensor(f, c,input_name=None):
tensor = numpy_helper.from_array(c)
if input_name:
tensor.name = input_name
body = tensor.SerializeToString()
f.write(body)
def generate_abs_op_test(type, X, top_test_folder):
for is_raw in [True, False]:
if is_raw:
@ -35,7 +49,7 @@ def generate_abs_op_test(type, X, top_test_folder):
node_def = helper.make_node('Abs', inputs=['X'], outputs=['Y'])
# Create the graph (GraphProto)
graph_def = helper.make_graph( [node_def], 'test-model', [X_INFO], [Y], [tensor_x])
graph_def = helper.make_graph([node_def], 'test-model', [X_INFO], [Y], [tensor_x])
# Create the model (ModelProto)
model_def = helper.make_model(graph_def, producer_name='onnx-example')
#final_model = onnx.utils.polish_model(model_def)
@ -59,7 +73,7 @@ def generate_size_op_test(type, X, test_folder):
node_def = helper.make_node('Size', inputs=['X'], outputs=['Y'])
# Create the graph (GraphProto)
graph_def = helper.make_graph( [node_def], 'test-model', [X_INFO], [Y], [tensor_x])
graph_def = helper.make_graph([node_def], 'test-model', [X_INFO], [Y], [tensor_x])
# Create the model (ModelProto)
model_def = helper.make_model(graph_def, producer_name='onnx-example')
final_model = onnx.utils.polish_model(model_def)
@ -69,6 +83,28 @@ def generate_size_op_test(type, X, test_folder):
with open(os.path.join(data_dir,"output_0.pb"),"wb") as f:
f.write(expected_output_tensor.SerializeToString())
def generate_reducesum_op_test(X, test_folder):
type = onnx.mapping.NP_TYPE_TO_TENSOR_TYPE[X.dtype]
data_dir = os.path.join(test_folder,"test_data_0")
os.makedirs(data_dir, exist_ok=True)
# Create one output (ValueInfoProto)
Y = helper.make_tensor_value_info('Y', type, [])
X_INFO = helper.make_tensor_value_info('X', type, X.shape)
tensor_x = onnx.helper.make_tensor(name='X', data_type=type, dims=X.shape, vals=X.ravel(),raw=False)
# Create a node (NodeProto)
node_def = helper.make_node('ReduceSum', inputs=['X'], outputs=['Y'], keepdims=0)
# Create the graph (GraphProto)
graph_def = helper.make_graph([node_def], 'test-model', [X_INFO], [Y], [tensor_x])
# Create the model (ModelProto)
model_def = helper.make_model(graph_def, producer_name='onnx-example')
final_model = onnx.utils.polish_model(model_def)
onnx.save(final_model, os.path.join(test_folder, 'model.onnx'))
expected_output_array = np.sum(X)
expected_output_tensor = numpy_helper.from_array(expected_output_array)
with open(os.path.join(data_dir,"output_0.pb"),"wb") as f:
f.write(expected_output_tensor.SerializeToString())
def test_abs(output_dir):
generate_abs_op_test(TensorProto.FLOAT, np.random.randn(3, 4, 5).astype(np.float32), os.path.join(output_dir,'test_abs_float'))
generate_abs_op_test(TensorProto.DOUBLE, np.random.randn(3, 4, 5).astype(np.float64), os.path.join(output_dir,'test_abs_double'))
@ -83,14 +119,48 @@ def test_abs(output_dir):
number_info = np.iinfo(np.uint64)
generate_abs_op_test(TensorProto.UINT64, np.uint64([0, 1, 20, number_info.max]), os.path.join(output_dir, 'test_abs_uint64'))
def test_reducesum(output_dir):
generate_reducesum_op_test(np.random.randn(3, 4, 5).astype(np.float32), os.path.join(output_dir, 'test_reducesum_random'))
def test_size(output_dir):
generate_size_op_test(TensorProto.FLOAT, np.random.randn(100, 3000, 10).astype(np.float32), os.path.join(output_dir,'test_size_float'))
generate_size_op_test(TensorProto.STRING, np.array(['abc', 'xy'], dtype=np.bytes_), os.path.join(output_dir,'test_size_string'))
def gen_cdist_test(output_dir, dtype, M, N, K):
for mode in ['euclidean', 'sqeuclidean']:
test_folder = os.path.join(output_dir,"test_cdist_%s_%s_%d_%d_%d" % (dtype.__name__, mode, M,N,K))
data_dir = os.path.join(test_folder, "test_data_0")
os.makedirs(data_dir, exist_ok=True)
a = np.random.randn(M, K).astype(dtype)
b = np.random.randn(N, K).astype(dtype)
type = onnx.mapping.NP_TYPE_TO_TENSOR_TYPE[a.dtype]
c = distance.cdist(a, b, mode).astype(dtype)
node_def = helper.make_node('CDist', inputs=['A','B'], outputs=['C'],domain="com.microsoft", metric=mode)
graph_def = helper.make_graph([node_def], 'test-model', [helper.make_tensor_value_info('A', type, a.shape),
helper.make_tensor_value_info('B', type, b.shape)],
[helper.make_tensor_value_info('C', type, c.shape)])
model_def = helper.make_model(graph_def, producer_name='onnx-example',opset_imports=[helper.make_opsetid("com.microsoft", 1)])
onnx.save(model_def, os.path.join(test_folder, 'model.onnx'))
with open(os.path.join(data_dir,"input_0.pb"),"wb") as f:
write_tensor(f, a, "A")
with open(os.path.join(data_dir,"input_1.pb"),"wb") as f:
write_tensor(f, b, "B")
with open(os.path.join(data_dir,"output_0.pb"),"wb") as f:
write_tensor(f, c, "C")
write_config(test_folder)
def test_cdist(output_dir):
for dtype in [np.float32, np.float64] :
gen_cdist_test(output_dir, dtype, 1000, 2000, 500)
gen_cdist_test(output_dir, dtype, 1000, 2000, 1)
gen_cdist_test(output_dir, dtype, 1, 1, 1)
args = parse_arguments()
os.makedirs(args.output_dir,exist_ok=True)
test_abs(args.output_dir)
test_size(args.output_dir)
test_reducesum(args.output_dir)
test_cdist(args.output_dir)

View file

@ -85,6 +85,7 @@ fi
if [ $DEVICE_TYPE = "Normal" ]; then
/usr/bin/python${PYTHON_VER} -m pip install --upgrade --force-reinstall sympy==1.1.1
fi
/usr/bin/python${PYTHON_VER} -m pip install --upgrade scipy
rm -rf /var/lib/apt/lists/*
if [ $DEVICE_TYPE = "Normal" ]; then