From 1e340d509f229120eb3aaa98c91028dc3c0d3305 Mon Sep 17 00:00:00 2001 From: Gertjan van den Burg Date: Mon, 25 Aug 2014 14:38:03 +0200 Subject: rename msvmmaj to gensvm --- src/gensvm_kernel.c | 412 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 412 insertions(+) create mode 100644 src/gensvm_kernel.c (limited to 'src/gensvm_kernel.c') diff --git a/src/gensvm_kernel.c b/src/gensvm_kernel.c new file mode 100644 index 0000000..55cfa03 --- /dev/null +++ b/src/gensvm_kernel.c @@ -0,0 +1,412 @@ +/** + * @file gensvm_kernel.c + * @author Gertjan van den Burg + * @date October 18, 2013 + * @brief Defines main functions for use of kernels in GenSVM. + * + * @details + * Functions for constructing different kernels using user-supplied + * parameters. Also contains the functions for decomposing the + * kernel matrix using several decomposition methods. + * + */ + +#include + +#include "gensvm.h" +#include "gensvm_kernel.h" +#include "gensvm_lapack.h" +#include "gensvm_matrix.h" +#include "util.h" + +/** + * @brief Create the kernel matrix + * + * Create a kernel matrix based on the specified kerneltype. Kernel parameters + * are assumed to be specified in the model. + * + * @param[in] model GenModel specifying the parameters + * @param[in] data GenData specifying the data. + * + */ +void gensvm_make_kernel(struct GenModel *model, struct GenData *data) +{ + long i, j; + // Determine if a kernel needs to be computed. This is not the case if + // a LINEAR kernel is requested in the model, or if the requested + // kernel is already in the data. + + if (model->kerneltype == K_LINEAR) { + data->J = Calloc(double, data->m+1); + for (i=1; im+1; i++) { + matrix_set(data->J, 1, i, 0, 1.0); + } + return; + } + + /* + switch (model->kerneltype) { + case K_LINEAR: + // if data has another kernel, free that matrix and + // assign Z to RAW + if (data->kerneltype != K_LINEAR) { + free(data->Z); + data->Z = data->RAW; + } + data->J = Calloc(double, data->m+1); + for (i=1; im+1; i++) { + matrix_set(data->J, 1, i, 0, 1.0); + } + return; + case K_POLY: + // if data has another kernel, we need to recalculate + if (data->kerneltype != K_POLY) { + break; + } + // if it is poly, we only recalculate if the kernel + // parameters differ + if (data->kernelparam[0] == model->kernelparam[0] && + data->kernelparam[1] == model->kernelparam[1] && + data->kernelparam[2] == model->kernelparam[2]) + // < do something with J ? + return; + case K_RBF: + if (data->kerneltype != K_RBF) + break; + if (data->kernelparam[0] == model->kernelparam[0]) + // < do something with J ? + return; + case K_SIGMOID: + if (data->kerneltype != K_SIGMOID) + break; + if (data->kernelparam[0] == model->kernelparam[0] && + data->kernelparam[1] == model->kernelparam[1]) + // < do something with J ? + return; + } + */ + long n = data->n; + double value; + double *x1, *x2; + double *K = Calloc(double, n*n); + + for (i=0; iRAW[i*(data->m+1)+1]; + x2 = &data->RAW[j*(data->m+1)+1]; + if (model->kerneltype == K_POLY) + value = gensvm_compute_poly(x1, x2, + model->kernelparam, data->m); + else if (model->kerneltype == K_RBF) + value = gensvm_compute_rbf(x1, x2, + model->kernelparam, data->m); + else if (model->kerneltype == K_SIGMOID) + value = gensvm_compute_sigmoid(x1, x2, + model->kernelparam, data->m); + else { + fprintf(stderr, "Unknown kernel type in " + "gensvm_make_kernel\n"); + exit(1); + } + matrix_set(K, n, i, j, value); + matrix_set(K, n, j, i, value); + } + } + + double *P = NULL; + double *Sigma = NULL; + long num_eigen = gensvm_make_eigen(K, n, &P, &Sigma); + //printf("num eigen: %li\n", num_eigen); + data->m = num_eigen; + + // copy eigendecomp to data + data->Z = Calloc(double, n*(num_eigen+1)); + for (i=0; iZ, num_eigen+1, i, j, value); + } + matrix_set(data->Z, num_eigen+1, i, 0, 1.0); + } + + // Set the regularization matrix (change if not full rank used) + if (data->J != NULL) + free(data->J); + data->J = Calloc(double, data->m+1); + for (i=1; im+1; i++) { + value = 1.0/matrix_get(Sigma, 1, i-1, 0); + matrix_set(data->J, 1, i, 0, value); + } + + // let data know what it's made of + data->kerneltype = model->kerneltype; + free(data->kernelparam); + switch (model->kerneltype) { + case K_LINEAR: + break; + case K_POLY: + data->kernelparam = Calloc(double, 3); + data->kernelparam[0] = model->kernelparam[0]; + data->kernelparam[1] = model->kernelparam[1]; + data->kernelparam[2] = model->kernelparam[2]; + break; + case K_RBF: + data->kernelparam = Calloc(double, 1); + data->kernelparam[0] = model->kernelparam[0]; + break; + case K_SIGMOID: + data->kernelparam = Calloc(double, 2); + data->kernelparam[0] = model->kernelparam[0]; + data->kernelparam[1] = model->kernelparam[1]; + } + free(K); + free(Sigma); + free(P); +} + +/** + * @brief Find the (reduced) eigendecomposition of a kernel matrix. + * + * @details. + * tbd + * + */ +long gensvm_make_eigen(double *K, long n, double **P, double **Sigma) +{ + int M, status, LWORK, *IWORK, *IFAIL; + long i, j, num_eigen, cutoff_idx; + double max_eigen, abstol, *WORK; + + double *tempSigma = Malloc(double, n); + double *tempP = Malloc(double, n*n); + + IWORK = Malloc(int, 5*n); + IFAIL = Malloc(int, n); + + // highest precision eigenvalues, may reduce for speed + abstol = 2.0*dlamch('S'); + + // first perform a workspace query to determine optimal size of the + // WORK array. + WORK = Malloc(double, 1); + status = dsyevx( + 'V', + 'A', + 'U', + n, + K, + n, + 0, + 0, + 0, + 0, + abstol, + &M, + tempSigma, + tempP, + n, + WORK, + -1, + IWORK, + IFAIL); + LWORK = WORK[0]; + + // allocate the requested memory for the eigendecomposition + WORK = (double *)realloc(WORK, LWORK*sizeof(double)); + status = dsyevx( + 'V', + 'A', + 'U', + n, + K, + n, + 0, + 0, + 0, + 0, + abstol, + &M, + tempSigma, + tempP, + n, + WORK, + LWORK, + IWORK, + IFAIL); + + if (status != 0) { + fprintf(stderr, "Nonzero exit status from dsyevx. Exiting..."); + exit(1); + } + + // Select the desired number of eigenvalues, depending on their size. + // dsyevx sorts eigenvalues in ascending order. + // + max_eigen = tempSigma[n-1]; + cutoff_idx = 0; + + for (i=0; i 1e-10 ) { + cutoff_idx = i; + break; + } + + num_eigen = n - cutoff_idx; + + *Sigma = Calloc(double, num_eigen); + + for (i=0; in-1-num_eigen; j--) { + for (i=0; in; + long n_test = data_test->n; + long m = data_test->m; + double value; + double *x1, *x2; + + *K2 = Calloc(double, n_test*n_train); + + //printf("Training RAW\n"); + //print_matrix(data_train->RAW, n_train, m+1); + + //printf("Testing RAW\n"); + //print_matrix(data_test->RAW, n_test, m+1); + + for (i=0; iRAW[i*(m+1)+1]; + x2 = &data_train->RAW[j*(m+1)+1]; + if (model->kerneltype == K_POLY) + value = gensvm_compute_poly(x1, x2, + model->kernelparam, + m); + else if (model->kerneltype == K_RBF) + value = gensvm_compute_rbf(x1, x2, + model->kernelparam, + m); + else if (model->kerneltype == K_SIGMOID) + value = gensvm_compute_sigmoid(x1, x2, + model->kernelparam, + m); + else { + fprintf(stderr, "Unknown kernel type in " + "gensvm_make_crosskernel\n"); + exit(1); + } + matrix_set((*K2), n_train, i, j, value); + } + } + + //printf("cross K2:\n"); + //print_matrix((*K2), n_test, n_train); + +} + +/** + * @brief Compute the RBF kernel between two vectors + * + * @details + * The RBF kernel is computed between two vectors. This kernel is defined as + * @f[ + * k(x_1, x_2) = \exp( -\gamma \| x_1 - x_2 \|^2 ) + * @f] + * where @f$ \gamma @f$ is a kernel parameter specified. + * + * @param[in] x1 first vector + * @param[in] x2 second vector + * @param[in] kernelparam array of kernel parameters (gamma is first + * element) + * @param[in] n length of the vectors x1 and x2 + * @returns kernel evaluation + */ +double gensvm_compute_rbf(double *x1, double *x2, double *kernelparam, long n) +{ + long i; + double value = 0.0; + + for (i=0; i