From 5ceb4e711109fd7802e3e1d0cdff3589221efc57 Mon Sep 17 00:00:00 2001 From: Richard Townsend Date: Sun, 6 Sep 2020 10:01:07 +0100 Subject: [PATCH] linear_models: fix cgo issues, upgrade to liblinear 2.14 Requires an additional step to install: - cd /tmp && - wget https://github.com/cjlin1/liblinear/archive/v241.tar.gz - tar xvf v241.tar.gz - cd liblinear-241 - make lib - sudo install -vm644 linear.h /usr/include - sudo install -vm755 liblinear.so.4 /usr/lib - sudo ln -sfv liblinear.so.4 /usr/lib/liblinear.so --- .travis.yml | 7 +- linear_models/blas.h | 25 - linear_models/blasp.h | 430 ---- linear_models/daxpy.c | 49 - linear_models/ddot.c | 50 - linear_models/dnrm2.c | 62 - linear_models/doc.go | 5 + linear_models/dscal.c | 44 - linear_models/integration.cpp | 109 + linear_models/integration.h | 24 + linear_models/liblinear.go | 176 +- linear_models/linear.cpp | 2837 --------------------------- linear_models/linear.h | 21 +- linear_models/linear_models_test.go | 7 +- linear_models/logistic.go | 1 + linear_models/logistic_test.go | 9 +- linear_models/test.csv | 4 +- linear_models/tmp | Bin 484 -> 0 bytes linear_models/train.csv | 25 +- 19 files changed, 308 insertions(+), 3577 deletions(-) delete mode 100644 linear_models/blas.h delete mode 100644 linear_models/blasp.h delete mode 100644 linear_models/daxpy.c delete mode 100644 linear_models/ddot.c delete mode 100644 linear_models/dnrm2.c delete mode 100644 linear_models/dscal.c create mode 100644 linear_models/integration.cpp create mode 100644 linear_models/integration.h delete mode 100644 linear_models/linear.cpp delete mode 100644 linear_models/tmp diff --git a/.travis.yml b/.travis.yml index 72f9346..4c48d3e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,13 +2,14 @@ language: go go: - 1.13.x - 1.14.x +arch: + - amd64 + - arm64 env: - # Temporary workaround for Go 1.6+ - - GODEBUG=cgocheck=0 before_install: - sudo apt-get update -qq - sudo apt-get install -qq libatlas-base-dev - - cd /tmp && wget http://www.csie.ntu.edu.tw/~cjlin/liblinear/oldfiles/liblinear-1.94.tar.gz && tar xf liblinear-1.94.tar.gz && cd liblinear-1.94 && make lib && sudo install -vm644 linear.h /usr/include && sudo install -vm755 liblinear.so.1 /usr/lib && sudo ln -sfv liblinear.so.1 /usr/lib/liblinear.so + - cd /tmp && wget https://github.com/cjlin1/liblinear/archive/v241.tar.gz && tar xvf v241.tar.gz && cd liblinear-241 && make lib && sudo install -vm644 linear.h /usr/include && sudo install -vm755 liblinear.so.4 /usr/lib && sudo ln -sfv liblinear.so.4 /usr/lib/liblinear.so - cd $TRAVIS_BUILD_DIR install: - go get github.com/smartystreets/goconvey/convey diff --git a/linear_models/blas.h b/linear_models/blas.h deleted file mode 100644 index 558893a..0000000 --- a/linear_models/blas.h +++ /dev/null @@ -1,25 +0,0 @@ -/* blas.h -- C header file for BLAS Ver 1.0 */ -/* Jesse Bennett March 23, 2000 */ - -/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." - - - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ - -#ifndef BLAS_INCLUDE -#define BLAS_INCLUDE - -/* Data types specific to BLAS implementation */ -typedef struct { float r, i; } fcomplex; -typedef struct { double r, i; } dcomplex; -typedef int blasbool; - -#include "blasp.h" /* Prototypes for all BLAS functions */ - -#define FALSE 0 -#define TRUE 1 - -/* Macro functions */ -#define MIN(a,b) ((a) <= (b) ? (a) : (b)) -#define MAX(a,b) ((a) >= (b) ? (a) : (b)) - -#endif diff --git a/linear_models/blasp.h b/linear_models/blasp.h deleted file mode 100644 index 745836d..0000000 --- a/linear_models/blasp.h +++ /dev/null @@ -1,430 +0,0 @@ -/* blasp.h -- C prototypes for BLAS Ver 1.0 */ -/* Jesse Bennett March 23, 2000 */ - -/* Functions listed in alphabetical order */ - -#ifdef F2C_COMPAT - -void cdotc_(fcomplex *dotval, int *n, fcomplex *cx, int *incx, - fcomplex *cy, int *incy); - -void cdotu_(fcomplex *dotval, int *n, fcomplex *cx, int *incx, - fcomplex *cy, int *incy); - -double sasum_(int *n, float *sx, int *incx); - -double scasum_(int *n, fcomplex *cx, int *incx); - -double scnrm2_(int *n, fcomplex *x, int *incx); - -double sdot_(int *n, float *sx, int *incx, float *sy, int *incy); - -double snrm2_(int *n, float *x, int *incx); - -void zdotc_(dcomplex *dotval, int *n, dcomplex *cx, int *incx, - dcomplex *cy, int *incy); - -void zdotu_(dcomplex *dotval, int *n, dcomplex *cx, int *incx, - dcomplex *cy, int *incy); - -#else - -fcomplex cdotc_(int *n, fcomplex *cx, int *incx, fcomplex *cy, int *incy); - -fcomplex cdotu_(int *n, fcomplex *cx, int *incx, fcomplex *cy, int *incy); - -float sasum_(int *n, float *sx, int *incx); - -float scasum_(int *n, fcomplex *cx, int *incx); - -float scnrm2_(int *n, fcomplex *x, int *incx); - -float sdot_(int *n, float *sx, int *incx, float *sy, int *incy); - -float snrm2_(int *n, float *x, int *incx); - -dcomplex zdotc_(int *n, dcomplex *cx, int *incx, dcomplex *cy, int *incy); - -dcomplex zdotu_(int *n, dcomplex *cx, int *incx, dcomplex *cy, int *incy); - -#endif - -/* Remaining functions listed in alphabetical order */ - -int caxpy_(int *n, fcomplex *ca, fcomplex *cx, int *incx, fcomplex *cy, - int *incy); - -int ccopy_(int *n, fcomplex *cx, int *incx, fcomplex *cy, int *incy); - -int cgbmv_(char *trans, int *m, int *n, int *kl, int *ku, - fcomplex *alpha, fcomplex *a, int *lda, fcomplex *x, int *incx, - fcomplex *beta, fcomplex *y, int *incy); - -int cgemm_(char *transa, char *transb, int *m, int *n, int *k, - fcomplex *alpha, fcomplex *a, int *lda, fcomplex *b, int *ldb, - fcomplex *beta, fcomplex *c, int *ldc); - -int cgemv_(char *trans, int *m, int *n, fcomplex *alpha, fcomplex *a, - int *lda, fcomplex *x, int *incx, fcomplex *beta, fcomplex *y, - int *incy); - -int cgerc_(int *m, int *n, fcomplex *alpha, fcomplex *x, int *incx, - fcomplex *y, int *incy, fcomplex *a, int *lda); - -int cgeru_(int *m, int *n, fcomplex *alpha, fcomplex *x, int *incx, - fcomplex *y, int *incy, fcomplex *a, int *lda); - -int chbmv_(char *uplo, int *n, int *k, fcomplex *alpha, fcomplex *a, - int *lda, fcomplex *x, int *incx, fcomplex *beta, fcomplex *y, - int *incy); - -int chemm_(char *side, char *uplo, int *m, int *n, fcomplex *alpha, - fcomplex *a, int *lda, fcomplex *b, int *ldb, fcomplex *beta, - fcomplex *c, int *ldc); - -int chemv_(char *uplo, int *n, fcomplex *alpha, fcomplex *a, int *lda, - fcomplex *x, int *incx, fcomplex *beta, fcomplex *y, int *incy); - -int cher_(char *uplo, int *n, float *alpha, fcomplex *x, int *incx, - fcomplex *a, int *lda); - -int cher2_(char *uplo, int *n, fcomplex *alpha, fcomplex *x, int *incx, - fcomplex *y, int *incy, fcomplex *a, int *lda); - -int cher2k_(char *uplo, char *trans, int *n, int *k, fcomplex *alpha, - fcomplex *a, int *lda, fcomplex *b, int *ldb, float *beta, - fcomplex *c, int *ldc); - -int cherk_(char *uplo, char *trans, int *n, int *k, float *alpha, - fcomplex *a, int *lda, float *beta, fcomplex *c, int *ldc); - -int chpmv_(char *uplo, int *n, fcomplex *alpha, fcomplex *ap, fcomplex *x, - int *incx, fcomplex *beta, fcomplex *y, int *incy); - -int chpr_(char *uplo, int *n, float *alpha, fcomplex *x, int *incx, - fcomplex *ap); - -int chpr2_(char *uplo, int *n, fcomplex *alpha, fcomplex *x, int *incx, - fcomplex *y, int *incy, fcomplex *ap); - -int crotg_(fcomplex *ca, fcomplex *cb, float *c, fcomplex *s); - -int cscal_(int *n, fcomplex *ca, fcomplex *cx, int *incx); - -int csscal_(int *n, float *sa, fcomplex *cx, int *incx); - -int cswap_(int *n, fcomplex *cx, int *incx, fcomplex *cy, int *incy); - -int csymm_(char *side, char *uplo, int *m, int *n, fcomplex *alpha, - fcomplex *a, int *lda, fcomplex *b, int *ldb, fcomplex *beta, - fcomplex *c, int *ldc); - -int csyr2k_(char *uplo, char *trans, int *n, int *k, fcomplex *alpha, - fcomplex *a, int *lda, fcomplex *b, int *ldb, fcomplex *beta, - fcomplex *c, int *ldc); - -int csyrk_(char *uplo, char *trans, int *n, int *k, fcomplex *alpha, - fcomplex *a, int *lda, fcomplex *beta, fcomplex *c, int *ldc); - -int ctbmv_(char *uplo, char *trans, char *diag, int *n, int *k, - fcomplex *a, int *lda, fcomplex *x, int *incx); - -int ctbsv_(char *uplo, char *trans, char *diag, int *n, int *k, - fcomplex *a, int *lda, fcomplex *x, int *incx); - -int ctpmv_(char *uplo, char *trans, char *diag, int *n, fcomplex *ap, - fcomplex *x, int *incx); - -int ctpsv_(char *uplo, char *trans, char *diag, int *n, fcomplex *ap, - fcomplex *x, int *incx); - -int ctrmm_(char *side, char *uplo, char *transa, char *diag, int *m, - int *n, fcomplex *alpha, fcomplex *a, int *lda, fcomplex *b, - int *ldb); - -int ctrmv_(char *uplo, char *trans, char *diag, int *n, fcomplex *a, - int *lda, fcomplex *x, int *incx); - -int ctrsm_(char *side, char *uplo, char *transa, char *diag, int *m, - int *n, fcomplex *alpha, fcomplex *a, int *lda, fcomplex *b, - int *ldb); - -int ctrsv_(char *uplo, char *trans, char *diag, int *n, fcomplex *a, - int *lda, fcomplex *x, int *incx); - -int daxpy_(int *n, double *sa, double *sx, int *incx, double *sy, - int *incy); - -int dcopy_(int *n, double *sx, int *incx, double *sy, int *incy); - -int dgbmv_(char *trans, int *m, int *n, int *kl, int *ku, - double *alpha, double *a, int *lda, double *x, int *incx, - double *beta, double *y, int *incy); - -int dgemm_(char *transa, char *transb, int *m, int *n, int *k, - double *alpha, double *a, int *lda, double *b, int *ldb, - double *beta, double *c, int *ldc); - -int dgemv_(char *trans, int *m, int *n, double *alpha, double *a, - int *lda, double *x, int *incx, double *beta, double *y, - int *incy); - -int dger_(int *m, int *n, double *alpha, double *x, int *incx, - double *y, int *incy, double *a, int *lda); - -int drot_(int *n, double *sx, int *incx, double *sy, int *incy, - double *c, double *s); - -int drotg_(double *sa, double *sb, double *c, double *s); - -int dsbmv_(char *uplo, int *n, int *k, double *alpha, double *a, - int *lda, double *x, int *incx, double *beta, double *y, - int *incy); - -int dscal_(int *n, double *sa, double *sx, int *incx); - -int dspmv_(char *uplo, int *n, double *alpha, double *ap, double *x, - int *incx, double *beta, double *y, int *incy); - -int dspr_(char *uplo, int *n, double *alpha, double *x, int *incx, - double *ap); - -int dspr2_(char *uplo, int *n, double *alpha, double *x, int *incx, - double *y, int *incy, double *ap); - -int dswap_(int *n, double *sx, int *incx, double *sy, int *incy); - -int dsymm_(char *side, char *uplo, int *m, int *n, double *alpha, - double *a, int *lda, double *b, int *ldb, double *beta, - double *c, int *ldc); - -int dsymv_(char *uplo, int *n, double *alpha, double *a, int *lda, - double *x, int *incx, double *beta, double *y, int *incy); - -int dsyr_(char *uplo, int *n, double *alpha, double *x, int *incx, - double *a, int *lda); - -int dsyr2_(char *uplo, int *n, double *alpha, double *x, int *incx, - double *y, int *incy, double *a, int *lda); - -int dsyr2k_(char *uplo, char *trans, int *n, int *k, double *alpha, - double *a, int *lda, double *b, int *ldb, double *beta, - double *c, int *ldc); - -int dsyrk_(char *uplo, char *trans, int *n, int *k, double *alpha, - double *a, int *lda, double *beta, double *c, int *ldc); - -int dtbmv_(char *uplo, char *trans, char *diag, int *n, int *k, - double *a, int *lda, double *x, int *incx); - -int dtbsv_(char *uplo, char *trans, char *diag, int *n, int *k, - double *a, int *lda, double *x, int *incx); - -int dtpmv_(char *uplo, char *trans, char *diag, int *n, double *ap, - double *x, int *incx); - -int dtpsv_(char *uplo, char *trans, char *diag, int *n, double *ap, - double *x, int *incx); - -int dtrmm_(char *side, char *uplo, char *transa, char *diag, int *m, - int *n, double *alpha, double *a, int *lda, double *b, - int *ldb); - -int dtrmv_(char *uplo, char *trans, char *diag, int *n, double *a, - int *lda, double *x, int *incx); - -int dtrsm_(char *side, char *uplo, char *transa, char *diag, int *m, - int *n, double *alpha, double *a, int *lda, double *b, - int *ldb); - -int dtrsv_(char *uplo, char *trans, char *diag, int *n, double *a, - int *lda, double *x, int *incx); - - -int saxpy_(int *n, float *sa, float *sx, int *incx, float *sy, int *incy); - -int scopy_(int *n, float *sx, int *incx, float *sy, int *incy); - -int sgbmv_(char *trans, int *m, int *n, int *kl, int *ku, - float *alpha, float *a, int *lda, float *x, int *incx, - float *beta, float *y, int *incy); - -int sgemm_(char *transa, char *transb, int *m, int *n, int *k, - float *alpha, float *a, int *lda, float *b, int *ldb, - float *beta, float *c, int *ldc); - -int sgemv_(char *trans, int *m, int *n, float *alpha, float *a, - int *lda, float *x, int *incx, float *beta, float *y, - int *incy); - -int sger_(int *m, int *n, float *alpha, float *x, int *incx, - float *y, int *incy, float *a, int *lda); - -int srot_(int *n, float *sx, int *incx, float *sy, int *incy, - float *c, float *s); - -int srotg_(float *sa, float *sb, float *c, float *s); - -int ssbmv_(char *uplo, int *n, int *k, float *alpha, float *a, - int *lda, float *x, int *incx, float *beta, float *y, - int *incy); - -int sscal_(int *n, float *sa, float *sx, int *incx); - -int sspmv_(char *uplo, int *n, float *alpha, float *ap, float *x, - int *incx, float *beta, float *y, int *incy); - -int sspr_(char *uplo, int *n, float *alpha, float *x, int *incx, - float *ap); - -int sspr2_(char *uplo, int *n, float *alpha, float *x, int *incx, - float *y, int *incy, float *ap); - -int sswap_(int *n, float *sx, int *incx, float *sy, int *incy); - -int ssymm_(char *side, char *uplo, int *m, int *n, float *alpha, - float *a, int *lda, float *b, int *ldb, float *beta, - float *c, int *ldc); - -int ssymv_(char *uplo, int *n, float *alpha, float *a, int *lda, - float *x, int *incx, float *beta, float *y, int *incy); - -int ssyr_(char *uplo, int *n, float *alpha, float *x, int *incx, - float *a, int *lda); - -int ssyr2_(char *uplo, int *n, float *alpha, float *x, int *incx, - float *y, int *incy, float *a, int *lda); - -int ssyr2k_(char *uplo, char *trans, int *n, int *k, float *alpha, - float *a, int *lda, float *b, int *ldb, float *beta, - float *c, int *ldc); - -int ssyrk_(char *uplo, char *trans, int *n, int *k, float *alpha, - float *a, int *lda, float *beta, float *c, int *ldc); - -int stbmv_(char *uplo, char *trans, char *diag, int *n, int *k, - float *a, int *lda, float *x, int *incx); - -int stbsv_(char *uplo, char *trans, char *diag, int *n, int *k, - float *a, int *lda, float *x, int *incx); - -int stpmv_(char *uplo, char *trans, char *diag, int *n, float *ap, - float *x, int *incx); - -int stpsv_(char *uplo, char *trans, char *diag, int *n, float *ap, - float *x, int *incx); - -int strmm_(char *side, char *uplo, char *transa, char *diag, int *m, - int *n, float *alpha, float *a, int *lda, float *b, - int *ldb); - -int strmv_(char *uplo, char *trans, char *diag, int *n, float *a, - int *lda, float *x, int *incx); - -int strsm_(char *side, char *uplo, char *transa, char *diag, int *m, - int *n, float *alpha, float *a, int *lda, float *b, - int *ldb); - -int strsv_(char *uplo, char *trans, char *diag, int *n, float *a, - int *lda, float *x, int *incx); - -int zaxpy_(int *n, dcomplex *ca, dcomplex *cx, int *incx, dcomplex *cy, - int *incy); - -int zcopy_(int *n, dcomplex *cx, int *incx, dcomplex *cy, int *incy); - -int zdscal_(int *n, double *sa, dcomplex *cx, int *incx); - -int zgbmv_(char *trans, int *m, int *n, int *kl, int *ku, - dcomplex *alpha, dcomplex *a, int *lda, dcomplex *x, int *incx, - dcomplex *beta, dcomplex *y, int *incy); - -int zgemm_(char *transa, char *transb, int *m, int *n, int *k, - dcomplex *alpha, dcomplex *a, int *lda, dcomplex *b, int *ldb, - dcomplex *beta, dcomplex *c, int *ldc); - -int zgemv_(char *trans, int *m, int *n, dcomplex *alpha, dcomplex *a, - int *lda, dcomplex *x, int *incx, dcomplex *beta, dcomplex *y, - int *incy); - -int zgerc_(int *m, int *n, dcomplex *alpha, dcomplex *x, int *incx, - dcomplex *y, int *incy, dcomplex *a, int *lda); - -int zgeru_(int *m, int *n, dcomplex *alpha, dcomplex *x, int *incx, - dcomplex *y, int *incy, dcomplex *a, int *lda); - -int zhbmv_(char *uplo, int *n, int *k, dcomplex *alpha, dcomplex *a, - int *lda, dcomplex *x, int *incx, dcomplex *beta, dcomplex *y, - int *incy); - -int zhemm_(char *side, char *uplo, int *m, int *n, dcomplex *alpha, - dcomplex *a, int *lda, dcomplex *b, int *ldb, dcomplex *beta, - dcomplex *c, int *ldc); - -int zhemv_(char *uplo, int *n, dcomplex *alpha, dcomplex *a, int *lda, - dcomplex *x, int *incx, dcomplex *beta, dcomplex *y, int *incy); - -int zher_(char *uplo, int *n, double *alpha, dcomplex *x, int *incx, - dcomplex *a, int *lda); - -int zher2_(char *uplo, int *n, dcomplex *alpha, dcomplex *x, int *incx, - dcomplex *y, int *incy, dcomplex *a, int *lda); - -int zher2k_(char *uplo, char *trans, int *n, int *k, dcomplex *alpha, - dcomplex *a, int *lda, dcomplex *b, int *ldb, double *beta, - dcomplex *c, int *ldc); - -int zherk_(char *uplo, char *trans, int *n, int *k, double *alpha, - dcomplex *a, int *lda, double *beta, dcomplex *c, int *ldc); - -int zhpmv_(char *uplo, int *n, dcomplex *alpha, dcomplex *ap, dcomplex *x, - int *incx, dcomplex *beta, dcomplex *y, int *incy); - -int zhpr_(char *uplo, int *n, double *alpha, dcomplex *x, int *incx, - dcomplex *ap); - -int zhpr2_(char *uplo, int *n, dcomplex *alpha, dcomplex *x, int *incx, - dcomplex *y, int *incy, dcomplex *ap); - -int zrotg_(dcomplex *ca, dcomplex *cb, double *c, dcomplex *s); - -int zscal_(int *n, dcomplex *ca, dcomplex *cx, int *incx); - -int zswap_(int *n, dcomplex *cx, int *incx, dcomplex *cy, int *incy); - -int zsymm_(char *side, char *uplo, int *m, int *n, dcomplex *alpha, - dcomplex *a, int *lda, dcomplex *b, int *ldb, dcomplex *beta, - dcomplex *c, int *ldc); - -int zsyr2k_(char *uplo, char *trans, int *n, int *k, dcomplex *alpha, - dcomplex *a, int *lda, dcomplex *b, int *ldb, dcomplex *beta, - dcomplex *c, int *ldc); - -int zsyrk_(char *uplo, char *trans, int *n, int *k, dcomplex *alpha, - dcomplex *a, int *lda, dcomplex *beta, dcomplex *c, int *ldc); - -int ztbmv_(char *uplo, char *trans, char *diag, int *n, int *k, - dcomplex *a, int *lda, dcomplex *x, int *incx); - -int ztbsv_(char *uplo, char *trans, char *diag, int *n, int *k, - dcomplex *a, int *lda, dcomplex *x, int *incx); - -int ztpmv_(char *uplo, char *trans, char *diag, int *n, dcomplex *ap, - dcomplex *x, int *incx); - -int ztpsv_(char *uplo, char *trans, char *diag, int *n, dcomplex *ap, - dcomplex *x, int *incx); - -int ztrmm_(char *side, char *uplo, char *transa, char *diag, int *m, - int *n, dcomplex *alpha, dcomplex *a, int *lda, dcomplex *b, - int *ldb); - -int ztrmv_(char *uplo, char *trans, char *diag, int *n, dcomplex *a, - int *lda, dcomplex *x, int *incx); - -int ztrsm_(char *side, char *uplo, char *transa, char *diag, int *m, - int *n, dcomplex *alpha, dcomplex *a, int *lda, dcomplex *b, - int *ldb); - -int ztrsv_(char *uplo, char *trans, char *diag, int *n, dcomplex *a, - int *lda, dcomplex *x, int *incx); diff --git a/linear_models/daxpy.c b/linear_models/daxpy.c deleted file mode 100644 index 58f345a..0000000 --- a/linear_models/daxpy.c +++ /dev/null @@ -1,49 +0,0 @@ -#include "blas.h" - -int daxpy_(int *n, double *sa, double *sx, int *incx, double *sy, - int *incy) -{ - long int i, m, ix, iy, nn, iincx, iincy; - register double ssa; - - /* constant times a vector plus a vector. - uses unrolled loop for increments equal to one. - jack dongarra, linpack, 3/11/78. - modified 12/3/93, array(1) declarations changed to array(*) */ - - /* Dereference inputs */ - nn = *n; - ssa = *sa; - iincx = *incx; - iincy = *incy; - - if( nn > 0 && ssa != 0.0 ) - { - if (iincx == 1 && iincy == 1) /* code for both increments equal to 1 */ - { - m = nn-3; - for (i = 0; i < m; i += 4) - { - sy[i] += ssa * sx[i]; - sy[i+1] += ssa * sx[i+1]; - sy[i+2] += ssa * sx[i+2]; - sy[i+3] += ssa * sx[i+3]; - } - for ( ; i < nn; ++i) /* clean-up loop */ - sy[i] += ssa * sx[i]; - } - else /* code for unequal increments or equal increments not equal to 1 */ - { - ix = iincx >= 0 ? 0 : (1 - nn) * iincx; - iy = iincy >= 0 ? 0 : (1 - nn) * iincy; - for (i = 0; i < nn; i++) - { - sy[iy] += ssa * sx[ix]; - ix += iincx; - iy += iincy; - } - } - } - - return 0; -} /* daxpy_ */ diff --git a/linear_models/ddot.c b/linear_models/ddot.c deleted file mode 100644 index a64a280..0000000 --- a/linear_models/ddot.c +++ /dev/null @@ -1,50 +0,0 @@ -#include "blas.h" - -double ddot_(int *n, double *sx, int *incx, double *sy, int *incy) -{ - long int i, m, nn, iincx, iincy; - double stemp; - long int ix, iy; - - /* forms the dot product of two vectors. - uses unrolled loops for increments equal to one. - jack dongarra, linpack, 3/11/78. - modified 12/3/93, array(1) declarations changed to array(*) */ - - /* Dereference inputs */ - nn = *n; - iincx = *incx; - iincy = *incy; - - stemp = 0.0; - if (nn > 0) - { - if (iincx == 1 && iincy == 1) /* code for both increments equal to 1 */ - { - m = nn-4; - for (i = 0; i < m; i += 5) - stemp += sx[i] * sy[i] + sx[i+1] * sy[i+1] + sx[i+2] * sy[i+2] + - sx[i+3] * sy[i+3] + sx[i+4] * sy[i+4]; - - for ( ; i < nn; i++) /* clean-up loop */ - stemp += sx[i] * sy[i]; - } - else /* code for unequal increments or equal increments not equal to 1 */ - { - ix = 0; - iy = 0; - if (iincx < 0) - ix = (1 - nn) * iincx; - if (iincy < 0) - iy = (1 - nn) * iincy; - for (i = 0; i < nn; i++) - { - stemp += sx[ix] * sy[iy]; - ix += iincx; - iy += iincy; - } - } - } - - return stemp; -} /* ddot_ */ diff --git a/linear_models/dnrm2.c b/linear_models/dnrm2.c deleted file mode 100644 index e50cdf7..0000000 --- a/linear_models/dnrm2.c +++ /dev/null @@ -1,62 +0,0 @@ -#include /* Needed for fabs() and sqrt() */ -#include "blas.h" - -double dnrm2_(int *n, double *x, int *incx) -{ - long int ix, nn, iincx; - double norm, scale, absxi, ssq, temp; - -/* DNRM2 returns the euclidean norm of a vector via the function - name, so that - - DNRM2 := sqrt( x'*x ) - - -- This version written on 25-October-1982. - Modified on 14-October-1993 to inline the call to SLASSQ. - Sven Hammarling, Nag Ltd. */ - - /* Dereference inputs */ - nn = *n; - iincx = *incx; - - if( nn > 0 && iincx > 0 ) - { - if (nn == 1) - { - norm = fabs(x[0]); - } - else - { - scale = 0.0; - ssq = 1.0; - - /* The following loop is equivalent to this call to the LAPACK - auxiliary routine: CALL SLASSQ( N, X, INCX, SCALE, SSQ ) */ - - for (ix=(nn-1)*iincx; ix>=0; ix-=iincx) - { - if (x[ix] != 0.0) - { - absxi = fabs(x[ix]); - if (scale < absxi) - { - temp = scale / absxi; - ssq = ssq * (temp * temp) + 1.0; - scale = absxi; - } - else - { - temp = absxi / scale; - ssq += temp * temp; - } - } - } - norm = scale * sqrt(ssq); - } - } - else - norm = 0.0; - - return norm; - -} /* dnrm2_ */ diff --git a/linear_models/doc.go b/linear_models/doc.go index 53b040e..56c4e2f 100644 --- a/linear_models/doc.go +++ b/linear_models/doc.go @@ -1,5 +1,10 @@ /* Package linear_models implements linear and logistic regression models. + +train.csv, test.csv are truncated versions of +https://www.kaggle.com/uciml/pima-indians-diabetes-database + + */ package linear_models diff --git a/linear_models/dscal.c b/linear_models/dscal.c deleted file mode 100644 index a0eca0c..0000000 --- a/linear_models/dscal.c +++ /dev/null @@ -1,44 +0,0 @@ -#include "blas.h" - -int dscal_(int *n, double *sa, double *sx, int *incx) -{ - long int i, m, nincx, nn, iincx; - double ssa; - - /* scales a vector by a constant. - uses unrolled loops for increment equal to 1. - jack dongarra, linpack, 3/11/78. - modified 3/93 to return if incx .le. 0. - modified 12/3/93, array(1) declarations changed to array(*) */ - - /* Dereference inputs */ - nn = *n; - iincx = *incx; - ssa = *sa; - - if (nn > 0 && iincx > 0) - { - if (iincx == 1) /* code for increment equal to 1 */ - { - m = nn-4; - for (i = 0; i < m; i += 5) - { - sx[i] = ssa * sx[i]; - sx[i+1] = ssa * sx[i+1]; - sx[i+2] = ssa * sx[i+2]; - sx[i+3] = ssa * sx[i+3]; - sx[i+4] = ssa * sx[i+4]; - } - for ( ; i < nn; ++i) /* clean-up loop */ - sx[i] = ssa * sx[i]; - } - else /* code for increment not equal to 1 */ - { - nincx = nn * iincx; - for (i = 0; i < nincx; i += iincx) - sx[i] = ssa * sx[i]; - } - } - - return 0; -} /* dscal_ */ diff --git a/linear_models/integration.cpp b/linear_models/integration.cpp new file mode 100644 index 0000000..f973875 --- /dev/null +++ b/linear_models/integration.cpp @@ -0,0 +1,109 @@ +/* + * This file contains functions related to creating + freeing + * objects on behalf of the go runtime + */ + +#include "linear.h" +#include + +extern "C" { + +/* NOTE: the Golang versions of the structures must call the corresponding + * Free functions via runtime.SetFinalize */ +/* CreateCProblem allocates a new struct problem outside of Golang's + * garbage collection. */ +struct problem *CreateCProblem() { + auto ret = new problem(); + *ret = {}; // < Clear all fields + return ret; +} + +/* CreateCModel allocates a new struct model outside of Golang's + * garbage collection. */ +struct model *CreateCModel() { + auto ret = new model(); + *ret = {}; // < Clear all fields + return ret; +} + +/* CreateCParameter allocates a new struct parameter outside of + * Golang's garbage collection.*/ +struct parameter *CreateCParameter() { + auto ret = new parameter(); + *ret = {}; + return ret; +} + +/* Free's a previously allocated problem and all its data */ +void FreeCProblem(struct problem *p) { + if (p->y != nullptr) { + free(p->y); + p->y = nullptr; + } + if (p->x != nullptr) { + // l is the total count of rows in the problem + // n is the number of values in each row + for (int i = 0; i < p->l; i++) { + if (p->x[i] != nullptr) { + free(p->x[i]); + p->x[i] = nullptr; + } + } + free(p->x); + p->x = nullptr; + } + delete p; +} + +/* free's a model with libsvm's internal routines */ +void FreeCModel(struct model *m) { + free_model_content(m); + delete m; +} + +/* free's a parameter via libsvm */ +void FreeCParameter(struct parameter *p) { + destroy_param(p); + delete p; +} + +/* Allocates a vector of doubles for storing target values + * outside of Go's garbage collection */ +int AllocateLabelsForProblem (struct problem *p, int numValues) { + p->y = reinterpret_cast(malloc(sizeof(double) * numValues)); + return p->y == nullptr; +} + +/* Utility method used to set the target value for a particular + * input row */ +void AssignLabelForProblem(struct problem *p, int i, double d) { + p->y[i] = d; +} + +/* Returns a feature node for a particular row and column. */ +struct feature_node *GetFeatureNodeForIndex(struct problem *p, int i, int j) { + return &(p->x[i][j]); +} + +/* Allocates a buffer of input rows and the values to fill them. */ +int AllocateFeatureNodesForProblem(struct problem *p, + int numSamples, int numValues, int *rowLengths) { + + numValues++; // Extend for terminating element + p->x = reinterpret_cast( + calloc(numSamples, sizeof(struct feature_node *)) + ); + if (p->x == nullptr) { + return -1; + } + // Allocate all of the feature nodes, then riffle them together + struct feature_node *feature_nodes = reinterpret_cast(calloc(numValues, sizeof(struct feature_node))); + + int cur = 0; + for (int i = 0; i < numSamples; cur += rowLengths[i], i++) { + p->x[i] = feature_nodes + cur; + } + return 0; +} + +} /* extern "C" */ diff --git a/linear_models/integration.h b/linear_models/integration.h new file mode 100644 index 0000000..5c28a32 --- /dev/null +++ b/linear_models/integration.h @@ -0,0 +1,24 @@ +#ifndef _H_INTEGRATION_ +#define _H_INTEGRATION_ + +#include "linear.h" + +struct problem *CreateCProblem(); +void FreeCProblem(struct problem*); +struct model *CreateCModel(); +void FreeCModel(struct model*); +struct parameter *CreateCParameter(); +void FreeCParameter(struct parameter*); +// Allocates memory outside of golang for describing feature +// vectors. First pointer is the C-managed liblinear problem, +// second parameter is the number of rows we're training on +// third parameter is the total number of non-zero elements +// (including bias and null terminators) that we need to allocate +// and final parameter is an array describing the number of +// nodes in each row. +int AllocateFeatureNodesForProblem(struct problem*, int, int, int*); +int AllocateLabelsForProblem(struct problem *, int); +void AssignLabelForProblem(struct problem *, int, double); +struct feature_node *GetFeatureNodeForIndex(struct problem *, int, int); + +#endif diff --git a/linear_models/liblinear.go b/linear_models/liblinear.go index 47e7456..4d542b3 100644 --- a/linear_models/liblinear.go +++ b/linear_models/liblinear.go @@ -1,22 +1,47 @@ package linear_models /* -#include "linear.h" +#include "integration.h" +#cgo CFLAGS: +#cgo CXXFLAGS: -std=c++11 -g -O0 +#cgo LDFLAGS: -g -llinear */ import "C" -import "fmt" -import "unsafe" +import ( + "fmt" + "runtime" +) +// Problem wraps a libsvm problem struct which describes a classification/ +// regression problem. No externally-accessible fields. type Problem struct { - c_prob C.struct_problem + c_prob *C.struct_problem } +// Free releases resources associated with a libsvm problem. +func (p *Problem) Free() { + C.FreeCProblem(p.c_prob) +} + +// Parameter encasulates all the possible libsvm training options. +// TODO: make user control of these more extensive. type Parameter struct { - c_param C.struct_parameter + c_param *C.struct_parameter } +// Free releases resources associated with a Parameter. +func (p *Parameter) Free() { + C.FreeCParameter(p.c_param) +} + +// Model encapsulates a trained libsvm model. type Model struct { - c_model unsafe.Pointer + c_model *C.struct_model +} + +// Free releases resources associated with a trained libsvm model. +func (m *Model) Free() { + C.FreeCModel(m.c_model) } const ( @@ -30,8 +55,14 @@ const ( L2R_LR_DUAL = C.L2R_LR_DUAL ) +// NewParameter creates a libsvm parameter structure, which controls +// various aspects of libsvm training. +// For more information on what these parameters do, consult the +// "`train` usage" section of +// https://github.com/cjlin1/liblinear/blob/master/README func NewParameter(solver_type int, C float64, eps float64) *Parameter { - param := Parameter{} + param := &Parameter{C.CreateCParameter()} + runtime.SetFinalizer(param, (*Parameter).Free) param.c_param.solver_type = C.int(solver_type) param.c_param.eps = C.double(eps) param.c_param.C = C.double(C) @@ -39,30 +70,37 @@ func NewParameter(solver_type int, C float64, eps float64) *Parameter { param.c_param.weight_label = nil param.c_param.weight = nil - return ¶m + return param } +// NewProblem creates input to libsvm which describes a particular +// regression/classification problem. It requires an array of float values +// and an array of y values. func NewProblem(X [][]float64, y []float64, bias float64) *Problem { - prob := Problem{} + prob := &Problem{C.CreateCProblem()} + runtime.SetFinalizer(prob, (*Problem).Free) prob.c_prob.l = C.int(len(X)) prob.c_prob.n = C.int(len(X[0]) + 1) - prob.c_prob.x = convert_features(X, bias) - c_y := make([]C.double, len(y)) + convert_features(prob, X, bias) + C.AllocateLabelsForProblem(prob.c_prob, C.int(len(y))) for i := 0; i < len(y); i++ { - c_y[i] = C.double(y[i]) + C.AssignLabelForProblem(prob.c_prob, C.int(i), C.double(y[i])) } - prob.c_prob.y = &c_y[0] + // Should not go out of scope until the Problem struct + // is cleaned up. prob.c_prob.bias = C.double(-1) - return &prob + return prob } +// Train invokes libsvm and returns a trained model. func Train(prob *Problem, param *Parameter) *Model { libLinearHookPrintFunc() // Sets up logging - tmpCProb := &prob.c_prob - tmpCParam := ¶m.c_param - return &Model{unsafe.Pointer(C.train(tmpCProb, tmpCParam))} + out := C.train(prob.c_prob, param.c_param) + m := &Model{out} + runtime.SetFinalizer(m, (*Model).Free) + return m } func Export(model *Model, filePath string) error { @@ -74,82 +112,100 @@ func Export(model *Model, filePath string) error { } func Load(model *Model, filePath string) error { - model.c_model = unsafe.Pointer(C.load_model(C.CString(filePath))) + model.c_model = C.load_model(C.CString(filePath)) if model.c_model == nil { return fmt.Errorf("Something went wrong") } return nil + } +// Predict takes a row of float values corresponding to a particular +// input and returns the regression result. func Predict(model *Model, x []float64) float64 { - c_x := convert_vector(x, 0) - c_y := C.predict((*C.struct_model)(model.c_model), c_x) - y := float64(c_y) + cX := convertVector(x, 0) + cY := C.predict((*C.struct_model)(model.c_model), &cX[0]) + y := float64(cY) return y } -func convert_vector(x []float64, bias float64) *C.struct_feature_node { - n_ele := 0 + +// convertVector is an internal function used for converting +// dense float64 vectors into the sparse input that libsvm accepts. +func convertVector(x []float64, bias float64) []C.struct_feature_node { + // Count the number of non-zero elements + nElements := 0 for i := 0; i < len(x); i++ { if x[i] > 0 { - n_ele++ + nElements++ } } - n_ele += 2 + // Add one at the end for the -1 terminator + nElements++ + if bias >= 0 { + // And one for the bias, if we have it + nElements++ + } - c_x := make([]C.struct_feature_node, n_ele) + cX := make([]C.struct_feature_node, nElements) j := 0 for i := 0; i < len(x); i++ { if x[i] > 0 { - c_x[j].index = C.int(i + 1) - c_x[j].value = C.double(x[i]) + cX[j].index = C.int(i + 1) + cX[j].value = C.double(x[i]) j++ } } - if bias > 0 { - c_x[j].index = C.int(0) - c_x[j].value = C.double(0) + if bias >= 0 { + cX[j].index = C.int(0) + cX[j].value = C.double(0) j++ } - c_x[j].index = C.int(-1) - return &c_x[0] + cX[j].index = C.int(-1) + return cX } -func convert_features(X [][]float64, bias float64) **C.struct_feature_node { - n_samples := len(X) - n_elements := 0 - for i := 0; i < n_samples; i++ { +// convert_features is an internal function used for converting +// dense 2D arrays of float values into the sparse format libsvm accepts. +func convert_features(prob *Problem, X [][]float64, bias float64) { + + nonZeroRowElements := make([]C.int, len(X)) + totalElements := 0 + + for i := 0; i < len(X); i++ { + // For each row of input data, we count how many non-zero things are in the row + nonZeroElementsInRow := 1 // Initially one, because we need the -1 null terminator for j := 0; j < len(X[i]); j++ { if X[i][j] != 0.0 { - n_elements++ + nonZeroElementsInRow++ + } + if bias >= 0 { + nonZeroElementsInRow++ } - n_elements++ //for bias } + nonZeroRowElements[i] = C.int(nonZeroElementsInRow) + totalElements += nonZeroElementsInRow } + // Allocate one feature vector for each row, total number + C.AllocateFeatureNodesForProblem(prob.c_prob, C.int(len(X)), C.int(totalElements), &nonZeroRowElements[0]) - x_space := make([]C.struct_feature_node, n_elements+n_samples) - - cursor := 0 - x := make([]*C.struct_feature_node, n_samples) - var c_x **C.struct_feature_node - - for i := 0; i < n_samples; i++ { - x[i] = &x_space[cursor] - + for i := 0; i < len(X); i++ { + nonZeroElementCounter := 0 for j := 0; j < len(X[i]); j++ { if X[i][j] != 0.0 { - x_space[cursor].index = C.int(j + 1) - x_space[cursor].value = C.double(X[i][j]) - cursor++ - } - if bias > 0 { - x_space[cursor].index = C.int(0) - x_space[cursor].value = C.double(bias) - cursor++ + xSpace := C.GetFeatureNodeForIndex(prob.c_prob, C.int(i), C.int(nonZeroElementCounter)) + xSpace.index = C.int(j + 1) + xSpace.value = C.double(X[i][j]) + nonZeroElementCounter++ } } - x_space[cursor].index = C.int(-1) - cursor++ + if bias >= 0 { + xSpace := C.GetFeatureNodeForIndex(prob.c_prob, C.int(i), C.int(nonZeroElementCounter)) + xSpace.index = C.int(len(X[i]) + 1) + xSpace.value = C.double(bias) + nonZeroElementCounter++ + } + xSpace := C.GetFeatureNodeForIndex(prob.c_prob, C.int(i), C.int(nonZeroElementCounter)) + xSpace.index = C.int(-1) + xSpace.value = C.double(0) } - c_x = &x[0] - return c_x } diff --git a/linear_models/linear.cpp b/linear_models/linear.cpp deleted file mode 100644 index e435db7..0000000 --- a/linear_models/linear.cpp +++ /dev/null @@ -1,2837 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include "linear.h" -#include "tron.h" -typedef signed char schar; -template static inline void swap(T& x, T& y) { T t=x; x=y; y=t; } -#ifndef min -template static inline T min(T x,T y) { return (x static inline T max(T x,T y) { return (x>y)?x:y; } -#endif -template static inline void clone(T*& dst, S* src, int n) -{ - dst = new T[n]; - memcpy((void *)dst,(void *)src,sizeof(T)*n); -} -#define Malloc(type,n) (type *)malloc((n)*sizeof(type)) -#define INF HUGE_VAL - -static void print_string_stdout(const char *s) -{ - fputs(s,stdout); - fflush(stdout); -} - -static void (*liblinear_print_string) (const char *) = &print_string_stdout; - -#if 1 -static void info(const char *fmt,...) -{ - char buf[BUFSIZ]; - va_list ap; - va_start(ap,fmt); - vsprintf(buf,fmt,ap); - va_end(ap); - (*liblinear_print_string)(buf); -} -#else -static void info(const char *fmt,...) {} -#endif - -class l2r_lr_fun: public function -{ -public: - l2r_lr_fun(const problem *prob, double *C); - ~l2r_lr_fun(); - - double fun(double *w); - void grad(double *w, double *g); - void Hv(double *s, double *Hs); - - int get_nr_variable(void); - -private: - void Xv(double *v, double *Xv); - void XTv(double *v, double *XTv); - - double *C; - double *z; - double *D; - const problem *prob; -}; - -l2r_lr_fun::l2r_lr_fun(const problem *prob, double *C) -{ - int l=prob->l; - - this->prob = prob; - - z = new double[l]; - D = new double[l]; - this->C = C; -} - -l2r_lr_fun::~l2r_lr_fun() -{ - delete[] z; - delete[] D; -} - - -double l2r_lr_fun::fun(double *w) -{ - int i; - double f=0; - double *y=prob->y; - int l=prob->l; - int w_size=get_nr_variable(); - - Xv(w, z); - - for(i=0;i= 0) - f += C[i]*log(1 + exp(-yz)); - else - f += C[i]*(-yz+log(1 + exp(yz))); - } - - return(f); -} - -void l2r_lr_fun::grad(double *w, double *g) -{ - int i; - double *y=prob->y; - int l=prob->l; - int w_size=get_nr_variable(); - - for(i=0;in; -} - -void l2r_lr_fun::Hv(double *s, double *Hs) -{ - int i; - int l=prob->l; - int w_size=get_nr_variable(); - double *wa = new double[l]; - - Xv(s, wa); - for(i=0;il; - feature_node **x=prob->x; - - for(i=0;iindex!=-1) - { - Xv[i]+=v[s->index-1]*s->value; - s++; - } - } -} - -void l2r_lr_fun::XTv(double *v, double *XTv) -{ - int i; - int l=prob->l; - int w_size=get_nr_variable(); - feature_node **x=prob->x; - - for(i=0;iindex!=-1) - { - XTv[s->index-1]+=v[i]*s->value; - s++; - } - } -} - -class l2r_l2_svc_fun: public function -{ -public: - l2r_l2_svc_fun(const problem *prob, double *C); - ~l2r_l2_svc_fun(); - - double fun(double *w); - void grad(double *w, double *g); - void Hv(double *s, double *Hs); - - int get_nr_variable(void); - -protected: - void Xv(double *v, double *Xv); - void subXv(double *v, double *Xv); - void subXTv(double *v, double *XTv); - - double *C; - double *z; - double *D; - int *I; - int sizeI; - const problem *prob; -}; - -l2r_l2_svc_fun::l2r_l2_svc_fun(const problem *prob, double *C) -{ - int l=prob->l; - - this->prob = prob; - - z = new double[l]; - D = new double[l]; - I = new int[l]; - this->C = C; -} - -l2r_l2_svc_fun::~l2r_l2_svc_fun() -{ - delete[] z; - delete[] D; - delete[] I; -} - -double l2r_l2_svc_fun::fun(double *w) -{ - int i; - double f=0; - double *y=prob->y; - int l=prob->l; - int w_size=get_nr_variable(); - - Xv(w, z); - - for(i=0;i 0) - f += C[i]*d*d; - } - - return(f); -} - -void l2r_l2_svc_fun::grad(double *w, double *g) -{ - int i; - double *y=prob->y; - int l=prob->l; - int w_size=get_nr_variable(); - - sizeI = 0; - for (i=0;in; -} - -void l2r_l2_svc_fun::Hv(double *s, double *Hs) -{ - int i; - int w_size=get_nr_variable(); - double *wa = new double[sizeI]; - - subXv(s, wa); - for(i=0;il; - feature_node **x=prob->x; - - for(i=0;iindex!=-1) - { - Xv[i]+=v[s->index-1]*s->value; - s++; - } - } -} - -void l2r_l2_svc_fun::subXv(double *v, double *Xv) -{ - int i; - feature_node **x=prob->x; - - for(i=0;iindex!=-1) - { - Xv[i]+=v[s->index-1]*s->value; - s++; - } - } -} - -void l2r_l2_svc_fun::subXTv(double *v, double *XTv) -{ - int i; - int w_size=get_nr_variable(); - feature_node **x=prob->x; - - for(i=0;iindex!=-1) - { - XTv[s->index-1]+=v[i]*s->value; - s++; - } - } -} - -class l2r_l2_svr_fun: public l2r_l2_svc_fun -{ -public: - l2r_l2_svr_fun(const problem *prob, double *C, double p); - - double fun(double *w); - void grad(double *w, double *g); - -private: - double p; -}; - -l2r_l2_svr_fun::l2r_l2_svr_fun(const problem *prob, double *C, double p): - l2r_l2_svc_fun(prob, C) -{ - this->p = p; -} - -double l2r_l2_svr_fun::fun(double *w) -{ - int i; - double f=0; - double *y=prob->y; - int l=prob->l; - int w_size=get_nr_variable(); - double d; - - Xv(w, z); - - for(i=0;i p) - f += C[i]*(d-p)*(d-p); - } - - return(f); -} - -void l2r_l2_svr_fun::grad(double *w, double *g) -{ - int i; - double *y=prob->y; - int l=prob->l; - int w_size=get_nr_variable(); - double d; - - sizeI = 0; - for(i=0;i p) - { - z[sizeI] = C[i]*(d-p); - I[sizeI] = i; - sizeI++; - } - - } - subXTv(z, g); - - for(i=0;iy[i]) -// To support weights for instances, use GETI(i) (i) - -class Solver_MCSVM_CS -{ - public: - Solver_MCSVM_CS(const problem *prob, int nr_class, double *C, double eps=0.1, int max_iter=100000); - ~Solver_MCSVM_CS(); - void Solve(double *w); - private: - void solve_sub_problem(double A_i, int yi, double C_yi, int active_i, double *alpha_new); - bool be_shrunk(int i, int m, int yi, double alpha_i, double minG); - double *B, *C, *G; - int w_size, l; - int nr_class; - int max_iter; - double eps; - const problem *prob; -}; - -Solver_MCSVM_CS::Solver_MCSVM_CS(const problem *prob, int nr_class, double *weighted_C, double eps, int max_iter) -{ - this->w_size = prob->n; - this->l = prob->l; - this->nr_class = nr_class; - this->eps = eps; - this->max_iter = max_iter; - this->prob = prob; - this->B = new double[nr_class]; - this->G = new double[nr_class]; - this->C = weighted_C; -} - -Solver_MCSVM_CS::~Solver_MCSVM_CS() -{ - delete[] B; - delete[] G; -} - -int compare_double(const void *a, const void *b) -{ - if(*(double *)a > *(double *)b) - return -1; - if(*(double *)a < *(double *)b) - return 1; - return 0; -} - -void Solver_MCSVM_CS::solve_sub_problem(double A_i, int yi, double C_yi, int active_i, double *alpha_new) -{ - int r; - double *D; - - clone(D, B, active_i); - if(yi < active_i) - D[yi] += A_i*C_yi; - qsort(D, active_i, sizeof(double), compare_double); - - double beta = D[0] - A_i*C_yi; - for(r=1;ry[i] == m - // alpha[i*nr_class+m] <= 0 if prob->y[i] != m - // If initial alpha isn't zero, uncomment the for loop below to initialize w - for(i=0;ix[i]; - QD[i] = 0; - while(xi->index != -1) - { - double val = xi->value; - QD[i] += val*val; - - // Uncomment the for loop if initial alpha isn't zero - // for(m=0; mindex-1)*nr_class+m] += alpha[i*nr_class+m]*val; - xi++; - } - active_size_i[i] = nr_class; - y_index[i] = (int)prob->y[i]; - index[i] = i; - } - - while(iter < max_iter) - { - double stopping = -INF; - for(i=0;i 0) - { - for(m=0;mx[i]; - while(xi->index!= -1) - { - double *w_i = &w[(xi->index-1)*nr_class]; - for(m=0;mvalue); - xi++; - } - - double minG = INF; - double maxG = -INF; - for(m=0;m maxG) - maxG = G[m]; - } - if(y_index[i] < active_size_i[i]) - if(alpha_i[(int) prob->y[i]] < C[GETI(i)] && G[y_index[i]] < minG) - minG = G[y_index[i]]; - - for(m=0;mm) - { - if(!be_shrunk(i, active_size_i[i], y_index[i], - alpha_i[alpha_index_i[active_size_i[i]]], minG)) - { - swap(alpha_index_i[m], alpha_index_i[active_size_i[i]]); - swap(G[m], G[active_size_i[i]]); - if(y_index[i] == active_size_i[i]) - y_index[i] = m; - else if(y_index[i] == m) - y_index[i] = active_size_i[i]; - break; - } - active_size_i[i]--; - } - } - } - - if(active_size_i[i] <= 1) - { - active_size--; - swap(index[s], index[active_size]); - s--; - continue; - } - - if(maxG-minG <= 1e-12) - continue; - else - stopping = max(maxG - minG, stopping); - - for(m=0;m= 1e-12) - { - d_ind[nz_d] = alpha_index_i[m]; - d_val[nz_d] = d; - nz_d++; - } - } - - xi = prob->x[i]; - while(xi->index != -1) - { - double *w_i = &w[(xi->index-1)*nr_class]; - for(m=0;mvalue; - xi++; - } - } - } - - iter++; - if(iter % 10 == 0) - { - info("."); - } - - if(stopping < eps_shrink) - { - if(stopping < eps && start_from_all == true) - break; - else - { - active_size = l; - for(i=0;i= max_iter) - info("\nWARNING: reaching max number of iterations\n"); - - // calculate objective value - double v = 0; - int nSV = 0; - for(i=0;i 0) - nSV++; - } - for(i=0;iy[i]]; - info("Objective value = %lf\n",v); - info("nSV = %d\n",nSV); - - delete [] alpha; - delete [] alpha_new; - delete [] index; - delete [] QD; - delete [] d_ind; - delete [] d_val; - delete [] alpha_index; - delete [] y_index; - delete [] active_size_i; -} - -// A coordinate descent algorithm for -// L1-loss and L2-loss SVM dual problems -// -// min_\alpha 0.5(\alpha^T (Q + D)\alpha) - e^T \alpha, -// s.t. 0 <= \alpha_i <= upper_bound_i, -// -// where Qij = yi yj xi^T xj and -// D is a diagonal matrix -// -// In L1-SVM case: -// upper_bound_i = Cp if y_i = 1 -// upper_bound_i = Cn if y_i = -1 -// D_ii = 0 -// In L2-SVM case: -// upper_bound_i = INF -// D_ii = 1/(2*Cp) if y_i = 1 -// D_ii = 1/(2*Cn) if y_i = -1 -// -// Given: -// x, y, Cp, Cn -// eps is the stopping tolerance -// -// solution will be put in w -// -// See Algorithm 3 of Hsieh et al., ICML 2008 - -#undef GETI -#define GETI(i) (y[i]+1) -// To support weights for instances, use GETI(i) (i) - -static void solve_l2r_l1l2_svc( - const problem *prob, double *w, double eps, - double Cp, double Cn, int solver_type) -{ - int l = prob->l; - int w_size = prob->n; - int i, s, iter = 0; - double C, d, G; - double *QD = new double[l]; - int max_iter = 1000; - int *index = new int[l]; - double *alpha = new double[l]; - schar *y = new schar[l]; - int active_size = l; - - // PG: projected gradient, for shrinking and stopping - double PG; - double PGmax_old = INF; - double PGmin_old = -INF; - double PGmax_new, PGmin_new; - - // default solver_type: L2R_L2LOSS_SVC_DUAL - double diag[3] = {0.5/Cn, 0, 0.5/Cp}; - double upper_bound[3] = {INF, 0, INF}; - if(solver_type == L2R_L1LOSS_SVC_DUAL) - { - diag[0] = 0; - diag[2] = 0; - upper_bound[0] = Cn; - upper_bound[2] = Cp; - } - - for(i=0; iy[i] > 0) - { - y[i] = +1; - } - else - { - y[i] = -1; - } - } - - // Initial alpha can be set here. Note that - // 0 <= alpha[i] <= upper_bound[GETI(i)] - for(i=0; ix[i]; - while (xi->index != -1) - { - double val = xi->value; - QD[i] += val*val; - w[xi->index-1] += y[i]*alpha[i]*val; - xi++; - } - index[i] = i; - } - - while (iter < max_iter) - { - PGmax_new = -INF; - PGmin_new = INF; - - for (i=0; ix[i]; - while(xi->index!= -1) - { - G += w[xi->index-1]*(xi->value); - xi++; - } - G = G*yi-1; - - C = upper_bound[GETI(i)]; - G += alpha[i]*diag[GETI(i)]; - - PG = 0; - if (alpha[i] == 0) - { - if (G > PGmax_old) - { - active_size--; - swap(index[s], index[active_size]); - s--; - continue; - } - else if (G < 0) - PG = G; - } - else if (alpha[i] == C) - { - if (G < PGmin_old) - { - active_size--; - swap(index[s], index[active_size]); - s--; - continue; - } - else if (G > 0) - PG = G; - } - else - PG = G; - - PGmax_new = max(PGmax_new, PG); - PGmin_new = min(PGmin_new, PG); - - if(fabs(PG) > 1.0e-12) - { - double alpha_old = alpha[i]; - alpha[i] = min(max(alpha[i] - G/QD[i], 0.0), C); - d = (alpha[i] - alpha_old)*yi; - xi = prob->x[i]; - while (xi->index != -1) - { - w[xi->index-1] += d*xi->value; - xi++; - } - } - } - - iter++; - if(iter % 10 == 0) - info("."); - - if(PGmax_new - PGmin_new <= eps) - { - if(active_size == l) - break; - else - { - active_size = l; - info("*"); - PGmax_old = INF; - PGmin_old = -INF; - continue; - } - } - PGmax_old = PGmax_new; - PGmin_old = PGmin_new; - if (PGmax_old <= 0) - PGmax_old = INF; - if (PGmin_old >= 0) - PGmin_old = -INF; - } - - info("\noptimization finished, #iter = %d\n",iter); - if (iter >= max_iter) - info("\nWARNING: reaching max number of iterations\nUsing -s 2 may be faster (also see FAQ)\n\n"); - - // calculate objective value - - double v = 0; - int nSV = 0; - for(i=0; i 0) - ++nSV; - } - info("Objective value = %lf\n",v/2); - info("nSV = %d\n",nSV); - - delete [] QD; - delete [] alpha; - delete [] y; - delete [] index; -} - - -// A coordinate descent algorithm for -// L1-loss and L2-loss epsilon-SVR dual problem -// -// min_\beta 0.5\beta^T (Q + diag(lambda)) \beta - p \sum_{i=1}^l|\beta_i| + \sum_{i=1}^l yi\beta_i, -// s.t. -upper_bound_i <= \beta_i <= upper_bound_i, -// -// where Qij = xi^T xj and -// D is a diagonal matrix -// -// In L1-SVM case: -// upper_bound_i = C -// lambda_i = 0 -// In L2-SVM case: -// upper_bound_i = INF -// lambda_i = 1/(2*C) -// -// Given: -// x, y, p, C -// eps is the stopping tolerance -// -// solution will be put in w -// -// See Algorithm 4 of Ho and Lin, 2012 - -#undef GETI -#define GETI(i) (0) -// To support weights for instances, use GETI(i) (i) - -static void solve_l2r_l1l2_svr( - const problem *prob, double *w, const parameter *param, - int solver_type) -{ - int l = prob->l; - double C = param->C; - double p = param->p; - int w_size = prob->n; - double eps = param->eps; - int i, s, iter = 0; - int max_iter = 1000; - int active_size = l; - int *index = new int[l]; - - double d, G, H; - double Gmax_old = INF; - double Gmax_new, Gnorm1_new; - double Gnorm1_init; - double *beta = new double[l]; - double *QD = new double[l]; - double *y = prob->y; - - // L2R_L2LOSS_SVR_DUAL - double lambda[1], upper_bound[1]; - lambda[0] = 0.5/C; - upper_bound[0] = INF; - - if(solver_type == L2R_L1LOSS_SVR_DUAL) - { - lambda[0] = 0; - upper_bound[0] = C; - } - - // Initial beta can be set here. Note that - // -upper_bound <= beta[i] <= upper_bound - for(i=0; ix[i]; - while(xi->index != -1) - { - double val = xi->value; - QD[i] += val*val; - w[xi->index-1] += beta[i]*val; - xi++; - } - - index[i] = i; - } - - - while(iter < max_iter) - { - Gmax_new = 0; - Gnorm1_new = 0; - - for(i=0; ix[i]; - while(xi->index != -1) - { - int ind = xi->index-1; - double val = xi->value; - G += val*w[ind]; - xi++; - } - - double Gp = G+p; - double Gn = G-p; - double violation = 0; - if(beta[i] == 0) - { - if(Gp < 0) - violation = -Gp; - else if(Gn > 0) - violation = Gn; - else if(Gp>Gmax_old && Gn<-Gmax_old) - { - active_size--; - swap(index[s], index[active_size]); - s--; - continue; - } - } - else if(beta[i] >= upper_bound[GETI(i)]) - { - if(Gp > 0) - violation = Gp; - else if(Gp < -Gmax_old) - { - active_size--; - swap(index[s], index[active_size]); - s--; - continue; - } - } - else if(beta[i] <= -upper_bound[GETI(i)]) - { - if(Gn < 0) - violation = -Gn; - else if(Gn > Gmax_old) - { - active_size--; - swap(index[s], index[active_size]); - s--; - continue; - } - } - else if(beta[i] > 0) - violation = fabs(Gp); - else - violation = fabs(Gn); - - Gmax_new = max(Gmax_new, violation); - Gnorm1_new += violation; - - // obtain Newton direction d - if(Gp < H*beta[i]) - d = -Gp/H; - else if(Gn > H*beta[i]) - d = -Gn/H; - else - d = -beta[i]; - - if(fabs(d) < 1.0e-12) - continue; - - double beta_old = beta[i]; - beta[i] = min(max(beta[i]+d, -upper_bound[GETI(i)]), upper_bound[GETI(i)]); - d = beta[i]-beta_old; - - if(d != 0) - { - xi = prob->x[i]; - while(xi->index != -1) - { - w[xi->index-1] += d*xi->value; - xi++; - } - } - } - - if(iter == 0) - Gnorm1_init = Gnorm1_new; - iter++; - if(iter % 10 == 0) - info("."); - - if(Gnorm1_new <= eps*Gnorm1_init) - { - if(active_size == l) - break; - else - { - active_size = l; - info("*"); - Gmax_old = INF; - continue; - } - } - - Gmax_old = Gmax_new; - } - - info("\noptimization finished, #iter = %d\n", iter); - if(iter >= max_iter) - info("\nWARNING: reaching max number of iterations\nUsing -s 11 may be faster\n\n"); - - // calculate objective value - double v = 0; - int nSV = 0; - for(i=0; il; - int w_size = prob->n; - int i, s, iter = 0; - double *xTx = new double[l]; - int max_iter = 1000; - int *index = new int[l]; - double *alpha = new double[2*l]; // store alpha and C - alpha - schar *y = new schar[l]; - int max_inner_iter = 100; // for inner Newton - double innereps = 1e-2; - double innereps_min = min(1e-8, eps); - double upper_bound[3] = {Cn, 0, Cp}; - - for(i=0; iy[i] > 0) - { - y[i] = +1; - } - else - { - y[i] = -1; - } - } - - // Initial alpha can be set here. Note that - // 0 < alpha[i] < upper_bound[GETI(i)] - // alpha[2*i] + alpha[2*i+1] = upper_bound[GETI(i)] - for(i=0; ix[i]; - while (xi->index != -1) - { - double val = xi->value; - xTx[i] += val*val; - w[xi->index-1] += y[i]*alpha[2*i]*val; - xi++; - } - index[i] = i; - } - - while (iter < max_iter) - { - for (i=0; ix[i]; - while (xi->index != -1) - { - ywTx += w[xi->index-1]*xi->value; - xi++; - } - ywTx *= y[i]; - double a = xisq, b = ywTx; - - // Decide to minimize g_1(z) or g_2(z) - int ind1 = 2*i, ind2 = 2*i+1, sign = 1; - if(0.5*a*(alpha[ind2]-alpha[ind1])+b < 0) - { - ind1 = 2*i+1; - ind2 = 2*i; - sign = -1; - } - - // g_t(z) = z*log(z) + (C-z)*log(C-z) + 0.5a(z-alpha_old)^2 + sign*b(z-alpha_old) - double alpha_old = alpha[ind1]; - double z = alpha_old; - if(C - z < 0.5 * C) - z = 0.1*z; - double gp = a*(z-alpha_old)+sign*b+log(z/(C-z)); - Gmax = max(Gmax, fabs(gp)); - - // Newton method on the sub-problem - const double eta = 0.1; // xi in the paper - int inner_iter = 0; - while (inner_iter <= max_inner_iter) - { - if(fabs(gp) < innereps) - break; - double gpp = a + C/(C-z)/z; - double tmpz = z - gp/gpp; - if(tmpz <= 0) - z *= eta; - else // tmpz in (0, C) - z = tmpz; - gp = a*(z-alpha_old)+sign*b+log(z/(C-z)); - newton_iter++; - inner_iter++; - } - - if(inner_iter > 0) // update w - { - alpha[ind1] = z; - alpha[ind2] = C-z; - xi = prob->x[i]; - while (xi->index != -1) - { - w[xi->index-1] += sign*(z-alpha_old)*yi*xi->value; - xi++; - } - } - } - - iter++; - if(iter % 10 == 0) - info("."); - - if(Gmax < eps) - break; - - if(newton_iter <= l/10) - innereps = max(innereps_min, 0.1*innereps); - - } - - info("\noptimization finished, #iter = %d\n",iter); - if (iter >= max_iter) - info("\nWARNING: reaching max number of iterations\nUsing -s 0 may be faster (also see FAQ)\n\n"); - - // calculate objective value - - double v = 0; - for(i=0; il; - int w_size = prob_col->n; - int j, s, iter = 0; - int max_iter = 1000; - int active_size = w_size; - int max_num_linesearch = 20; - - double sigma = 0.01; - double d, G_loss, G, H; - double Gmax_old = INF; - double Gmax_new, Gnorm1_new; - double Gnorm1_init; - double d_old, d_diff; - double loss_old, loss_new; - double appxcond, cond; - - int *index = new int[w_size]; - schar *y = new schar[l]; - double *b = new double[l]; // b = 1-ywTx - double *xj_sq = new double[w_size]; - feature_node *x; - - double C[3] = {Cn,0,Cp}; - - // Initial w can be set here. - for(j=0; jy[j] > 0) - y[j] = 1; - else - y[j] = -1; - } - for(j=0; jx[j]; - while(x->index != -1) - { - int ind = x->index-1; - x->value *= y[ind]; // x->value stores yi*xij - double val = x->value; - b[ind] -= w[j]*val; - xj_sq[j] += C[GETI(ind)]*val*val; - x++; - } - } - - while(iter < max_iter) - { - Gmax_new = 0; - Gnorm1_new = 0; - - for(j=0; jx[j]; - while(x->index != -1) - { - int ind = x->index-1; - if(b[ind] > 0) - { - double val = x->value; - double tmp = C[GETI(ind)]*val; - G_loss -= tmp*b[ind]; - H += tmp*val; - } - x++; - } - G_loss *= 2; - - G = G_loss; - H *= 2; - H = max(H, 1e-12); - - double Gp = G+1; - double Gn = G-1; - double violation = 0; - if(w[j] == 0) - { - if(Gp < 0) - violation = -Gp; - else if(Gn > 0) - violation = Gn; - else if(Gp>Gmax_old/l && Gn<-Gmax_old/l) - { - active_size--; - swap(index[s], index[active_size]); - s--; - continue; - } - } - else if(w[j] > 0) - violation = fabs(Gp); - else - violation = fabs(Gn); - - Gmax_new = max(Gmax_new, violation); - Gnorm1_new += violation; - - // obtain Newton direction d - if(Gp < H*w[j]) - d = -Gp/H; - else if(Gn > H*w[j]) - d = -Gn/H; - else - d = -w[j]; - - if(fabs(d) < 1.0e-12) - continue; - - double delta = fabs(w[j]+d)-fabs(w[j]) + G*d; - d_old = 0; - int num_linesearch; - for(num_linesearch=0; num_linesearch < max_num_linesearch; num_linesearch++) - { - d_diff = d_old - d; - cond = fabs(w[j]+d)-fabs(w[j]) - sigma*delta; - - appxcond = xj_sq[j]*d*d + G_loss*d + cond; - if(appxcond <= 0) - { - x = prob_col->x[j]; - while(x->index != -1) - { - b[x->index-1] += d_diff*x->value; - x++; - } - break; - } - - if(num_linesearch == 0) - { - loss_old = 0; - loss_new = 0; - x = prob_col->x[j]; - while(x->index != -1) - { - int ind = x->index-1; - if(b[ind] > 0) - loss_old += C[GETI(ind)]*b[ind]*b[ind]; - double b_new = b[ind] + d_diff*x->value; - b[ind] = b_new; - if(b_new > 0) - loss_new += C[GETI(ind)]*b_new*b_new; - x++; - } - } - else - { - loss_new = 0; - x = prob_col->x[j]; - while(x->index != -1) - { - int ind = x->index-1; - double b_new = b[ind] + d_diff*x->value; - b[ind] = b_new; - if(b_new > 0) - loss_new += C[GETI(ind)]*b_new*b_new; - x++; - } - } - - cond = cond + loss_new - loss_old; - if(cond <= 0) - break; - else - { - d_old = d; - d *= 0.5; - delta *= 0.5; - } - } - - w[j] += d; - - // recompute b[] if line search takes too many steps - if(num_linesearch >= max_num_linesearch) - { - info("#"); - for(int i=0; ix[i]; - while(x->index != -1) - { - b[x->index-1] -= w[i]*x->value; - x++; - } - } - } - } - - if(iter == 0) - Gnorm1_init = Gnorm1_new; - iter++; - if(iter % 10 == 0) - info("."); - - if(Gnorm1_new <= eps*Gnorm1_init) - { - if(active_size == w_size) - break; - else - { - active_size = w_size; - info("*"); - Gmax_old = INF; - continue; - } - } - - Gmax_old = Gmax_new; - } - - info("\noptimization finished, #iter = %d\n", iter); - if(iter >= max_iter) - info("\nWARNING: reaching max number of iterations\n"); - - // calculate objective value - - double v = 0; - int nnz = 0; - for(j=0; jx[j]; - while(x->index != -1) - { - x->value *= prob_col->y[x->index-1]; // restore x->value - x++; - } - if(w[j] != 0) - { - v += fabs(w[j]); - nnz++; - } - } - for(j=0; j 0) - v += C[GETI(j)]*b[j]*b[j]; - - info("Objective value = %lf\n", v); - info("#nonzeros/#features = %d/%d\n", nnz, w_size); - - delete [] index; - delete [] y; - delete [] b; - delete [] xj_sq; -} - -// A coordinate descent algorithm for -// L1-regularized logistic regression problems -// -// min_w \sum |wj| + C \sum log(1+exp(-yi w^T xi)), -// -// Given: -// x, y, Cp, Cn -// eps is the stopping tolerance -// -// solution will be put in w -// -// See Yuan et al. (2011) and appendix of LIBLINEAR paper, Fan et al. (2008) - -#undef GETI -#define GETI(i) (y[i]+1) -// To support weights for instances, use GETI(i) (i) - -static void solve_l1r_lr( - const problem *prob_col, double *w, double eps, - double Cp, double Cn) -{ - int l = prob_col->l; - int w_size = prob_col->n; - int j, s, newton_iter=0, iter=0; - int max_newton_iter = 100; - int max_iter = 1000; - int max_num_linesearch = 20; - int active_size; - int QP_active_size; - - double nu = 1e-12; - double inner_eps = 1; - double sigma = 0.01; - double w_norm, w_norm_new; - double z, G, H; - double Gnorm1_init; - double Gmax_old = INF; - double Gmax_new, Gnorm1_new; - double QP_Gmax_old = INF; - double QP_Gmax_new, QP_Gnorm1_new; - double delta, negsum_xTd, cond; - - int *index = new int[w_size]; - schar *y = new schar[l]; - double *Hdiag = new double[w_size]; - double *Grad = new double[w_size]; - double *wpd = new double[w_size]; - double *xjneg_sum = new double[w_size]; - double *xTd = new double[l]; - double *exp_wTx = new double[l]; - double *exp_wTx_new = new double[l]; - double *tau = new double[l]; - double *D = new double[l]; - feature_node *x; - - double C[3] = {Cn,0,Cp}; - - // Initial w can be set here. - for(j=0; jy[j] > 0) - y[j] = 1; - else - y[j] = -1; - - exp_wTx[j] = 0; - } - - w_norm = 0; - for(j=0; jx[j]; - while(x->index != -1) - { - int ind = x->index-1; - double val = x->value; - exp_wTx[ind] += w[j]*val; - if(y[ind] == -1) - xjneg_sum[j] += C[GETI(ind)]*val; - x++; - } - } - for(j=0; jx[j]; - while(x->index != -1) - { - int ind = x->index-1; - Hdiag[j] += x->value*x->value*D[ind]; - tmp += x->value*tau[ind]; - x++; - } - Grad[j] = -tmp + xjneg_sum[j]; - - double Gp = Grad[j]+1; - double Gn = Grad[j]-1; - double violation = 0; - if(w[j] == 0) - { - if(Gp < 0) - violation = -Gp; - else if(Gn > 0) - violation = Gn; - //outer-level shrinking - else if(Gp>Gmax_old/l && Gn<-Gmax_old/l) - { - active_size--; - swap(index[s], index[active_size]); - s--; - continue; - } - } - else if(w[j] > 0) - violation = fabs(Gp); - else - violation = fabs(Gn); - - Gmax_new = max(Gmax_new, violation); - Gnorm1_new += violation; - } - - if(newton_iter == 0) - Gnorm1_init = Gnorm1_new; - - if(Gnorm1_new <= eps*Gnorm1_init) - break; - - iter = 0; - QP_Gmax_old = INF; - QP_active_size = active_size; - - for(int i=0; ix[j]; - G = Grad[j] + (wpd[j]-w[j])*nu; - while(x->index != -1) - { - int ind = x->index-1; - G += x->value*D[ind]*xTd[ind]; - x++; - } - - double Gp = G+1; - double Gn = G-1; - double violation = 0; - if(wpd[j] == 0) - { - if(Gp < 0) - violation = -Gp; - else if(Gn > 0) - violation = Gn; - //inner-level shrinking - else if(Gp>QP_Gmax_old/l && Gn<-QP_Gmax_old/l) - { - QP_active_size--; - swap(index[s], index[QP_active_size]); - s--; - continue; - } - } - else if(wpd[j] > 0) - violation = fabs(Gp); - else - violation = fabs(Gn); - - QP_Gmax_new = max(QP_Gmax_new, violation); - QP_Gnorm1_new += violation; - - // obtain solution of one-variable problem - if(Gp < H*wpd[j]) - z = -Gp/H; - else if(Gn > H*wpd[j]) - z = -Gn/H; - else - z = -wpd[j]; - - if(fabs(z) < 1.0e-12) - continue; - z = min(max(z,-10.0),10.0); - - wpd[j] += z; - - x = prob_col->x[j]; - while(x->index != -1) - { - int ind = x->index-1; - xTd[ind] += x->value*z; - x++; - } - } - - iter++; - - if(QP_Gnorm1_new <= inner_eps*Gnorm1_init) - { - //inner stopping - if(QP_active_size == active_size) - break; - //active set reactivation - else - { - QP_active_size = active_size; - QP_Gmax_old = INF; - continue; - } - } - - QP_Gmax_old = QP_Gmax_new; - } - - if(iter >= max_iter) - info("WARNING: reaching max number of inner iterations\n"); - - delta = 0; - w_norm_new = 0; - for(j=0; j= max_num_linesearch) - { - for(int i=0; ix[i]; - while(x->index != -1) - { - exp_wTx[x->index-1] += w[i]*x->value; - x++; - } - } - - for(int i=0; i= max_newton_iter) - info("WARNING: reaching max number of iterations\n"); - - // calculate objective value - - double v = 0; - int nnz = 0; - for(j=0; jl; - int n = prob->n; - long int nnz = 0; - long int *col_ptr = new long int [n+1]; - feature_node *x_space; - prob_col->l = l; - prob_col->n = n; - prob_col->y = new double[l]; - prob_col->x = new feature_node*[n]; - - for(i=0; iy[i] = prob->y[i]; - - for(i=0; ix[i]; - while(x->index != -1) - { - nnz++; - col_ptr[x->index]++; - x++; - } - } - for(i=1; ix[i] = &x_space[col_ptr[i]]; - - for(i=0; ix[i]; - while(x->index != -1) - { - int ind = x->index-1; - x_space[col_ptr[ind]].index = i+1; // starts from 1 - x_space[col_ptr[ind]].value = x->value; - col_ptr[ind]++; - x++; - } - } - for(i=0; il; - int max_nr_class = 16; - int nr_class = 0; - int *label = Malloc(int,max_nr_class); - int *count = Malloc(int,max_nr_class); - int *data_label = Malloc(int,l); - int i; - - for(i=0;iy[i]; - int j; - for(j=0;jeps; - int pos = 0; - int neg = 0; - for(int i=0;il;i++) - if(prob->y[i] > 0) - pos++; - neg = prob->l - pos; - - double primal_solver_tol = eps*max(min(pos,neg), 1)/prob->l; - - function *fun_obj=NULL; - switch(param->solver_type) - { - case L2R_LR: - { - double *C = new double[prob->l]; - for(int i = 0; i < prob->l; i++) - { - if(prob->y[i] > 0) - C[i] = Cp; - else - C[i] = Cn; - } - fun_obj=new l2r_lr_fun(prob, C); - TRON tron_obj(fun_obj, primal_solver_tol); - tron_obj.set_print_string(liblinear_print_string); - tron_obj.tron(w); - delete fun_obj; - delete[] C; - break; - } - case L2R_L2LOSS_SVC: - { - double *C = new double[prob->l]; - for(int i = 0; i < prob->l; i++) - { - if(prob->y[i] > 0) - C[i] = Cp; - else - C[i] = Cn; - } - fun_obj=new l2r_l2_svc_fun(prob, C); - TRON tron_obj(fun_obj, primal_solver_tol); - tron_obj.set_print_string(liblinear_print_string); - tron_obj.tron(w); - delete fun_obj; - delete[] C; - break; - } - case L2R_L2LOSS_SVC_DUAL: - solve_l2r_l1l2_svc(prob, w, eps, Cp, Cn, L2R_L2LOSS_SVC_DUAL); - break; - case L2R_L1LOSS_SVC_DUAL: - solve_l2r_l1l2_svc(prob, w, eps, Cp, Cn, L2R_L1LOSS_SVC_DUAL); - break; - case L1R_L2LOSS_SVC: - { - problem prob_col; - feature_node *x_space = NULL; - transpose(prob, &x_space ,&prob_col); - solve_l1r_l2_svc(&prob_col, w, primal_solver_tol, Cp, Cn); - delete [] prob_col.y; - delete [] prob_col.x; - delete [] x_space; - break; - } - case L1R_LR: - { - problem prob_col; - feature_node *x_space = NULL; - transpose(prob, &x_space ,&prob_col); - solve_l1r_lr(&prob_col, w, primal_solver_tol, Cp, Cn); - delete [] prob_col.y; - delete [] prob_col.x; - delete [] x_space; - break; - } - case L2R_LR_DUAL: - solve_l2r_lr_dual(prob, w, eps, Cp, Cn); - break; - case L2R_L2LOSS_SVR: - { - double *C = new double[prob->l]; - for(int i = 0; i < prob->l; i++) - C[i] = param->C; - - fun_obj=new l2r_l2_svr_fun(prob, C, param->p); - TRON tron_obj(fun_obj, param->eps); - tron_obj.set_print_string(liblinear_print_string); - tron_obj.tron(w); - delete fun_obj; - delete[] C; - break; - - } - case L2R_L1LOSS_SVR_DUAL: - solve_l2r_l1l2_svr(prob, w, param, L2R_L1LOSS_SVR_DUAL); - break; - case L2R_L2LOSS_SVR_DUAL: - solve_l2r_l1l2_svr(prob, w, param, L2R_L2LOSS_SVR_DUAL); - break; - default: - fprintf(stderr, "ERROR: unknown solver_type\n"); - break; - } -} - -// -// Interface functions -// -model* train(const problem *prob, const parameter *param) -{ - int i,j; - int l = prob->l; - int n = prob->n; - int w_size = prob->n; - - model *model_ = Malloc(model,1); - - if(prob->bias>=0) - model_->nr_feature=n-1; - else - model_->nr_feature=n; - model_->param = *param; - model_->bias = prob->bias; - - if(param->solver_type == L2R_L2LOSS_SVR || - param->solver_type == L2R_L1LOSS_SVR_DUAL || - param->solver_type == L2R_L2LOSS_SVR_DUAL) - { - model_->w = Malloc(double, w_size); - model_->nr_class = 2; - model_->label = NULL; - train_one(prob, param, &model_->w[0], 0, 0); - } - else - { - int nr_class; - int *label = NULL; - int *start = NULL; - int *count = NULL; - int *perm = Malloc(int,l); - - // group training data of the same class - group_classes(prob,&nr_class,&label,&start,&count,perm); - - model_->nr_class=nr_class; - model_->label = Malloc(int,nr_class); - for(i=0;ilabel[i] = label[i]; - - // calculate weighted C - double *weighted_C = Malloc(double, nr_class); - for(i=0;iC; - for(i=0;inr_weight;i++) - { - for(j=0;jweight_label[i] == label[j]) - break; - if(j == nr_class) - fprintf(stderr,"WARNING: class label %d specified in weight is not found\n", param->weight_label[i]); - else - weighted_C[j] *= param->weight[i]; - } - - // constructing the subproblem - feature_node **x = Malloc(feature_node *,l); - - for(i=0;ix[perm[i]]; - } - - int k; - problem sub_prob; - sub_prob.l = l; - sub_prob.n = n; - sub_prob.x = Malloc(feature_node *,sub_prob.l); - sub_prob.y = Malloc(double,sub_prob.l); - - for(k=0; ksolver_type == MCSVM_CS) - { - model_->w=Malloc(double, n*nr_class); - for(i=0;ieps); - Solver.Solve(model_->w); - } - else - { - if(nr_class == 2) - { - model_->w=Malloc(double, w_size); - - int e0 = start[0]+count[0]; - k=0; - for(; kw[0], weighted_C[0], weighted_C[1]); - } - else - { - model_->w=Malloc(double, w_size*nr_class); - double *w=Malloc(double, w_size); - for(i=0;iC); - - for(int j=0;jw[j*nr_class+i] = w[j]; - } - free(w); - } - - } - - free(x); - free(label); - free(start); - free(count); - free(perm); - free(sub_prob.x); - free(sub_prob.y); - free(weighted_C); - } - return model_; -} - -void cross_validation(const problem *prob, const parameter *param, int nr_fold, double *target) -{ - int i; - int *fold_start; - int l = prob->l; - int *perm = Malloc(int,l); - if (nr_fold > l) - { - nr_fold = l; - fprintf(stderr,"WARNING: # folds > # data. Will use # folds = # data instead (i.e., leave-one-out cross validation)\n"); - } - fold_start = Malloc(int,nr_fold+1); - for(i=0;ibias; - subprob.n = prob->n; - subprob.l = l-(end-begin); - subprob.x = Malloc(struct feature_node*,subprob.l); - subprob.y = Malloc(double,subprob.l); - - k=0; - for(j=0;jx[perm[j]]; - subprob.y[k] = prob->y[perm[j]]; - ++k; - } - for(j=end;jx[perm[j]]; - subprob.y[k] = prob->y[perm[j]]; - ++k; - } - struct model *submodel = train(&subprob,param); - for(j=begin;jx[perm[j]]); - free_and_destroy_model(&submodel); - free(subprob.x); - free(subprob.y); - } - free(fold_start); - free(perm); -} - -double predict_values(const struct model *model_, const struct feature_node *x, double *dec_values) -{ - int idx; - int n; - if(model_->bias>=0) - n=model_->nr_feature+1; - else - n=model_->nr_feature; - double *w=model_->w; - int nr_class=model_->nr_class; - int i; - int nr_w; - if(nr_class==2 && model_->param.solver_type != MCSVM_CS) - nr_w = 1; - else - nr_w = nr_class; - - const feature_node *lx=x; - for(i=0;iindex)!=-1; lx++) - { - // the dimension of testing data may exceed that of training - if(idx<=n) - for(i=0;ivalue; - } - - if(nr_class==2) - { - if(model_->param.solver_type == L2R_L2LOSS_SVR || - model_->param.solver_type == L2R_L1LOSS_SVR_DUAL || - model_->param.solver_type == L2R_L2LOSS_SVR_DUAL) - return dec_values[0]; - else - return (dec_values[0]>0)?model_->label[0]:model_->label[1]; - } - else - { - int dec_max_idx = 0; - for(i=1;i dec_values[dec_max_idx]) - dec_max_idx = i; - } - return model_->label[dec_max_idx]; - } -} - -double predict(const model *model_, const feature_node *x) -{ - double *dec_values = Malloc(double, model_->nr_class); - double label=predict_values(model_, x, dec_values); - free(dec_values); - return label; -} - -double predict_probability(const struct model *model_, const struct feature_node *x, double* prob_estimates) -{ - if(check_probability_model(model_)) - { - int i; - int nr_class=model_->nr_class; - int nr_w; - if(nr_class==2) - nr_w = 1; - else - nr_w = nr_class; - - double label=predict_values(model_, x, prob_estimates); - for(i=0;inr_feature; - int n; - const parameter& param = model_->param; - - if(model_->bias>=0) - n=nr_feature+1; - else - n=nr_feature; - int w_size = n; - FILE *fp = fopen(model_file_name,"w"); - if(fp==NULL) return -1; - - char *old_locale = strdup(setlocale(LC_ALL, NULL)); - setlocale(LC_ALL, "C"); - - int nr_w; - if(model_->nr_class==2 && model_->param.solver_type != MCSVM_CS) - nr_w=1; - else - nr_w=model_->nr_class; - - fprintf(fp, "solver_type %s\n", solver_type_table[param.solver_type]); - fprintf(fp, "nr_class %d\n", model_->nr_class); - - if(model_->label) - { - fprintf(fp, "label"); - for(i=0; inr_class; i++) - fprintf(fp, " %d", model_->label[i]); - fprintf(fp, "\n"); - } - - fprintf(fp, "nr_feature %d\n", nr_feature); - - fprintf(fp, "bias %.16g\n", model_->bias); - - fprintf(fp, "w\n"); - for(i=0; iw[i*nr_w+j]); - fprintf(fp, "\n"); - } - - setlocale(LC_ALL, old_locale); - free(old_locale); - - if (ferror(fp) != 0 || fclose(fp) != 0) return -1; - else return 0; -} - -struct model *load_model(const char *model_file_name) -{ - FILE *fp = fopen(model_file_name,"r"); - if(fp==NULL) return NULL; - - int i; - int nr_feature; - int n; - int nr_class; - double bias; - model *model_ = Malloc(model,1); - parameter& param = model_->param; - - model_->label = NULL; - - char *old_locale = strdup(setlocale(LC_ALL, NULL)); - setlocale(LC_ALL, "C"); - - char cmd[81]; - while(1) - { - fscanf(fp,"%80s",cmd); - if(strcmp(cmd,"solver_type")==0) - { - fscanf(fp,"%80s",cmd); - int i; - for(i=0;solver_type_table[i];i++) - { - if(strcmp(solver_type_table[i],cmd)==0) - { - param.solver_type=i; - break; - } - } - if(solver_type_table[i] == NULL) - { - fprintf(stderr,"unknown solver type.\n"); - - setlocale(LC_ALL, old_locale); - free(model_->label); - free(model_); - free(old_locale); - return NULL; - } - } - else if(strcmp(cmd,"nr_class")==0) - { - fscanf(fp,"%d",&nr_class); - model_->nr_class=nr_class; - } - else if(strcmp(cmd,"nr_feature")==0) - { - fscanf(fp,"%d",&nr_feature); - model_->nr_feature=nr_feature; - } - else if(strcmp(cmd,"bias")==0) - { - fscanf(fp,"%lf",&bias); - model_->bias=bias; - } - else if(strcmp(cmd,"w")==0) - { - break; - } - else if(strcmp(cmd,"label")==0) - { - int nr_class = model_->nr_class; - model_->label = Malloc(int,nr_class); - for(int i=0;ilabel[i]); - } - else - { - fprintf(stderr,"unknown text in model file: [%s]\n",cmd); - setlocale(LC_ALL, old_locale); - free(model_->label); - free(model_); - free(old_locale); - return NULL; - } - } - - nr_feature=model_->nr_feature; - if(model_->bias>=0) - n=nr_feature+1; - else - n=nr_feature; - int w_size = n; - int nr_w; - if(nr_class==2 && param.solver_type != MCSVM_CS) - nr_w = 1; - else - nr_w = nr_class; - - model_->w=Malloc(double, w_size*nr_w); - for(i=0; iw[i*nr_w+j]); - fscanf(fp, "\n"); - } - - setlocale(LC_ALL, old_locale); - free(old_locale); - - if (ferror(fp) != 0 || fclose(fp) != 0) return NULL; - - return model_; -} - -int get_nr_feature(const model *model_) -{ - return model_->nr_feature; -} - -int get_nr_class(const model *model_) -{ - return model_->nr_class; -} - -void get_labels(const model *model_, int* label) -{ - if (model_->label != NULL) - for(int i=0;inr_class;i++) - label[i] = model_->label[i]; -} - -void free_model_content(struct model *model_ptr) -{ - if(model_ptr->w != NULL) - free(model_ptr->w); - if(model_ptr->label != NULL) - free(model_ptr->label); -} - -void free_and_destroy_model(struct model **model_ptr_ptr) -{ - struct model *model_ptr = *model_ptr_ptr; - if(model_ptr != NULL) - { - free_model_content(model_ptr); - free(model_ptr); - } -} - -void destroy_param(parameter* param) -{ - if(param->weight_label != NULL) - free(param->weight_label); - if(param->weight != NULL) - free(param->weight); -} - -const char *check_parameter(const problem *prob, const parameter *param) -{ - if(param->eps <= 0) - return "eps <= 0"; - - if(param->C <= 0) - return "C <= 0"; - - if(param->p < 0) - return "p < 0"; - - if(param->solver_type != L2R_LR - && param->solver_type != L2R_L2LOSS_SVC_DUAL - && param->solver_type != L2R_L2LOSS_SVC - && param->solver_type != L2R_L1LOSS_SVC_DUAL - && param->solver_type != MCSVM_CS - && param->solver_type != L1R_L2LOSS_SVC - && param->solver_type != L1R_LR - && param->solver_type != L2R_LR_DUAL - && param->solver_type != L2R_L2LOSS_SVR - && param->solver_type != L2R_L2LOSS_SVR_DUAL - && param->solver_type != L2R_L1LOSS_SVR_DUAL) - return "unknown solver type"; - - return NULL; -} - -int check_probability_model(const struct model *model_) -{ - return (model_->param.solver_type==L2R_LR || - model_->param.solver_type==L2R_LR_DUAL || - model_->param.solver_type==L1R_LR); -} - -void set_print_string_function(void (*print_func)(const char*)) -{ - if (print_func == NULL) - liblinear_print_string = &print_string_stdout; - else - liblinear_print_string = print_func; -} - diff --git a/linear_models/linear.h b/linear_models/linear.h index 22a3567..79816c0 100644 --- a/linear_models/linear.h +++ b/linear_models/linear.h @@ -1,10 +1,14 @@ #ifndef _LIBLINEAR_H #define _LIBLINEAR_H +#define LIBLINEAR_VERSION 241 + #ifdef __cplusplus extern "C" { #endif +extern int liblinear_version; + struct feature_node { int index; @@ -16,10 +20,10 @@ struct problem int l, n; double *y; struct feature_node **x; - double bias; /* < 0 if no bias term */ + double bias; /* < 0 if no bias term */ }; -enum { L2R_LR, L2R_L2LOSS_SVC_DUAL, L2R_L2LOSS_SVC, L2R_L1LOSS_SVC_DUAL, MCSVM_CS, L1R_L2LOSS_SVC, L1R_LR, L2R_LR_DUAL, L2R_L2LOSS_SVR = 11, L2R_L2LOSS_SVR_DUAL, L2R_L1LOSS_SVR_DUAL }; /* solver_type */ +enum { L2R_LR, L2R_L2LOSS_SVC_DUAL, L2R_L2LOSS_SVC, L2R_L1LOSS_SVC_DUAL, MCSVM_CS, L1R_L2LOSS_SVC, L1R_LR, L2R_LR_DUAL, L2R_L2LOSS_SVR = 11, L2R_L2LOSS_SVR_DUAL, L2R_L1LOSS_SVR_DUAL, ONECLASS_SVM = 21 }; /* solver_type */ struct parameter { @@ -32,6 +36,9 @@ struct parameter int *weight_label; double* weight; double p; + double nu; + double *init_sol; + int regularize_bias; }; struct model @@ -42,10 +49,12 @@ struct model double *w; int *label; /* label of each class */ double bias; + double rho; /* one-class SVM only */ }; struct model* train(const struct problem *prob, const struct parameter *param); void cross_validation(const struct problem *prob, const struct parameter *param, int nr_fold, double *target); +void find_parameters(const struct problem *prob, const struct parameter *param, int nr_fold, double start_C, double start_p, double *best_C, double *best_p, double *best_score); double predict_values(const struct model *model_, const struct feature_node *x, double* dec_values); double predict(const struct model *model_, const struct feature_node *x); @@ -57,6 +66,9 @@ struct model *load_model(const char *model_file_name); int get_nr_feature(const struct model *model_); int get_nr_class(const struct model *model_); void get_labels(const struct model *model_, int* label); +double get_decfun_coef(const struct model *model_, int feat_idx, int label_idx); +double get_decfun_bias(const struct model *model_, int label_idx); +double get_decfun_rho(const struct model *model_); void free_model_content(struct model *model_ptr); void free_and_destroy_model(struct model **model_ptr_ptr); @@ -64,11 +76,12 @@ void destroy_param(struct parameter *param); const char *check_parameter(const struct problem *prob, const struct parameter *param); int check_probability_model(const struct model *model); +int check_regression_model(const struct model *model); +int check_oneclass_model(const struct model *model); void set_print_string_function(void (*print_func) (const char*)); #ifdef __cplusplus } #endif -#endif /* _LIBLINEAR_H */ - +#endif /* _LIBLINEAR_H */ \ No newline at end of file diff --git a/linear_models/linear_models_test.go b/linear_models/linear_models_test.go index 2a6cc41..f4e6615 100644 --- a/linear_models/linear_models_test.go +++ b/linear_models/linear_models_test.go @@ -1,9 +1,10 @@ package linear_models import ( + "testing" + "github.com/sjwhitworth/golearn/base" . "github.com/smartystreets/goconvey/convey" - "testing" ) func TestLogisticRegression(t *testing.T) { @@ -24,14 +25,14 @@ func TestLogisticRegression(t *testing.T) { Z, err := lr.Predict(Y) So(err, ShouldEqual, nil) Convey("The result should be 1", func() { - So(Z.RowString(0), ShouldEqual, "1.0") + So(Z.RowString(0), ShouldEqual, "1") }) }) Convey("When predicting the label of second vector", func() { Z, err := lr.Predict(Y) So(err, ShouldEqual, nil) Convey("The result should be -1", func() { - So(Z.RowString(1), ShouldEqual, "-1.0") + So(Z.RowString(1), ShouldEqual, "0") }) }) }) diff --git a/linear_models/logistic.go b/linear_models/logistic.go index 96c3206..14ff0d2 100644 --- a/linear_models/logistic.go +++ b/linear_models/logistic.go @@ -3,6 +3,7 @@ package linear_models import ( "errors" "fmt" + "github.com/sjwhitworth/golearn/base" ) diff --git a/linear_models/logistic_test.go b/linear_models/logistic_test.go index a07f961..dae9993 100644 --- a/linear_models/logistic_test.go +++ b/linear_models/logistic_test.go @@ -1,9 +1,10 @@ package linear_models import ( + "testing" + "github.com/sjwhitworth/golearn/base" . "github.com/smartystreets/goconvey/convey" - "testing" ) func TestLogistic(t *testing.T) { @@ -23,14 +24,14 @@ func TestLogistic(t *testing.T) { Z, err := lr.Predict(Y) So(err, ShouldEqual, nil) Convey("The result should be 1", func() { - So(Z.RowString(0), ShouldEqual, "-1.0") + So(Z.RowString(0), ShouldEqual, "1") }) }) Convey("When predicting the label of second vector", func() { Z, err := lr.Predict(Y) So(err, ShouldEqual, nil) - Convey("The result should be -1", func() { - So(Z.RowString(1), ShouldEqual, "-1.0") + Convey("The result should be 2", func() { + So(Z.RowString(1), ShouldEqual, "0") }) }) So((*lr).String(), ShouldEqual, "LogisticRegression") diff --git a/linear_models/test.csv b/linear_models/test.csv index d2ac5e5..1d82045 100644 --- a/linear_models/test.csv +++ b/linear_models/test.csv @@ -1,2 +1,2 @@ -1.0,1.0,0.0,0.0,1.0 -0.0,0.0,1.0,1.0,-1.0 \ No newline at end of file +6,148,72,35,0,33.6,0.627,50,1 +1,85,66,29,0,26.6,0.351,31,0 \ No newline at end of file diff --git a/linear_models/tmp b/linear_models/tmp deleted file mode 100644 index 4560a48d7ba89141002965bdd54fed1b020adaaf..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 484 zcmV`}MMWHrvPG zp_|THfkj6aI=GGnvYm(ir-sQ1ln&Hr=()#$kmTV93aVzgXJFk2$$ z@zESQNC5x<0RR8;^>=agp|Kr=>whEE{BK}xW@b8?|LFmW!MQ)aq_QAY!N({l-Y1AF zuP7c|n=2S|Gh<|r5{7%3Qlgws+JOG=AU6^yx(G82mxbPc)6xpWQm3=PbUj7`i; zEi4TUEG>;q6^IoBX*M!6F)_0=HZ->|GPJZXH&h^2jLSd)1i1z}1cx>N00030{~E=k aU=)m800000|Nj910RR7Vw`K|e6aWAL0^^JT diff --git a/linear_models/train.csv b/linear_models/train.csv index 6c0a1f1..40895d4 100644 --- a/linear_models/train.csv +++ b/linear_models/train.csv @@ -1,4 +1,21 @@ -0.0, 0.0, 0.0, 1.0, -1.0 -0.0, 0.0, 1.0, 0.0, -1.0 -0.0, 1.0, 0.0, 0.0, 1.0 -1.0, 0.0, 0.0, 0.0, 1.0 \ No newline at end of file +6,148,72,35,0,33.6,0.627,50,1 +1,85,66,29,0,26.6,0.351,31,0 +8,183,64,0,0,23.3,0.672,32,1 +1,89,66,23,94,28.1,0.167,21,0 +0,137,40,35,168,43.1,2.288,33,1 +5,116,74,0,0,25.6,0.201,30,0 +3,78,50,32,88,31,0.248,26,1 +10,115,0,0,0,35.3,0.134,29,0 +2,197,70,45,543,30.5,0.158,53,1 +8,125,96,0,0,0,0.232,54,1 +4,110,92,0,0,37.6,0.191,30,0 +10,168,74,0,0,38,0.537,34,1 +10,139,80,0,0,27.1,1.441,57,0 +1,189,60,23,846,30.1,0.398,59,1 +5,166,72,19,175,25.8,0.587,51,1 +7,100,0,0,0,30,0.484,32,1 +0,118,84,47,230,45.8,0.551,31,1 +7,107,74,0,0,29.6,0.254,31,1 +1,103,30,38,83,43.3,0.183,33,0 +1,115,70,30,96,34.6,0.529,32,1 +3,126,88,41,235,39.3,0.704,27,0