From 4ec06cd5cba3ee888a819bc2f30903955f325179 Mon Sep 17 00:00:00 2001 From: jim0421 Date: Sun, 24 May 2020 09:12:34 +0800 Subject: [PATCH] ENH scipy blas for svm kernel function (#16530) --- doc/whats_new/v0.24.rst | 8 + sklearn/ensemble/tests/test_bagging.py | 100 +++++++----- sklearn/svm/_libsvm.pxi | 15 +- sklearn/svm/_libsvm.pyx | 24 ++- sklearn/svm/_libsvm_sparse.pyx | 40 +++-- sklearn/svm/src/libsvm/LIBSVM_CHANGES | 2 +- .../svm/src/libsvm/_svm_cython_blas_helpers.h | 9 + sklearn/svm/src/libsvm/libsvm_helper.c | 16 +- sklearn/svm/src/libsvm/libsvm_sparse_helper.c | 15 +- sklearn/svm/src/libsvm/svm.cpp | 154 +++++++++--------- sklearn/svm/src/libsvm/svm.h | 21 +-- 11 files changed, 233 insertions(+), 171 deletions(-) create mode 100644 sklearn/svm/src/libsvm/_svm_cython_blas_helpers.h diff --git a/doc/whats_new/v0.24.rst b/doc/whats_new/v0.24.rst index d33be203a3087..8e87fc4a261a1 100644 --- a/doc/whats_new/v0.24.rst +++ b/doc/whats_new/v0.24.rst @@ -97,6 +97,14 @@ Changelog and by validating data out of loops. :pr:`17038` by :user:`Wenbo Zhao `. +:mod:`sklearn.svm` +.................... +- |Enhancement| invoke scipy blas api for svm kernel function in ``fit``, + ``predict`` and related methods of :class:`svm.SVC`, :class:`svm.NuSVC`, + :class:`svm.SVR`, :class:`svm.NuSVR`, :class:`OneClassSVM`. + :pr:`16530` by :user:`Shuhua Fan `. + + Code and Documentation Contributors ----------------------------------- diff --git a/sklearn/ensemble/tests/test_bagging.py b/sklearn/ensemble/tests/test_bagging.py index f87c8652a4368..8b557f50ef5e0 100644 --- a/sklearn/ensemble/tests/test_bagging.py +++ b/sklearn/ensemble/tests/test_bagging.py @@ -4,9 +4,11 @@ # Author: Gilles Louppe # License: BSD 3 clause +from itertools import product import numpy as np import joblib +import pytest from sklearn.base import BaseEstimator @@ -30,7 +32,7 @@ from sklearn.model_selection import train_test_split from sklearn.datasets import load_diabetes, load_iris, make_hastie_10_2 from sklearn.utils import check_random_state -from sklearn.preprocessing import FunctionTransformer +from sklearn.preprocessing import FunctionTransformer, scale from scipy.sparse import csc_matrix, csr_matrix @@ -74,7 +76,32 @@ def test_classification(): **params).fit(X_train, y_train).predict(X_test) -def test_sparse_classification(): +@pytest.mark.parametrize( + 'sparse_format, params, method', + product( + [csc_matrix, csr_matrix], + [{ + "max_samples": 0.5, + "max_features": 2, + "bootstrap": True, + "bootstrap_features": True + }, { + "max_samples": 1.0, + "max_features": 4, + "bootstrap": True, + "bootstrap_features": True + }, { + "max_features": 2, + "bootstrap": False, + "bootstrap_features": True + }, { + "max_samples": 0.5, + "bootstrap": True, + "bootstrap_features": False + }], + ['predict', 'predict_proba', + 'predict_log_proba', 'decision_function'])) +def test_sparse_classification(sparse_format, params, method): # Check classification for various parameter settings on sparse input. class CustomSVC(SVC): @@ -86,52 +113,35 @@ def fit(self, X, y): return self rng = check_random_state(0) - X_train, X_test, y_train, y_test = train_test_split(iris.data, + X_train, X_test, y_train, y_test = train_test_split(scale(iris.data), iris.target, random_state=rng) - parameter_sets = [ - {"max_samples": 0.5, - "max_features": 2, - "bootstrap": True, - "bootstrap_features": True}, - {"max_samples": 1.0, - "max_features": 4, - "bootstrap": True, - "bootstrap_features": True}, - {"max_features": 2, - "bootstrap": False, - "bootstrap_features": True}, - {"max_samples": 0.5, - "bootstrap": True, - "bootstrap_features": False}, - ] - for sparse_format in [csc_matrix, csr_matrix]: - X_train_sparse = sparse_format(X_train) - X_test_sparse = sparse_format(X_test) - for params in parameter_sets: - for f in ['predict', 'predict_proba', 'predict_log_proba', 'decision_function']: - # Trained on sparse format - sparse_classifier = BaggingClassifier( - base_estimator=CustomSVC(decision_function_shape='ovr'), - random_state=1, - **params - ).fit(X_train_sparse, y_train) - sparse_results = getattr(sparse_classifier, f)(X_test_sparse) - - # Trained on dense format - dense_classifier = BaggingClassifier( - base_estimator=CustomSVC(decision_function_shape='ovr'), - random_state=1, - **params - ).fit(X_train, y_train) - dense_results = getattr(dense_classifier, f)(X_test) - assert_array_almost_equal(sparse_results, dense_results) - - sparse_type = type(X_train_sparse) - types = [i.data_type_ for i in sparse_classifier.estimators_] - - assert all([t == sparse_type for t in types]) + X_train_sparse = sparse_format(X_train) + X_test_sparse = sparse_format(X_test) + # Trained on sparse format + sparse_classifier = BaggingClassifier( + base_estimator=CustomSVC(kernel="linear", + decision_function_shape='ovr'), + random_state=1, + **params + ).fit(X_train_sparse, y_train) + sparse_results = getattr(sparse_classifier, method)(X_test_sparse) + + # Trained on dense format + dense_classifier = BaggingClassifier( + base_estimator=CustomSVC(kernel="linear", + decision_function_shape='ovr'), + random_state=1, + **params + ).fit(X_train, y_train) + dense_results = getattr(dense_classifier, method)(X_test) + assert_array_almost_equal(sparse_results, dense_results) + + sparse_type = type(X_train_sparse) + types = [i.data_type_ for i in sparse_classifier.estimators_] + + assert all([t == sparse_type for t in types]) def test_regression(): diff --git a/sklearn/svm/_libsvm.pxi b/sklearn/svm/_libsvm.pxi index a3c8f1c33dd1e..ab7d212e3ba1a 100644 --- a/sklearn/svm/_libsvm.pxi +++ b/sklearn/svm/_libsvm.pxi @@ -1,5 +1,10 @@ ################################################################################ # Includes +cdef extern from "_svm_cython_blas_helpers.h": + ctypedef double (*dot_func)(int, double*, int, double*, int) + cdef struct BlasFunctions: + dot_func dot + cdef extern from "svm.h": cdef struct svm_node @@ -32,9 +37,9 @@ cdef extern from "svm.h": double *W # instance weights char *svm_check_parameter(svm_problem *, svm_parameter *) - svm_model *svm_train(svm_problem *, svm_parameter *, int *) nogil + svm_model *svm_train(svm_problem *, svm_parameter *, int *, BlasFunctions *) nogil void svm_free_and_destroy_model(svm_model** model_ptr_ptr) - void svm_cross_validation(svm_problem *, svm_parameter *, int nr_fold, double *target) nogil + void svm_cross_validation(svm_problem *, svm_parameter *, int nr_fold, double *target, BlasFunctions *) nogil cdef extern from "libsvm_helper.c": @@ -54,9 +59,9 @@ cdef extern from "libsvm_helper.c": void copy_intercept (char *, svm_model *, np.npy_intp *) void copy_SV (char *, svm_model *, np.npy_intp *) int copy_support (char *data, svm_model *model) - int copy_predict (char *, svm_model *, np.npy_intp *, char *) nogil - int copy_predict_proba (char *, svm_model *, np.npy_intp *, char *) nogil - int copy_predict_values(char *, svm_model *, np.npy_intp *, char *, int) nogil + int copy_predict (char *, svm_model *, np.npy_intp *, char *, BlasFunctions *) nogil + int copy_predict_proba (char *, svm_model *, np.npy_intp *, char *, BlasFunctions *) nogil + int copy_predict_values(char *, svm_model *, np.npy_intp *, char *, int, BlasFunctions *) nogil void copy_nSV (char *, svm_model *) void copy_probA (char *, svm_model *, np.npy_intp *) void copy_probB (char *, svm_model *, np.npy_intp *) diff --git a/sklearn/svm/_libsvm.pyx b/sklearn/svm/_libsvm.pyx index 079a791fef3b6..9488bda4ccf58 100644 --- a/sklearn/svm/_libsvm.pyx +++ b/sklearn/svm/_libsvm.pyx @@ -34,6 +34,7 @@ import warnings import numpy as np cimport numpy as np from libc.stdlib cimport free +from ..utils._cython_blas cimport _dot include "_libsvm.pxi" @@ -189,11 +190,12 @@ def fit( # for SVR: epsilon is called p in libsvm error_repl = error_msg.decode('utf-8').replace("p < 0", "epsilon < 0") raise ValueError(error_repl) - + cdef BlasFunctions blas_functions + blas_functions.dot = _dot[double] # this does the real work cdef int fit_status = 0 with nogil: - model = svm_train(&problem, ¶m, &fit_status) + model = svm_train(&problem, ¶m, &fit_status, &blas_functions) # from here until the end, we just copy the data returned by # svm_train @@ -352,12 +354,13 @@ def predict(np.ndarray[np.float64_t, ndim=2, mode='c'] X, model = set_model(¶m, nSV.shape[0], SV.data, SV.shape, support.data, support.shape, sv_coef.strides, sv_coef.data, intercept.data, nSV.data, probA.data, probB.data) - + cdef BlasFunctions blas_functions + blas_functions.dot = _dot[double] #TODO: use check_model try: dec_values = np.empty(X.shape[0]) with nogil: - rv = copy_predict(X.data, model, X.shape, dec_values.data) + rv = copy_predict(X.data, model, X.shape, dec_values.data, &blas_functions) if rv < 0: raise MemoryError("We've run out of memory") finally: @@ -457,10 +460,12 @@ def predict_proba( probA.data, probB.data) cdef np.npy_intp n_class = get_nr(model) + cdef BlasFunctions blas_functions + blas_functions.dot = _dot[double] try: dec_values = np.empty((X.shape[0], n_class), dtype=np.float64) with nogil: - rv = copy_predict_proba(X.data, model, X.shape, dec_values.data) + rv = copy_predict_proba(X.data, model, X.shape, dec_values.data, &blas_functions) if rv < 0: raise MemoryError("We've run out of memory") finally: @@ -561,11 +566,12 @@ def decision_function( else: n_class = get_nr(model) n_class = n_class * (n_class - 1) // 2 - + cdef BlasFunctions blas_functions + blas_functions.dot = _dot[double] try: dec_values = np.empty((X.shape[0], n_class), dtype=np.float64) with nogil: - rv = copy_predict_values(X.data, model, X.shape, dec_values.data, n_class) + rv = copy_predict_values(X.data, model, X.shape, dec_values.data, n_class, &blas_functions) if rv < 0: raise MemoryError("We've run out of memory") finally: @@ -704,10 +710,12 @@ def cross_validation( raise ValueError(error_msg) cdef np.ndarray[np.float64_t, ndim=1, mode='c'] target + cdef BlasFunctions blas_functions + blas_functions.dot = _dot[double] try: target = np.empty((X.shape[0]), dtype=np.float64) with nogil: - svm_cross_validation(&problem, ¶m, n_fold, target.data) + svm_cross_validation(&problem, ¶m, n_fold, target.data, &blas_functions) finally: free(problem.x) diff --git a/sklearn/svm/_libsvm_sparse.pyx b/sklearn/svm/_libsvm_sparse.pyx index 4b5070a64aad8..92d94b0c685a5 100644 --- a/sklearn/svm/_libsvm_sparse.pyx +++ b/sklearn/svm/_libsvm_sparse.pyx @@ -3,23 +3,27 @@ import numpy as np cimport numpy as np from scipy import sparse from ..exceptions import ConvergenceWarning - +from ..utils._cython_blas cimport _dot np.import_array() - cdef extern from *: ctypedef char* const_char_p "const char*" ################################################################################ # Includes +cdef extern from "_svm_cython_blas_helpers.h": + ctypedef double (*dot_func)(int, double*, int, double*, int) + cdef struct BlasFunctions: + dot_func dot + cdef extern from "svm.h": cdef struct svm_csr_node cdef struct svm_csr_model cdef struct svm_parameter cdef struct svm_csr_problem char *svm_csr_check_parameter(svm_csr_problem *, svm_parameter *) - svm_csr_model *svm_csr_train(svm_csr_problem *, svm_parameter *, int *) nogil + svm_csr_model *svm_csr_train(svm_csr_problem *, svm_parameter *, int *, BlasFunctions *) nogil void svm_csr_free_and_destroy_model(svm_csr_model** model_ptr_ptr) cdef extern from "libsvm_sparse_helper.c": @@ -39,18 +43,18 @@ cdef extern from "libsvm_sparse_helper.c": void copy_sv_coef (char *, svm_csr_model *) void copy_support (char *, svm_csr_model *) void copy_intercept (char *, svm_csr_model *, np.npy_intp *) - int copy_predict (char *, svm_csr_model *, np.npy_intp *, char *) + int copy_predict (char *, svm_csr_model *, np.npy_intp *, char *, BlasFunctions *) int csr_copy_predict_values (np.npy_intp *data_size, char *data, np.npy_intp *index_size, char *index, np.npy_intp *intptr_size, char *size, - svm_csr_model *model, char *dec_values, int nr_class) + svm_csr_model *model, char *dec_values, int nr_class, BlasFunctions *) int csr_copy_predict (np.npy_intp *data_size, char *data, np.npy_intp *index_size, char *index, np.npy_intp *intptr_size, char *size, - svm_csr_model *model, char *dec_values) nogil + svm_csr_model *model, char *dec_values, BlasFunctions *) nogil int csr_copy_predict_proba (np.npy_intp *data_size, char *data, np.npy_intp *index_size, char *index, np.npy_intp *intptr_size, char *size, - svm_csr_model *model, char *dec_values) nogil + svm_csr_model *model, char *dec_values, BlasFunctions *) nogil - int copy_predict_values(char *, svm_csr_model *, np.npy_intp *, char *, int) + int copy_predict_values(char *, svm_csr_model *, np.npy_intp *, char *, int, BlasFunctions *) int csr_copy_SV (char *values, np.npy_intp *n_indices, char *indices, np.npy_intp *n_indptr, char *indptr, svm_csr_model *model, int n_features) @@ -145,11 +149,12 @@ def libsvm_sparse_train ( int n_features, free_problem(problem) free_param(param) raise ValueError(error_msg) - + cdef BlasFunctions blas_functions + blas_functions.dot = _dot[double] # call svm_train, this does the real work cdef int fit_status = 0 with nogil: - model = svm_csr_train(problem, param, &fit_status) + model = svm_csr_train(problem, param, &fit_status, &blas_functions) cdef np.npy_intp SV_len = get_l(model) cdef np.npy_intp n_class = get_nr(model) @@ -275,11 +280,14 @@ def libsvm_sparse_predict (np.ndarray[np.float64_t, ndim=1, mode='c'] T_data, nSV.data, probA.data, probB.data) #TODO: use check_model dec_values = np.empty(T_indptr.shape[0]-1) + cdef BlasFunctions blas_functions + blas_functions.dot = _dot[double] with nogil: rv = csr_copy_predict(T_data.shape, T_data.data, T_indices.shape, T_indices.data, T_indptr.shape, T_indptr.data, - model, dec_values.data) + model, dec_values.data, + &blas_functions) if rv < 0: raise MemoryError("We've run out of memory") # free model and param @@ -331,11 +339,14 @@ def libsvm_sparse_predict_proba( cdef np.npy_intp n_class = get_nr(model) cdef int rv dec_values = np.empty((T_indptr.shape[0]-1, n_class), dtype=np.float64) + cdef BlasFunctions blas_functions + blas_functions.dot = _dot[double] with nogil: rv = csr_copy_predict_proba(T_data.shape, T_data.data, T_indices.shape, T_indices.data, T_indptr.shape, T_indptr.data, - model, dec_values.data) + model, dec_values.data, + &blas_functions) if rv < 0: raise MemoryError("We've run out of memory") # free model and param @@ -397,10 +408,13 @@ def libsvm_sparse_decision_function( n_class = n_class * (n_class - 1) // 2 dec_values = np.empty((T_indptr.shape[0] - 1, n_class), dtype=np.float64) + cdef BlasFunctions blas_functions + blas_functions.dot = _dot[double] if csr_copy_predict_values(T_data.shape, T_data.data, T_indices.shape, T_indices.data, T_indptr.shape, T_indptr.data, - model, dec_values.data, n_class) < 0: + model, dec_values.data, n_class, + &blas_functions) < 0: raise MemoryError("We've run out of memory") # free model and param free_model_SV(model) diff --git a/sklearn/svm/src/libsvm/LIBSVM_CHANGES b/sklearn/svm/src/libsvm/LIBSVM_CHANGES index c437720def7e1..bde6beaca2694 100644 --- a/sklearn/svm/src/libsvm/LIBSVM_CHANGES +++ b/sklearn/svm/src/libsvm/LIBSVM_CHANGES @@ -6,5 +6,5 @@ This is here mainly as checklist for incorporation of new versions of libsvm. * Add random_seed support and call to srand in fit function * Improved random number generator (fix on windows, enhancement on other platforms). See - + * invoke scipy blas api for svm kernel function to improve performance with speedup rate of 1.5X to 2X for dense data only. See The changes made with respect to upstream are detailed in the heading of svm.cpp diff --git a/sklearn/svm/src/libsvm/_svm_cython_blas_helpers.h b/sklearn/svm/src/libsvm/_svm_cython_blas_helpers.h new file mode 100644 index 0000000000000..057e08195e9a5 --- /dev/null +++ b/sklearn/svm/src/libsvm/_svm_cython_blas_helpers.h @@ -0,0 +1,9 @@ +#ifndef _SVM_CYTHON_BLAS_HELPERS_H +#define _SVM_CYTHON_BLAS_HELPERS_H + +typedef double (*dot_func)(int, double*, int, double*, int); +typedef struct BlasFunctions{ + dot_func dot; +} BlasFunctions; + +#endif diff --git a/sklearn/svm/src/libsvm/libsvm_helper.c b/sklearn/svm/src/libsvm/libsvm_helper.c index b61cfd2c51b58..17f328f9e7c4c 100644 --- a/sklearn/svm/src/libsvm/libsvm_helper.c +++ b/sklearn/svm/src/libsvm/libsvm_helper.c @@ -1,7 +1,7 @@ #include #include #include "svm.h" - +#include "_svm_cython_blas_helpers.h" /* * Some helper methods for libsvm bindings. * @@ -293,7 +293,7 @@ void copy_probB(char *data, struct svm_model *model, npy_intp * dims) * It will return -1 if we run out of memory. */ int copy_predict(char *predict, struct svm_model *model, npy_intp *predict_dims, - char *dec_values) + char *dec_values, BlasFunctions *blas_functions) { double *t = (double *) dec_values; struct svm_node *predict_nodes; @@ -304,7 +304,7 @@ int copy_predict(char *predict, struct svm_model *model, npy_intp *predict_dims, if (predict_nodes == NULL) return -1; for(i=0; i #include #include "svm.h" - +#include "_svm_cython_blas_helpers.h" /* * Convert scipy.sparse.csr to libsvm's sparse data structure @@ -246,7 +246,7 @@ npy_intp get_nonzero_SV (struct svm_csr_model *model) { */ int csr_copy_predict (npy_intp *data_size, char *data, npy_intp *index_size, char *index, npy_intp *intptr_size, char *intptr, struct svm_csr_model *model, - char *dec_values) { + char *dec_values, BlasFunctions *blas_functions) { double *t = (double *) dec_values; struct svm_csr_node **predict_nodes; npy_intp i; @@ -257,7 +257,7 @@ int csr_copy_predict (npy_intp *data_size, char *data, npy_intp *index_size, if (predict_nodes == NULL) return -1; for(i=0; i < intptr_size[0] - 1; ++i) { - *t = svm_csr_predict(model, predict_nodes[i]); + *t = svm_csr_predict(model, predict_nodes[i], blas_functions); free(predict_nodes[i]); ++t; } @@ -267,7 +267,7 @@ int csr_copy_predict (npy_intp *data_size, char *data, npy_intp *index_size, int csr_copy_predict_values (npy_intp *data_size, char *data, npy_intp *index_size, char *index, npy_intp *intptr_size, char *intptr, struct svm_csr_model *model, - char *dec_values, int nr_class) { + char *dec_values, int nr_class, BlasFunctions *blas_functions) { struct svm_csr_node **predict_nodes; npy_intp i; @@ -278,7 +278,8 @@ int csr_copy_predict_values (npy_intp *data_size, char *data, npy_intp *index_si return -1; for(i=0; i < intptr_size[0] - 1; ++i) { svm_csr_predict_values(model, predict_nodes[i], - ((double *) dec_values) + i*nr_class); + ((double *) dec_values) + i*nr_class, + blas_functions); free(predict_nodes[i]); } free(predict_nodes); @@ -288,7 +289,7 @@ int csr_copy_predict_values (npy_intp *data_size, char *data, npy_intp *index_si int csr_copy_predict_proba (npy_intp *data_size, char *data, npy_intp *index_size, char *index, npy_intp *intptr_size, char *intptr, struct svm_csr_model *model, - char *dec_values) { + char *dec_values, BlasFunctions *blas_functions) { struct svm_csr_node **predict_nodes; npy_intp i; @@ -301,7 +302,7 @@ int csr_copy_predict_proba (npy_intp *data_size, char *data, npy_intp *index_siz return -1; for(i=0; i < intptr_size[0] - 1; ++i) { svm_csr_predict_probability( - model, predict_nodes[i], ((double *) dec_values) + i*m); + model, predict_nodes[i], ((double *) dec_values) + i*m, blas_functions); free(predict_nodes[i]); } free(predict_nodes); diff --git a/sklearn/svm/src/libsvm/svm.cpp b/sklearn/svm/src/libsvm/svm.cpp index a5f735d8111cf..d209e35fc0a35 100644 --- a/sklearn/svm/src/libsvm/svm.cpp +++ b/sklearn/svm/src/libsvm/svm.cpp @@ -67,8 +67,10 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include #include #include "svm.h" +#include "_svm_cython_blas_helpers.h" #include "../newrand/newrand.h" + #ifndef _LIBSVM_CPP typedef float Qfloat; typedef signed char schar; @@ -289,14 +291,14 @@ class QMatrix { class Kernel: public QMatrix { public: #ifdef _DENSE_REP - Kernel(int l, PREFIX(node) * x, const svm_parameter& param); + Kernel(int l, PREFIX(node) * x, const svm_parameter& param, BlasFunctions *blas_functions); #else - Kernel(int l, PREFIX(node) * const * x, const svm_parameter& param); + Kernel(int l, PREFIX(node) * const * x, const svm_parameter& param, BlasFunctions *blas_functions); #endif virtual ~Kernel(); static double k_function(const PREFIX(node) *x, const PREFIX(node) *y, - const svm_parameter& param); + const svm_parameter& param, BlasFunctions *blas_functions); virtual Qfloat *get_Q(int column, int len) const = 0; virtual double *get_QD() const = 0; virtual void swap_index(int i, int j) const // no so const... @@ -315,6 +317,8 @@ class Kernel: public QMatrix { const PREFIX(node) **x; #endif double *x_square; + // scipy blas pointer + BlasFunctions *m_blas; // svm_parameter const int kernel_type; @@ -322,26 +326,26 @@ class Kernel: public QMatrix { const double gamma; const double coef0; - static double dot(const PREFIX(node) *px, const PREFIX(node) *py); + static double dot(const PREFIX(node) *px, const PREFIX(node) *py, BlasFunctions *blas_functions); #ifdef _DENSE_REP - static double dot(const PREFIX(node) &px, const PREFIX(node) &py); + static double dot(const PREFIX(node) &px, const PREFIX(node) &py, BlasFunctions *blas_functions); #endif double kernel_linear(int i, int j) const { - return dot(x[i],x[j]); + return dot(x[i],x[j],m_blas); } double kernel_poly(int i, int j) const { - return powi(gamma*dot(x[i],x[j])+coef0,degree); + return powi(gamma*dot(x[i],x[j],m_blas)+coef0,degree); } double kernel_rbf(int i, int j) const { - return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j]))); + return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j],m_blas))); } double kernel_sigmoid(int i, int j) const { - return tanh(gamma*dot(x[i],x[j])+coef0); + return tanh(gamma*dot(x[i],x[j],m_blas)+coef0); } double kernel_precomputed(int i, int j) const { @@ -354,13 +358,14 @@ class Kernel: public QMatrix { }; #ifdef _DENSE_REP -Kernel::Kernel(int l, PREFIX(node) * x_, const svm_parameter& param) +Kernel::Kernel(int l, PREFIX(node) * x_, const svm_parameter& param, BlasFunctions *blas_functions) #else -Kernel::Kernel(int l, PREFIX(node) * const * x_, const svm_parameter& param) +Kernel::Kernel(int l, PREFIX(node) * const * x_, const svm_parameter& param, BlasFunctions *blas_functions) #endif :kernel_type(param.kernel_type), degree(param.degree), gamma(param.gamma), coef0(param.coef0) { + m_blas = blas_functions; switch(kernel_type) { case LINEAR: @@ -386,7 +391,7 @@ Kernel::Kernel(int l, PREFIX(node) * const * x_, const svm_parameter& param) { x_square = new double[l]; for(int i=0;idim, py->dim); - for (int i = 0; i < dim; i++) - sum += (px->values)[i] * (py->values)[i]; + sum = blas_functions->dot(dim, px->values, 1, py->values, 1); return sum; } -double Kernel::dot(const PREFIX(node) &px, const PREFIX(node) &py) +double Kernel::dot(const PREFIX(node) &px, const PREFIX(node) &py, BlasFunctions *blas_functions) { double sum = 0; int dim = min(px.dim, py.dim); - for (int i = 0; i < dim; i++) - sum += px.values[i] * py.values[i]; + sum = blas_functions->dot(dim, px.values, 1, py.values, 1); return sum; } #else -double Kernel::dot(const PREFIX(node) *px, const PREFIX(node) *py) +double Kernel::dot(const PREFIX(node) *px, const PREFIX(node) *py, BlasFunctions *blas_functions) { double sum = 0; while(px->index != -1 && py->index != -1) @@ -443,24 +446,26 @@ double Kernel::dot(const PREFIX(node) *px, const PREFIX(node) *py) #endif double Kernel::k_function(const PREFIX(node) *x, const PREFIX(node) *y, - const svm_parameter& param) + const svm_parameter& param, BlasFunctions *blas_functions) { switch(param.kernel_type) { case LINEAR: - return dot(x,y); + return dot(x,y,blas_functions); case POLY: - return powi(param.gamma*dot(x,y)+param.coef0,param.degree); + return powi(param.gamma*dot(x,y,blas_functions)+param.coef0,param.degree); case RBF: { double sum = 0; #ifdef _DENSE_REP int dim = min(x->dim, y->dim), i; + double* m_array = (double*)malloc(sizeof(double)*dim); for (i = 0; i < dim; i++) { - double d = x->values[i] - y->values[i]; - sum += d*d; + m_array[i] = x->values[i] - y->values[i]; } + sum = blas_functions->dot(dim, m_array, 1, m_array, 1); + free(m_array); for (; i < x->dim; i++) sum += x->values[i] * x->values[i]; for (; i < y->dim; i++) @@ -505,7 +510,7 @@ double Kernel::k_function(const PREFIX(node) *x, const PREFIX(node) *y, return exp(-param.gamma*sum); } case SIGMOID: - return tanh(param.gamma*dot(x,y)+param.coef0); + return tanh(param.gamma*dot(x,y,blas_functions)+param.coef0); case PRECOMPUTED: //x: test (validation), y: SV { #ifdef _DENSE_REP @@ -518,7 +523,6 @@ double Kernel::k_function(const PREFIX(node) *x, const PREFIX(node) *y, return 0; // Unreachable } } - // An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918 // Solves: // @@ -1411,8 +1415,8 @@ double Solver_NU::calculate_rho() class SVC_Q: public Kernel { public: - SVC_Q(const PREFIX(problem)& prob, const svm_parameter& param, const schar *y_) - :Kernel(prob.l, prob.x, param) + SVC_Q(const PREFIX(problem)& prob, const svm_parameter& param, const schar *y_, BlasFunctions *blas_functions) + :Kernel(prob.l, prob.x, param, blas_functions) { clone(y,y_,prob.l); cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20))); @@ -1461,8 +1465,8 @@ class SVC_Q: public Kernel class ONE_CLASS_Q: public Kernel { public: - ONE_CLASS_Q(const PREFIX(problem)& prob, const svm_parameter& param) - :Kernel(prob.l, prob.x, param) + ONE_CLASS_Q(const PREFIX(problem)& prob, const svm_parameter& param, BlasFunctions *blas_functions) + :Kernel(prob.l, prob.x, param, blas_functions) { cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20))); QD = new double[prob.l]; @@ -1507,8 +1511,8 @@ class ONE_CLASS_Q: public Kernel class SVR_Q: public Kernel { public: - SVR_Q(const PREFIX(problem)& prob, const svm_parameter& param) - :Kernel(prob.l, prob.x, param) + SVR_Q(const PREFIX(problem)& prob, const svm_parameter& param, BlasFunctions *blas_functions) + :Kernel(prob.l, prob.x, param, blas_functions) { l = prob.l; cache = new Cache(l,(long int)(param.cache_size*(1<<20))); @@ -1584,7 +1588,7 @@ class SVR_Q: public Kernel // static void solve_c_svc( const PREFIX(problem) *prob, const svm_parameter* param, - double *alpha, Solver::SolutionInfo* si, double Cp, double Cn) + double *alpha, Solver::SolutionInfo* si, double Cp, double Cn, BlasFunctions *blas_functions) { int l = prob->l; double *minus_ones = new double[l]; @@ -1610,7 +1614,7 @@ static void solve_c_svc( } Solver s; - s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y, + s.Solve(l, SVC_Q(*prob,*param,y, blas_functions), minus_ones, y, alpha, C, param->eps, si, param->shrinking, param->max_iter); @@ -1633,7 +1637,7 @@ static void solve_c_svc( static void solve_nu_svc( const PREFIX(problem) *prob, const svm_parameter *param, - double *alpha, Solver::SolutionInfo* si) + double *alpha, Solver::SolutionInfo* si, BlasFunctions *blas_functions) { int i; int l = prob->l; @@ -1675,7 +1679,7 @@ static void solve_nu_svc( zeros[i] = 0; Solver_NU s; - s.Solve(l, SVC_Q(*prob,*param,y), zeros, y, + s.Solve(l, SVC_Q(*prob,*param,y,blas_functions), zeros, y, alpha, C, param->eps, si, param->shrinking, param->max_iter); double r = si->r; @@ -1697,7 +1701,7 @@ static void solve_nu_svc( static void solve_one_class( const PREFIX(problem) *prob, const svm_parameter *param, - double *alpha, Solver::SolutionInfo* si) + double *alpha, Solver::SolutionInfo* si, BlasFunctions *blas_functions) { int l = prob->l; double *zeros = new double[l]; @@ -1730,7 +1734,7 @@ static void solve_one_class( } Solver s; - s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones, + s.Solve(l, ONE_CLASS_Q(*prob,*param,blas_functions), zeros, ones, alpha, C, param->eps, si, param->shrinking, param->max_iter); delete[] C; @@ -1740,7 +1744,7 @@ static void solve_one_class( static void solve_epsilon_svr( const PREFIX(problem) *prob, const svm_parameter *param, - double *alpha, Solver::SolutionInfo* si) + double *alpha, Solver::SolutionInfo* si, BlasFunctions *blas_functions) { int l = prob->l; double *alpha2 = new double[2*l]; @@ -1763,7 +1767,7 @@ static void solve_epsilon_svr( } Solver s; - s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y, + s.Solve(2*l, SVR_Q(*prob,*param,blas_functions), linear_term, y, alpha2, C, param->eps, si, param->shrinking, param->max_iter); double sum_alpha = 0; @@ -1782,7 +1786,7 @@ static void solve_epsilon_svr( static void solve_nu_svr( const PREFIX(problem) *prob, const svm_parameter *param, - double *alpha, Solver::SolutionInfo* si) + double *alpha, Solver::SolutionInfo* si, BlasFunctions *blas_functions) { int l = prob->l; double *C = new double[2*l]; @@ -1812,7 +1816,7 @@ static void solve_nu_svr( } Solver_NU s; - s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y, + s.Solve(2*l, SVR_Q(*prob,*param,blas_functions), linear_term, y, alpha2, C, param->eps, si, param->shrinking, param->max_iter); info("epsilon = %f\n",-si->r); @@ -1837,7 +1841,7 @@ struct decision_function static decision_function svm_train_one( const PREFIX(problem) *prob, const svm_parameter *param, - double Cp, double Cn, int *status) + double Cp, double Cn, int *status, BlasFunctions *blas_functions) { double *alpha = Malloc(double,prob->l); Solver::SolutionInfo si; @@ -1845,23 +1849,23 @@ static decision_function svm_train_one( { case C_SVC: si.upper_bound = Malloc(double,prob->l); - solve_c_svc(prob,param,alpha,&si,Cp,Cn); + solve_c_svc(prob,param,alpha,&si,Cp,Cn,blas_functions); break; case NU_SVC: si.upper_bound = Malloc(double,prob->l); - solve_nu_svc(prob,param,alpha,&si); + solve_nu_svc(prob,param,alpha,&si,blas_functions); break; case ONE_CLASS: si.upper_bound = Malloc(double,prob->l); - solve_one_class(prob,param,alpha,&si); + solve_one_class(prob,param,alpha,&si,blas_functions); break; case EPSILON_SVR: si.upper_bound = Malloc(double,2*prob->l); - solve_epsilon_svr(prob,param,alpha,&si); + solve_epsilon_svr(prob,param,alpha,&si,blas_functions); break; case NU_SVR: si.upper_bound = Malloc(double,2*prob->l); - solve_nu_svr(prob,param,alpha,&si); + solve_nu_svr(prob,param,alpha,&si,blas_functions); break; } @@ -2092,7 +2096,7 @@ static void multiclass_probability(int k, double **r, double *p) // Cross-validation decision values for probability estimates static void svm_binary_svc_probability( const PREFIX(problem) *prob, const svm_parameter *param, - double Cp, double Cn, double& probA, double& probB, int * status) + double Cp, double Cn, double& probA, double& probB, int * status, BlasFunctions *blas_functions) { int i; int nr_fold = 5; @@ -2165,13 +2169,13 @@ static void svm_binary_svc_probability( subparam.weight_label[1]=-1; subparam.weight[0]=Cp; subparam.weight[1]=Cn; - struct PREFIX(model) *submodel = PREFIX(train)(&subprob,&subparam, status); + struct PREFIX(model) *submodel = PREFIX(train)(&subprob,&subparam, status, blas_functions); for(j=begin;jx+perm[j]),&(dec_values[perm[j]])); + PREFIX(predict_values)(submodel,(prob->x+perm[j]),&(dec_values[perm[j]]), blas_functions); #else - PREFIX(predict_values)(submodel,prob->x[perm[j]],&(dec_values[perm[j]])); + PREFIX(predict_values)(submodel,prob->x[perm[j]],&(dec_values[perm[j]]), blas_functions); #endif // ensure +1 -1 order; reason not using CV subroutine dec_values[perm[j]] *= submodel->label[0]; @@ -2190,7 +2194,7 @@ static void svm_binary_svc_probability( // Return parameter of a Laplace distribution static double svm_svr_probability( - const PREFIX(problem) *prob, const svm_parameter *param) + const PREFIX(problem) *prob, const svm_parameter *param, BlasFunctions *blas_functions) { int i; int nr_fold = 5; @@ -2201,7 +2205,7 @@ static double svm_svr_probability( newparam.probability = 0; newparam.random_seed = -1; // This is called from train, which already sets // the seed. - PREFIX(cross_validation)(prob,&newparam,nr_fold,ymv); + PREFIX(cross_validation)(prob,&newparam,nr_fold,ymv, blas_functions); for(i=0;il;i++) { ymv[i]=prob->y[i]-ymv[i]; @@ -2346,7 +2350,7 @@ static void remove_zero_weight(PREFIX(problem) *newprob, const PREFIX(problem) * // Interface functions // PREFIX(model) *PREFIX(train)(const PREFIX(problem) *prob, const svm_parameter *param, - int *status) + int *status, BlasFunctions *blas_functions) { PREFIX(problem) newprob; remove_zero_weight(&newprob, prob); @@ -2377,10 +2381,10 @@ PREFIX(model) *PREFIX(train)(const PREFIX(problem) *prob, const svm_parameter *p param->svm_type == NU_SVR)) { model->probA = Malloc(double,1); - model->probA[0] = NAMESPACE::svm_svr_probability(prob,param); + model->probA[0] = NAMESPACE::svm_svr_probability(prob,param,blas_functions); } - NAMESPACE::decision_function f = NAMESPACE::svm_train_one(prob,param,0,0, status); + NAMESPACE::decision_function f = NAMESPACE::svm_train_one(prob,param,0,0, status,blas_functions); model->rho = Malloc(double,1); model->rho[0] = f.rho; @@ -2495,9 +2499,9 @@ PREFIX(model) *PREFIX(train)(const PREFIX(problem) *prob, const svm_parameter *p } if(param->probability) - NAMESPACE::svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p], status); + NAMESPACE::svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p], status, blas_functions); - f[p] = NAMESPACE::svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j], status); + f[p] = NAMESPACE::svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j], status, blas_functions); for(k=0;k 0) nonzero[si+k] = true; @@ -2629,7 +2633,7 @@ PREFIX(model) *PREFIX(train)(const PREFIX(problem) *prob, const svm_parameter *p } // Stratified cross validation -void PREFIX(cross_validation)(const PREFIX(problem) *prob, const svm_parameter *param, int nr_fold, double *target) +void PREFIX(cross_validation)(const PREFIX(problem) *prob, const svm_parameter *param, int nr_fold, double *target, BlasFunctions *blas_functions) { int i; int *fold_start = Malloc(int,nr_fold+1); @@ -2736,25 +2740,25 @@ void PREFIX(cross_validation)(const PREFIX(problem) *prob, const svm_parameter * ++k; } int dummy_status = 0; // IGNORES TIMEOUT ERRORS - struct PREFIX(model) *submodel = PREFIX(train)(&subprob,param, &dummy_status); + struct PREFIX(model) *submodel = PREFIX(train)(&subprob,param, &dummy_status, blas_functions); if(param->probability && (param->svm_type == C_SVC || param->svm_type == NU_SVC)) { double *prob_estimates=Malloc(double, PREFIX(get_nr_class)(submodel)); for(j=begin;jx + perm[j]),prob_estimates); + target[perm[j]] = PREFIX(predict_probability)(submodel,(prob->x + perm[j]),prob_estimates, blas_functions); #else - target[perm[j]] = PREFIX(predict_probability)(submodel,prob->x[perm[j]],prob_estimates); + target[perm[j]] = PREFIX(predict_probability)(submodel,prob->x[perm[j]],prob_estimates, blas_functions); #endif free(prob_estimates); } else for(j=begin;jx+perm[j]); + target[perm[j]] = PREFIX(predict)(submodel,prob->x+perm[j],blas_functions); #else - target[perm[j]] = PREFIX(predict)(submodel,prob->x[perm[j]]); + target[perm[j]] = PREFIX(predict)(submodel,prob->x[perm[j]],blas_functions); #endif PREFIX(free_and_destroy_model)(&submodel); free(subprob.x); @@ -2795,7 +2799,7 @@ double PREFIX(get_svr_probability)(const PREFIX(model) *model) } } -double PREFIX(predict_values)(const PREFIX(model) *model, const PREFIX(node) *x, double* dec_values) +double PREFIX(predict_values)(const PREFIX(model) *model, const PREFIX(node) *x, double* dec_values, BlasFunctions *blas_functions) { int i; if(model->param.svm_type == ONE_CLASS || @@ -2807,9 +2811,9 @@ double PREFIX(predict_values)(const PREFIX(model) *model, const PREFIX(node) *x, for(i=0;il;i++) #ifdef _DENSE_REP - sum += sv_coef[i] * NAMESPACE::Kernel::k_function(x,model->SV+i,model->param); + sum += sv_coef[i] * NAMESPACE::Kernel::k_function(x,model->SV+i,model->param,blas_functions); #else - sum += sv_coef[i] * NAMESPACE::Kernel::k_function(x,model->SV[i],model->param); + sum += sv_coef[i] * NAMESPACE::Kernel::k_function(x,model->SV[i],model->param,blas_functions); #endif sum -= model->rho[0]; *dec_values = sum; @@ -2827,9 +2831,9 @@ double PREFIX(predict_values)(const PREFIX(model) *model, const PREFIX(node) *x, double *kvalue = Malloc(double,l); for(i=0;iSV+i,model->param); + kvalue[i] = NAMESPACE::Kernel::k_function(x,model->SV+i,model->param,blas_functions); #else - kvalue[i] = NAMESPACE::Kernel::k_function(x,model->SV[i],model->param); + kvalue[i] = NAMESPACE::Kernel::k_function(x,model->SV[i],model->param,blas_functions); #endif int *start = Malloc(int,nr_class); @@ -2880,7 +2884,7 @@ double PREFIX(predict_values)(const PREFIX(model) *model, const PREFIX(node) *x, } } -double PREFIX(predict)(const PREFIX(model) *model, const PREFIX(node) *x) +double PREFIX(predict)(const PREFIX(model) *model, const PREFIX(node) *x, BlasFunctions *blas_functions) { int nr_class = model->nr_class; double *dec_values; @@ -2890,13 +2894,13 @@ double PREFIX(predict)(const PREFIX(model) *model, const PREFIX(node) *x) dec_values = Malloc(double, 1); else dec_values = Malloc(double, nr_class*(nr_class-1)/2); - double pred_result = PREFIX(predict_values)(model, x, dec_values); + double pred_result = PREFIX(predict_values)(model, x, dec_values, blas_functions); free(dec_values); return pred_result; } double PREFIX(predict_probability)( - const PREFIX(model) *model, const PREFIX(node) *x, double *prob_estimates) + const PREFIX(model) *model, const PREFIX(node) *x, double *prob_estimates, BlasFunctions *blas_functions) { if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) && model->probA!=NULL && model->probB!=NULL) @@ -2904,7 +2908,7 @@ double PREFIX(predict_probability)( int i; int nr_class = model->nr_class; double *dec_values = Malloc(double, nr_class*(nr_class-1)/2); - PREFIX(predict_values)(model, x, dec_values); + PREFIX(predict_values)(model, x, dec_values, blas_functions); double min_prob=1e-7; double **pairwise_prob=Malloc(double *,nr_class); @@ -2931,7 +2935,7 @@ double PREFIX(predict_probability)( return model->label[prob_max_idx]; } else - return PREFIX(predict)(model, x); + return PREFIX(predict)(model, x, blas_functions); } diff --git a/sklearn/svm/src/libsvm/svm.h b/sklearn/svm/src/libsvm/svm.h index 4002a77c93ac4..0e509c61c37ed 100644 --- a/sklearn/svm/src/libsvm/svm.h +++ b/sklearn/svm/src/libsvm/svm.h @@ -6,6 +6,7 @@ #ifdef __cplusplus extern "C" { #endif +#include "_svm_cython_blas_helpers.h" struct svm_node { @@ -118,8 +119,8 @@ struct svm_csr_model }; -struct svm_model *svm_train(const struct svm_problem *prob, const struct svm_parameter *param, int *status); -void svm_cross_validation(const struct svm_problem *prob, const struct svm_parameter *param, int nr_fold, double *target); +struct svm_model *svm_train(const struct svm_problem *prob, const struct svm_parameter *param, int *status, BlasFunctions *blas_functions); +void svm_cross_validation(const struct svm_problem *prob, const struct svm_parameter *param, int nr_fold, double *target, BlasFunctions *blas_functions); int svm_save_model(const char *model_file_name, const struct svm_model *model); struct svm_model *svm_load_model(const char *model_file_name); @@ -129,9 +130,9 @@ int svm_get_nr_class(const struct svm_model *model); void svm_get_labels(const struct svm_model *model, int *label); double svm_get_svr_probability(const struct svm_model *model); -double svm_predict_values(const struct svm_model *model, const struct svm_node *x, double* dec_values); -double svm_predict(const struct svm_model *model, const struct svm_node *x); -double svm_predict_probability(const struct svm_model *model, const struct svm_node *x, double* prob_estimates); +double svm_predict_values(const struct svm_model *model, const struct svm_node *x, double* dec_values, BlasFunctions *blas_functions); +double svm_predict(const struct svm_model *model, const struct svm_node *x, BlasFunctions *blas_functions); +double svm_predict_probability(const struct svm_model *model, const struct svm_node *x, double* prob_estimates, BlasFunctions *blas_functions); void svm_free_model_content(struct svm_model *model_ptr); void svm_free_and_destroy_model(struct svm_model **model_ptr_ptr); @@ -144,17 +145,17 @@ void svm_set_print_string_function(void (*print_func)(const char *)); /* sparse version */ -struct svm_csr_model *svm_csr_train(const struct svm_csr_problem *prob, const struct svm_parameter *param, int *status); -void svm_csr_cross_validation(const struct svm_csr_problem *prob, const struct svm_parameter *param, int nr_fold, double *target); +struct svm_csr_model *svm_csr_train(const struct svm_csr_problem *prob, const struct svm_parameter *param, int *status, BlasFunctions *blas_functions); +void svm_csr_cross_validation(const struct svm_csr_problem *prob, const struct svm_parameter *param, int nr_fold, double *target, BlasFunctions *blas_functions); int svm_csr_get_svm_type(const struct svm_csr_model *model); int svm_csr_get_nr_class(const struct svm_csr_model *model); void svm_csr_get_labels(const struct svm_csr_model *model, int *label); double svm_csr_get_svr_probability(const struct svm_csr_model *model); -double svm_csr_predict_values(const struct svm_csr_model *model, const struct svm_csr_node *x, double* dec_values); -double svm_csr_predict(const struct svm_csr_model *model, const struct svm_csr_node *x); -double svm_csr_predict_probability(const struct svm_csr_model *model, const struct svm_csr_node *x, double* prob_estimates); +double svm_csr_predict_values(const struct svm_csr_model *model, const struct svm_csr_node *x, double* dec_values, BlasFunctions *blas_functions); +double svm_csr_predict(const struct svm_csr_model *model, const struct svm_csr_node *x, BlasFunctions *blas_functions); +double svm_csr_predict_probability(const struct svm_csr_model *model, const struct svm_csr_node *x, double* prob_estimates, BlasFunctions *blas_functions); void svm_csr_free_model_content(struct svm_csr_model *model_ptr); void svm_csr_free_and_destroy_model(struct svm_csr_model **model_ptr_ptr);