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/crossval.c | 24 +- src/gensvm_init.c | 282 +++++++++++++++++ src/gensvm_io.c | 295 +++++++++++++++++ src/gensvm_kernel.c | 412 ++++++++++++++++++++++++ src/gensvm_lapack.c | 135 ++++++++ src/gensvm_matrix.c | 38 +++ src/gensvm_pred.c | 215 +++++++++++++ src/gensvm_sv.c | 45 +++ src/gensvm_train.c | 534 +++++++++++++++++++++++++++++++ src/gensvm_train_dataset.c | 748 ++++++++++++++++++++++++++++++++++++++++++++ src/libGenSVM.c | 326 +++++++++++++++++++ src/libMSVMMaj.c | 326 ------------------- src/msvmmaj_init.c | 282 ----------------- src/msvmmaj_io.c | 295 ----------------- src/msvmmaj_kernel.c | 412 ------------------------ src/msvmmaj_lapack.c | 135 -------- src/msvmmaj_matrix.c | 38 --- src/msvmmaj_pred.c | 215 ------------- src/msvmmaj_sv.c | 45 --- src/msvmmaj_train.c | 534 ------------------------------- src/msvmmaj_train_dataset.c | 748 -------------------------------------------- src/predGenSVM.c | 180 +++++++++++ src/predMSVMMaj.c | 180 ----------- src/trainGenSVM.c | 244 +++++++++++++++ src/trainGenSVMdataset.c | 321 +++++++++++++++++++ src/trainMSVMMaj.c | 244 --------------- src/trainMSVMMajdataset.c | 321 ------------------- src/util.c | 24 +- 28 files changed, 3799 insertions(+), 3799 deletions(-) create mode 100644 src/gensvm_init.c create mode 100644 src/gensvm_io.c create mode 100644 src/gensvm_kernel.c create mode 100644 src/gensvm_lapack.c create mode 100644 src/gensvm_matrix.c create mode 100644 src/gensvm_pred.c create mode 100644 src/gensvm_sv.c create mode 100644 src/gensvm_train.c create mode 100644 src/gensvm_train_dataset.c create mode 100644 src/libGenSVM.c delete mode 100644 src/libMSVMMaj.c delete mode 100644 src/msvmmaj_init.c delete mode 100644 src/msvmmaj_io.c delete mode 100644 src/msvmmaj_kernel.c delete mode 100644 src/msvmmaj_lapack.c delete mode 100644 src/msvmmaj_matrix.c delete mode 100644 src/msvmmaj_pred.c delete mode 100644 src/msvmmaj_sv.c delete mode 100644 src/msvmmaj_train.c delete mode 100644 src/msvmmaj_train_dataset.c create mode 100644 src/predGenSVM.c delete mode 100644 src/predMSVMMaj.c create mode 100644 src/trainGenSVM.c create mode 100644 src/trainGenSVMdataset.c delete mode 100644 src/trainMSVMMaj.c delete mode 100644 src/trainMSVMMajdataset.c (limited to 'src') diff --git a/src/crossval.c b/src/crossval.c index 1b5a592..85f9341 100644 --- a/src/crossval.c +++ b/src/crossval.c @@ -6,16 +6,16 @@ * * @details * This file contains functions for performing cross validation. The funtion - * msvmmaj_make_cv_split() creates a cross validation vector for non-stratified - * cross validation. The function msvmmaj_get_tt_split() creates a train and + * gensvm_make_cv_split() creates a cross validation vector for non-stratified + * cross validation. The function gensvm_get_tt_split() creates a train and * test dataset from a given dataset and a pre-determined CV partition vector. * See individual function documentation for details. * */ #include "crossval.h" -#include "msvmmaj.h" -#include "msvmmaj_matrix.h" +#include "gensvm.h" +#include "gensvm_matrix.h" /** * @brief Create a cross validation split vector @@ -35,7 +35,7 @@ * for each observation on exit * */ -void msvmmaj_make_cv_split(long N, long folds, long *cv_idx) +void gensvm_make_cv_split(long N, long folds, long *cv_idx) { long i, j, idx; @@ -71,23 +71,23 @@ void msvmmaj_make_cv_split(long N, long folds, long *cv_idx) * @brief Create train and test datasets for a CV split * * @details - * Given a MajData structure for the full dataset, a previously created + * Given a GenData structure for the full dataset, a previously created * cross validation split vector and a fold index, a training and test dataset * are created. * - * @param[in] full_data a MajData structure for the entire + * @param[in] full_data a GenData structure for the entire * dataset - * @param[in,out] train_data an initialized MajData structure which + * @param[in,out] train_data an initialized GenData structure which * on exit contains the training dataset - * @param[in,out] test_data an initialized MajData structure which + * @param[in,out] test_data an initialized GenData structure which * on exit contains the test dataset * @param[in] cv_idx a vector of cv partitions created by - * msvmmaj_make_cv_split() + * gensvm_make_cv_split() * @param[in] fold_idx index of the fold which becomes the * test dataset */ -void msvmmaj_get_tt_split(struct MajData *full_data, struct MajData *train_data, - struct MajData *test_data, long *cv_idx, long fold_idx) +void gensvm_get_tt_split(struct GenData *full_data, struct GenData *train_data, + struct GenData *test_data, long *cv_idx, long fold_idx) { long i, j, k, l, test_n, train_n; diff --git a/src/gensvm_init.c b/src/gensvm_init.c new file mode 100644 index 0000000..b3f214e --- /dev/null +++ b/src/gensvm_init.c @@ -0,0 +1,282 @@ +/** + * @file gensvm_init.c + * @author Gertjan van den Burg + * @date January 7, 2014 + * @brief Functions for initializing model and data structures + * + * @details + * This file contains functions for initializing a GenModel instance + * and a GenData instance. In addition, default values for these + * structures are defined here (and only here). Functions for allocating + * memory for the model structure and freeing of the model and data structures + * are also included. + * + */ + +#include + +#include "gensvm.h" +#include "gensvm_init.h" + +/** + * @brief Initialize a GenModel structure + * + * @details + * A GenModel structure is initialized and the default value for the + * parameters are set. A pointer to the initialized model is returned. + * + * @returns initialized GenModel + */ +struct GenModel *gensvm_init_model() +{ + struct GenModel *model = Malloc(struct GenModel, 1); + + // set default values + model->p = 1.0; + model->lambda = pow(2, -8.0); + model->epsilon = 1e-6; + model->kappa = 0.0; + model->weight_idx = 1; + model->kerneltype = K_LINEAR; + model->kernelparam = NULL; + + model->W = NULL; + model->t = NULL; + model->V = NULL; + model->Vbar = NULL; + model->U = NULL; + model->UU = NULL; + model->Q = NULL; + model->H = NULL; + model->R = NULL; + model->rho = NULL; + model->data_file = NULL; + + return model; +} + +/** + * @brief Initialize a GenData structure + * + * @details + * A GenData structure is initialized and default values are set. + * A pointer to the initialized data is returned. + * + * @returns initialized GenData + * + */ +struct GenData *gensvm_init_data() +{ + struct GenData *data = Malloc(struct GenData, 1); + data->J = NULL; + data->y = NULL; + data->Z = NULL; + data->RAW = NULL; + + // set default values + data->kerneltype = K_LINEAR; + data->kernelparam = NULL; + + return data; +} + +/** + * @brief Allocate memory for a GenModel + * + * @details + * This function can be used to allocate the memory needed for a GenModel. All + * arrays in the model are specified and initialized to 0. + * + * @param[in] model GenModel to allocate + * + */ +void gensvm_allocate_model(struct GenModel *model) +{ + long n = model->n; + long m = model->m; + long K = model->K; + + model->W = Calloc(double, m*(K-1)); + if (model->W == NULL) { + fprintf(stderr, "Failed to allocate memory for W.\n"); + exit(1); + } + + model->t = Calloc(double, K-1); + if (model->t == NULL) { + fprintf(stderr, "Failed to allocate memory for t.\n"); + exit(1); + } + + model->V = Calloc(double, (m+1)*(K-1)); + if (model->V == NULL) { + fprintf(stderr, "Failed to allocate memory for V.\n"); + exit(1); + } + + model->Vbar = Calloc(double, (m+1)*(K-1)); + if (model->Vbar == NULL) { + fprintf(stderr, "Failed to allocate memory for Vbar.\n"); + exit(1); + } + + model->U = Calloc(double, K*(K-1)); + if (model->U == NULL) { + fprintf(stderr, "Failed to allocate memory for U.\n"); + exit(1); + } + + model->UU = Calloc(double, n*K*(K-1)); + if (model->UU == NULL) { + fprintf(stderr, "Failed to allocate memory for UU.\n"); + exit(1); + } + + model->Q = Calloc(double, n*K); + if (model->Q == NULL) { + fprintf(stderr, "Failed to allocate memory for Q.\n"); + exit(1); + } + + model->H = Calloc(double, n*K); + if (model->H == NULL) { + fprintf(stderr, "Failed to allocate memory for H.\n"); + exit(1); + } + + model->R = Calloc(double, n*K); + if (model->R == NULL) { + fprintf(stderr, "Failed to allocate memory for R.\n"); + exit(1); + } + + model->rho = Calloc(double, n); + if (model->rho == NULL) { + fprintf(stderr, "Failed to allocate memory for rho.\n"); + exit(1); + } +} + +/** + * @brief Reallocate memory for GenModel + * + * @details + * This function can be used to reallocate existing memory for a GenModel, + * upon a change in the model dimensions. This is used in combination with + * kernels. + * + * @param[in] model GenModel to reallocate + * @param[in] n new value of GenModel->n + * @param[in] m new value of GenModel->m + * + */ +void gensvm_reallocate_model(struct GenModel *model, long n, long m) +{ + long K = model->K; + + if (model->n == n && model->m == m) + return; + if (model->n != n) { + model->UU = (double *) realloc(model->UU, + n*K*(K-1)*sizeof(double)); + if (model->UU == NULL) { + fprintf(stderr, "Failed to reallocate UU\n"); + exit(1); + } + + model->Q = (double *) realloc(model->Q, n*K*sizeof(double)); + if (model->Q == NULL) { + fprintf(stderr, "Failed to reallocate Q\n"); + exit(1); + } + + model->H = (double *) realloc(model->H, n*K*sizeof(double)); + if (model->H == NULL) { + fprintf(stderr, "Failed to reallocate H\n"); + exit(1); + } + + model->R = (double *) realloc(model->R, n*K*sizeof(double)); + if (model->R == NULL) { + fprintf(stderr, "Failed to reallocate R\n"); + exit(1); + } + + model->rho = (double *) realloc(model->rho, n*sizeof(double)); + if (model->rho == NULL) { + fprintf(stderr, "Failed to reallocte rho\n"); + exit(1); + } + + model->n = n; + } + if (model->m != m) { + model->W = (double *) realloc(model->W, + m*(K-1)*sizeof(double)); + if (model->W == NULL) { + fprintf(stderr, "Failed to reallocate W\n"); + exit(1); + } + + model->V = (double *) realloc(model->V, + (m+1)*(K-1)*sizeof(double)); + if (model->V == NULL) { + fprintf(stderr, "Failed to reallocate V\n"); + exit(1); + } + + model->Vbar = (double *) realloc(model->Vbar, + (m+1)*(K-1)*sizeof(double)); + if (model->Vbar == NULL) { + fprintf(stderr, "Failed to reallocate Vbar\n"); + exit(1); + } + + model->m = m; + } +} + +/** + * @brief Free allocated GenModel struct + * + * @details + * Simply free a previously allocated GenModel by freeing all its component + * arrays. Note that the model struct itself is also freed here. + * + * @param[in] model GenModel to free + * + */ +void gensvm_free_model(struct GenModel *model) +{ + free(model->W); + free(model->t); + free(model->V); + free(model->Vbar); + free(model->U); + free(model->UU); + free(model->Q); + free(model->H); + free(model->rho); + free(model->R); + free(model->kernelparam); + + free(model); +} + +/** + * @brief Free allocated GenData struct + * + * @details + * Simply free a previously allocated GenData struct by freeing all its + * components. Note that the data struct itself is also freed here. + * + * @param[in] data GenData struct to free + * + */ +void gensvm_free_data(struct GenData *data) +{ + free(data->Z); + free(data->y); + free(data->J); + free(data); +} diff --git a/src/gensvm_io.c b/src/gensvm_io.c new file mode 100644 index 0000000..546ecd5 --- /dev/null +++ b/src/gensvm_io.c @@ -0,0 +1,295 @@ +/** + * @file gensvm_io.c + * @author Gertjan van den Burg + * @date January, 2014 + * @brief Functions for input and output of data and model files + * + * @details + * This file contains functions for reading and writing model files, and data + * files. + * + */ + +#include "gensvm.h" +#include "gensvm_io.h" +#include "gensvm_matrix.h" +#include "strutil.h" +#include "timer.h" + +/** + * @brief Read data from file + * + * @details + * Read the data from the data_file. The data matrix X is augmented + * with a column of ones, to get the matrix Z. The data is expected + * to follow a specific format, which is specified in the @ref spec_data_file. + * The class labels are corrected internally to correspond to the interval + * [1 .. K], where K is the total number of classes. + * + * @todo + * Make sure that this function allows datasets without class labels for + * testing. + * + * @param[in,out] dataset initialized GenData struct + * @param[in] data_file filename of the data file. + */ +void gensvm_read_data(struct GenData *dataset, char *data_file) +{ + FILE *fid; + long i, j; + long n, m; // dimensions of data + long nr = 0; // used to check consistency of data + double value; + long K = 0; + long min_y = 1000000; + + char buf[MAX_LINE_LENGTH]; + + if ((fid = fopen(data_file, "r")) == NULL) { + fprintf(stderr, "\nERROR: datafile %s could not be opened.\n", + data_file); + exit(0); + } + + // Read data dimensions + nr += fscanf(fid, "%ld", &n); + nr += fscanf(fid, "%ld", &m); + + // Allocate memory + dataset->RAW = Malloc(double, n*(m+1)); + + // Read first line of data + for (j=1; jRAW, n, 0, j, value); + } + + // Check if there is a label at the end of the line + if (fgets(buf, MAX_LINE_LENGTH, fid) == NULL) { + fprintf(stderr, "ERROR: No label found on first line.\n"); + exit(1); + } + if (sscanf(buf, "%lf", &value) > 0) { + dataset->y = Malloc(long, n); + dataset->y[0] = value; + } else if (dataset->y != NULL) { + free(dataset->y); + dataset->y = NULL; + } + + // Read the rest of the file + for (i=1; iRAW, m+1, i, j, value); + } + if (dataset->y != NULL) { + nr += fscanf(fid, "%lf", &value); + dataset->y[i] = (long) value; + K = maximum(K, value); + min_y = minimum(min_y, value); + } + } + fclose(fid); + + // Correct labels: must be in [1, K] + if (min_y == 0) { + for (i=0; iy[i]++; + K++; + } else if (min_y < 0 ) { + fprintf(stderr, "ERROR: wrong class labels in %s, minimum " + "value is: %ld\n", + data_file, min_y); + exit(0); + } + + if (nr < n * m) { + fprintf(stderr, "ERROR: not enough data found in %s\n", + data_file); + exit(0); + } + + // Set the column of ones + for (i=0; iRAW, m+1, i, 0, 1.0); + + dataset->n = n; + dataset->m = m; + dataset->K = K; + dataset->Z = dataset->RAW; +} + + +/** + * @brief Read model from file + * + * @details + * Read a GenModel from a model file. The GenModel struct must have been + * initalized elswhere. The model file is expected to follow the @ref + * spec_model_file. The easiest way to generate a model file is through + * gensvm_write_model(), which can for instance be used in trainGenSVM.c. + * + * @param[in,out] model initialized GenModel + * @param[in] model_filename filename of the model file + * + */ +void gensvm_read_model(struct GenModel *model, char *model_filename) +{ + long i, j, nr = 0; + FILE *fid; + char buffer[MAX_LINE_LENGTH]; + char data_filename[MAX_LINE_LENGTH]; + double value = 0; + + fid = fopen(model_filename, "r"); + if (fid == NULL) { + fprintf(stderr, "Error opening model file %s\n", + model_filename); + exit(1); + } + // skip the first four lines + for (i=0; i<4; i++) + next_line(fid, model_filename); + + // read all model variables + model->p = get_fmt_double(fid, model_filename, "p = %lf"); + model->lambda = get_fmt_double(fid, model_filename, "lambda = %lf"); + model->kappa = get_fmt_double(fid, model_filename, "kappa = %lf"); + model->epsilon = get_fmt_double(fid, model_filename, "epsilon = %lf"); + model->weight_idx = (int) get_fmt_long(fid, model_filename, + "weight_idx = %li"); + + // skip to data section + for (i=0; i<2; i++) + next_line(fid, model_filename); + + // read filename of data file + if (fgets(buffer, MAX_LINE_LENGTH, fid) == NULL) { + fprintf(stderr, "Error reading model file %s\n", + model_filename); + exit(1); + } + sscanf(buffer, "filename = %s\n", data_filename); + model->data_file = data_filename; + + // read all data variables + model->n = get_fmt_long(fid, model_filename, "n = %li\n"); + model->m = get_fmt_long(fid, model_filename, "m = %li\n"); + model->K = get_fmt_long(fid, model_filename, "K = %li\n"); + + // skip to output + for (i=0; i<2; i++) + next_line(fid, model_filename); + + // read the matrix V and check for consistency + model->V = Malloc(double, (model->m+1)*(model->K-1)); + for (i=0; im+1; i++) { + for (j=0; jK-1; j++) { + nr += fscanf(fid, "%lf ", &value); + matrix_set(model->V, model->K-1, i, j, value); + } + } + if (nr != (model->m+1)*(model->K-1)) { + fprintf(stderr, "Error reading model file %s. " + "Not enough elements of V found.\n", + model_filename); + exit(1); + } +} + +/** + * @brief Write model to file + * + * @details + * Write a GenModel to a file. The current time is specified in the file in + * UTC + offset. The model file further corresponds to the @ref + * spec_model_file. + * + * @param[in] model GenModel which contains an estimate for + * GenModel::V + * @param[in] output_filename the output file to write the model to + * + */ +void gensvm_write_model(struct GenModel *model, char *output_filename) +{ + FILE *fid; + long i, j; + char timestr[MAX_LINE_LENGTH]; + + // open output file + fid = fopen(output_filename, "w"); + if (fid == NULL) { + fprintf(stderr, "Error opening output file %s", + output_filename); + exit(1); + } + get_time_string(timestr); + + // Write output to file + fprintf(fid, "Output file for GenSVM (version %1.1f)\n", VERSION); + fprintf(fid, "Generated on: %s\n\n", timestr); + fprintf(fid, "Model:\n"); + fprintf(fid, "p = %15.16f\n", model->p); + fprintf(fid, "lambda = %15.16f\n", model->lambda); + fprintf(fid, "kappa = %15.16f\n", model->kappa); + fprintf(fid, "epsilon = %g\n", model->epsilon); + fprintf(fid, "weight_idx = %i\n", model->weight_idx); + fprintf(fid, "\n"); + fprintf(fid, "Data:\n"); + fprintf(fid, "filename = %s\n", model->data_file); + fprintf(fid, "n = %li\n", model->n); + fprintf(fid, "m = %li\n", model->m); + fprintf(fid, "K = %li\n", model->K); + fprintf(fid, "\n"); + fprintf(fid, "Output:\n"); + for (i=0; im+1; i++) { + for (j=0; jK-1; j++) { + fprintf(fid, "%+15.16f ", + matrix_get(model->V, + model->K-1, i, j)); + } + fprintf(fid, "\n"); + } + + fclose(fid); +} + +/** + * @brief Write predictions to file + * + * @details + * Write the given predictions to an output file, such that the resulting file + * corresponds to the @ref spec_data_file. + * + * @param[in] data GenData with the original instances + * @param[in] predy predictions of the class labels of the + * instances in the given GenData. Note that the + * order of the instances is assumed to be the + * same. + * @param[in] output_filename the file to which the predictions are written + * + */ +void gensvm_write_predictions(struct GenData *data, long *predy, + char *output_filename) +{ + long i, j; + FILE *fid; + + fid = fopen(output_filename, "w"); + if (fid == NULL) { + fprintf(stderr, "Error opening output file %s", + output_filename); + exit(1); + } + + for (i=0; in; i++) { + for (j=0; jm; j++) + fprintf(fid, "%f ", + matrix_get(data->Z, + data->m+1, i, j+1)); + fprintf(fid, "%li\n", predy[i]); + } + + fclose(fid); +} 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; i0: if i, the leading minor of A + * was not positive definite + * + * See the LAPACK documentation at: + * http://www.netlib.org/lapack/explore-html/dc/de9/group__double_p_osolve.html + */ +int dposv(char UPLO, int N, int NRHS, double *A, int LDA, double *B, + int LDB) +{ + extern void dposv_(char *UPLO, int *Np, int *NRHSp, double *A, + int *LDAp, double *B, int *LDBp, int *INFOp); + int INFO; + dposv_(&UPLO, &N, &NRHS, A, &LDA, B, &LDB, &INFO); + return INFO; +} + +/** + * @brief Solve a system of equations AX = B where A is symmetric. + * + * @details + * Solve a linear system of equations AX = B where A is symmetric. This + * function uses the external LAPACK routine dsysv. + * + * @param[in] UPLO which triangle of A is stored + * @param[in] N order of A + * @param[in] NRHS number of columns of B + * @param[in,out] A double precision array of size (LDA, N). On + * exit contains the block diagonal matrix D and + * the multipliers used to obtain the factor U or + * L from the factorization A = U*D*U**T or + * A = L*D*L**T. + * @param[in] LDA leading dimension of A + * @param[in] IPIV integer array containing the details of D + * @param[in,out] B double precision array of size (LDB, NRHS). On + * exit contains the N-by-NRHS matrix X + * @param[in] LDB leading dimension of B + * @param[out] WORK double precision array of size max(1,LWORK). On + * exit, WORK(1) contains the optimal LWORK + * @param[in] LWORK the length of WORK, can be used for determining + * the optimal blocksize for dsystrf. + * @returns info parameter which contains the status of the + * computation: + * - =0: success + * - <0: if -i, the i-th argument had an + * illegal value + * - >0: if i, D(i, i) is exactly zero, + * no solution can be computed. + * + * See the LAPACK documentation at: + * http://www.netlib.org/lapack/explore-html/d6/d0e/group__double_s_ysolve.html + */ +int dsysv(char UPLO, int N, int NRHS, double *A, int LDA, int *IPIV, + double *B, int LDB, double *WORK, int LWORK) +{ + extern void dsysv_(char *UPLO, int *Np, int *NRHSp, double *A, + int *LDAp, int *IPIV, double *B, int *LDBp, + double *WORK, int *LWORK, int *INFOp); + int INFO; + dsysv_(&UPLO, &N, &NRHS, A, &LDA, IPIV, B, &LDB, WORK, &LWORK, &INFO); + return INFO; +} + +/** + * @brief Compute the eigenvalues and optionally the eigenvectors of a + * symmetric matrix. + * + * @details + * See the LAPACK documentation at: + * http://www.netlib.org/lapack/explore-html/d2/d97/dsyevx_8f.html + * + * + */ +int dsyevx(char JOBZ, char RANGE, char UPLO, int N, double *A, int LDA, + double VL, double VU, int IL, int IU, double ABSTOL, int *M, + double *W, double *Z, int LDZ, double *WORK, int LWORK, + int *IWORK, int *IFAIL) +{ + extern void dsyevx_(char *JOBZ, char *RANGE, char *UPLO, int *Np, + double *A, int *LDAp, double *VLp, double *VUp, + int *ILp, int *IUp, double *ABSTOLp, int *M, + double *W, double *Z, int *LDZp, double *WORK, + int *LWORKp, int *IWORK, int *IFAIL, int *INFOp); + int INFO; + dsyevx_(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, + M, W, Z, &LDZ, WORK, &LWORK, IWORK, IFAIL, &INFO); + return INFO; +} + +/** + * @brief Determine double precision machine parameters. + * + * @details + * See the LAPACK documentation at: + * http://www.netlib.org/lapack/explore-html/d5/dd4/dlamch_8f.html + */ +double dlamch(char CMACH) +{ + extern double dlamch_(char *CMACH); + return dlamch_(&CMACH); +} diff --git a/src/gensvm_matrix.c b/src/gensvm_matrix.c new file mode 100644 index 0000000..43f284f --- /dev/null +++ b/src/gensvm_matrix.c @@ -0,0 +1,38 @@ +/** + * @file gensvm_matrix.c + * @author Gertjan van den Burg + * @date August, 2013 + * @brief Functions facilitating matrix access + * + * @details + * The functions contained in this file are used when + * accessing or writing to matrices. Seperate functions + * exist of adding and multiplying existing matrix + * elements, to ensure this is done in place. + * + */ + +#include "gensvm_matrix.h" +#include "util.h" + +/** + * @brief print a matrix + * + * @details + * Debug function to print a matrix + * + * @param[in] M matrix + * @param[in] rows number of rows of M + * @param[in] cols number of columns of M + */ +void print_matrix(double *M, long rows, long cols) +{ + long i, j; + + for (i=0; i + +#include "libGenSVM.h" +#include "gensvm.h" +#include "gensvm_kernel.h" +#include "gensvm_matrix.h" +#include "gensvm_pred.h" + +#include "util.h" // testing + +void gensvm_predict_labels(struct GenData *data_test, + struct GenData *data_train, struct GenModel *model, + long *predy) +{ + if (model->kerneltype == K_LINEAR) + gensvm_predict_labels_linear(data_test, model, predy); + else + gensvm_predict_labels_kernel(data_test, data_train, model, + predy); +} + +/** + * @brief Predict class labels of data given and output in predy + * + * @details + * The labels are predicted by mapping each instance in data to the + * simplex space using the matrix V in the given model. Next, for each + * instance the nearest simplex vertex is determined using an Euclidean + * norm. The nearest simplex vertex determines the predicted class label, + * which is recorded in predy. + * + * @param[in] data GenData to predict labels for + * @param[in] model GenModel with optimized V + * @param[out] predy pre-allocated vector to record predictions in + */ +void gensvm_predict_labels_linear(struct GenData *data, + struct GenModel *model, long *predy) +{ + long i, j, k, label; + double norm, min_dist; + + long n = data->n; // note that model->n is the size of the training sample. + long m = data->m; + long K = model->K; //data->K does not necessarily equal the original K. + + double *S = Calloc(double, K-1); + double *ZV = Calloc(double, n*(K-1)); + double *U = Calloc(double, K*(K-1)); + + // Get the simplex matrix + gensvm_simplex_gen(K, U); + + // Generate the simplex-space vectors + cblas_dgemm( + CblasRowGenor, + CblasNoTrans, + CblasNoTrans, + n, + K-1, + m+1, + 1.0, + data->Z, + m+1, + model->V, + K-1, + 0.0, + ZV, + K-1); + + // Calculate the distance to each of the vertices of the simplex. + // The closest vertex defines the class label. + for (i=0; in; + long n_test = data_test->n; + long r = model->m; + long K = model->K; + + double *K2 = NULL; + gensvm_make_crosskernel(model, data_train, data_test, &K2); + + double *S = Calloc(double, K-1); + double *ZV = Calloc(double, n_test*(r+1)); + double *KPS = Calloc(double, n_test*(r+1)); + double *U = Calloc(double, K*(K-1)); + + gensvm_simplex_gen(K, U); + + // were doing the computations explicitly since P is included in + // data_train->Z. Might want to look at this some more if it turns out + // to be slow. + + double value, rowvalue; + for (i=0; iZ, r+1, k, + j); + value += rowvalue; + } + value *= matrix_get(data_train->J, 1, j, 0); + matrix_set(KPS, r+1, i, j, value); + } + matrix_set(KPS, r+1, i, 0, 1.0); + } + + cblas_dgemm( + CblasRowGenor, + CblasNoTrans, + CblasNoTrans, + n_test, + K-1, + r+1, + 1.0, + KPS, + r+1, + model->V, + K-1, + 0.0, + ZV, + K-1); + + for (i=0; in; i++) + if (data->y[i] == predy[i]) + correct++; + + performance = ((double) correct)/((double) data->n)* 100.0; + + return performance; +} diff --git a/src/gensvm_sv.c b/src/gensvm_sv.c new file mode 100644 index 0000000..787b869 --- /dev/null +++ b/src/gensvm_sv.c @@ -0,0 +1,45 @@ +/** + * @file gensvm_sv.c + * @author Gertjan van den Burg + * @date May, 2014 + * @brief Calculate the number of support vectors + * + * @details + * The function in this file can be used to calculate the number of support + * vectors are left in a model. + * + */ + +#include "gensvm.h" +#include "gensvm_matrix.h" + +/** + * @brief Calculate the number of support vectors in a model + * + * @details + * If an object is correctly classified, the number of classes for which the + * error q is larger than 1, is K-1 (i.e., there is no error w.r.t. any of the + * other classes). All objects for which this is not the case are thus support + * vectors. + * + * @param[in] model GenModel with solution + * @param[in] data GenData to be used + * @return number of support vectors with this solution + * + */ +long gensvm_num_sv(struct GenModel *model, struct GenData *data) +{ + long i, j, num_correct, num_sv = 0; + double value; + + for (i=0; in; i++) { + num_correct = 0; + for (j=0; jK; j++) { + value = matrix_get(model->Q, data->K, i, j); + num_correct += (value > 1); + } + num_sv += (num_correct < data->K - 1); + } + + return num_sv; +} diff --git a/src/gensvm_train.c b/src/gensvm_train.c new file mode 100644 index 0000000..9deac80 --- /dev/null +++ b/src/gensvm_train.c @@ -0,0 +1,534 @@ +/** + * @file gensvm_train.c + * @author Gertjan van den Burg + * @date August 9, 2013 + * @brief Main functions for training the GenSVM solution. + * + * @details + * Contains update and loss functions used to actually find + * the optimal V. + * + */ + +#include +#include + +#include "libGenSVM.h" +#include "gensvm.h" +#include "gensvm_lapack.h" +#include "gensvm_matrix.h" +#include "gensvm_sv.h" +#include "gensvm_train.h" +#include "util.h" + +/** + * Maximum number of iterations of the algorithm. + */ +#define MAX_ITER 1000000000 + +/** + * @brief The main training loop for GenSVM + * + * @details + * This function is the main training function. This function + * handles the optimization of the model with the given model parameters, with + * the data given. On return the matrix GenModel::V contains the optimal + * weight matrix. + * + * In this function, step doubling is used in the majorization algorithm after + * a burn-in of 50 iterations. If the training is finished, GenModel::t and + * GenModel::W are extracted from GenModel::V. + * + * @param[in,out] model the GenModel to be trained. Contains optimal + * V on exit. + * @param[in] data the GenData to train the model with. + */ +void gensvm_optimize(struct GenModel *model, struct GenData *data) +{ + long i, j, it = 0; + double L, Lbar, value; + + long n = model->n; + long m = model->m; + long K = model->K; + + double *B = Calloc(double, n*(K-1)); + double *ZV = Calloc(double, n*(K-1)); + double *ZAZ = Calloc(double, (m+1)*(m+1)); + double *ZAZV = Calloc(double, (m+1)*(K-1)); + double *ZAZVT = Calloc(double, (m+1)*(K-1)); + + note("Starting main loop.\n"); + note("Dataset:\n"); + note("\tn = %i\n", n); + note("\tm = %i\n", m); + note("\tK = %i\n", K); + note("Parameters:\n"); + note("\tkappa = %f\n", model->kappa); + note("\tp = %f\n", model->p); + note("\tlambda = %15.16f\n", model->lambda); + note("\tepsilon = %g\n", model->epsilon); + note("\n"); + + gensvm_simplex_gen(model->K, model->U); + gensvm_simplex_diff(model, data); + gensvm_category_matrix(model, data); + + L = gensvm_get_loss(model, data, ZV); + Lbar = L + 2.0*model->epsilon*L; + + while ((it < MAX_ITER) && (Lbar - L)/L > model->epsilon) + { + // ensure V contains newest V and Vbar contains V from + // previous + gensvm_get_update(model, data, B, ZAZ, ZAZV, ZAZVT); + if (it > 50) + gensvm_step_doubling(model); + + Lbar = L; + L = gensvm_get_loss(model, data, ZV); + + if (it%100 == 0) + note("iter = %li, L = %15.16f, Lbar = %15.16f, " + "reldiff = %15.16f\n", it, L, Lbar, (Lbar - L)/L); + it++; + } + if (L > Lbar) + fprintf(stderr, "Negative step occurred in majorization.\n"); + + note("optimization finished, iter = %li, loss = %15.16f, " + "rel. diff. = %15.16f\n", it-1, L, + (Lbar - L)/L); + note("number of support vectors: %li\n", gensvm_num_sv(model, data)); + + model->training_error = (Lbar - L)/L; + + for (i=0; it[i] = matrix_get(model->V, K-1, 0, i); + for (i=1; iV, K-1, i, j); + matrix_set(model->W, K-1, i-1, j, value); + } + } + free(B); + free(ZV); + free(ZAZ); + free(ZAZV); + free(ZAZVT); +} + +/** + * @brief Calculate the current value of the loss function + * + * @details + * The current loss function value is calculated based on the matrix V in the + * given model. Note that the matrix ZV is passed explicitly to avoid having + * to reallocate memory at every step. + * + * @param[in] model GenModel structure which holds the current + * estimate V + * @param[in] data GenData structure + * @param[in,out] ZV pre-allocated matrix ZV which is updated on + * output + * @returns the current value of the loss function + */ +double gensvm_get_loss(struct GenModel *model, struct GenData *data, + double *ZV) +{ + long i, j; + long n = data->n; + long K = data->K; + long m = data->m; + + double value, rowvalue, loss = 0.0; + + gensvm_calculate_errors(model, data, ZV); + gensvm_calculate_huber(model); + + for (i=0; iH, K, i, j); + value = pow(value, model->p); + value *= matrix_get(model->R, K, i, j); + rowvalue += value; + } + rowvalue = pow(rowvalue, 1.0/(model->p)); + rowvalue *= model->rho[i]; + loss += rowvalue; + } + loss /= ((double) n); + + value = 0; + for (i=0; iV, K-1, i, j), 2.0); + } + value += data->J[i] * rowvalue; + } + loss += model->lambda * value; + + return loss; +} + +/** + * @brief Perform a single step of the majorization algorithm to update V + * + * @details + * This function contains the main update calculations of the algorithm. These + * calculations are necessary to find a new update V. The calculations exist of + * recalculating the majorization coefficients for all instances and all + * classes, and solving a linear system to find V. + * + * Because the function gensvm_get_update() is always called after a call to + * gensvm_get_loss() with the same GenModel::V, it is unnecessary to calculate + * the updated errors GenModel::Q and GenModel::H here too. This saves on + * computation time. + * + * In calculating the majorization coefficients we calculate the elements of a + * diagonal matrix A with elements + * @f[ + * A_{i, i} = \frac{1}{n} \rho_i \sum_{j \neq k} \left[ + * \varepsilon_i a_{ijk}^{(p)} + (1 - \varepsilon_i) \omega_i + * a_{ijk}^{(p)} \right], + * @f] + * where @f$ k = y_i @f$. + * Since this matrix is only used to calculate the matrix @f$ Z' A Z @f$, it is + * efficient to update a matrix ZAZ through consecutive rank 1 updates with + * a single element of A and the corresponding row of Z. The BLAS function + * dsyr is used for this. + * + * The B matrix is has rows + * @f[ + * \boldsymbol{\beta}_i' = \frac{1}{n} \rho_i \sum_{j \neq k} \left[ + * \varepsilon_i \left( b_{ijk}^{(1)} - a_{ijk}^{(1)} + * \overline{q}_i^{(kj)} \right) + (1 - \varepsilon_i) + * \omega_i \left( b_{ijk}^{(p)} - a_{ijk}^{(p)} + * \overline{q}_i^{(kj)} \right) \right] + * \boldsymbol{\delta}_{kj}' + * @f] + * This is also split into two cases, one for which @f$ \varepsilon_i = 1 @f$, + * and one for when it is 0. The 3D simplex difference matrix is used here, in + * the form of the @f$ \boldsymbol{\delta}_{kj}' @f$. + * + * Finally, the following system is solved + * @f[ + * (\textbf{Z}'\textbf{AZ} + \lambda \textbf{J})\textbf{V} = + * (\textbf{Z}'\textbf{AZ}\overline{\textbf{V}} + \textbf{Z}' + * \textbf{B}) + * @f] + * solving this system is done through dposv(). + * + * @todo + * Consider allocating IPIV and WORK at a higher level, they probably don't + * change much during the iterations. + * + * @param [in,out] model model to be updated + * @param [in] data data used in model + * @param [in] B pre-allocated matrix used for linear coefficients + * @param [in] ZAZ pre-allocated matrix used in system + * @param [in] ZAZV pre-allocated matrix used in system solving + * @param [in] ZAZVT pre-allocated matrix used in system solving + */ +void gensvm_get_update(struct GenModel *model, struct GenData *data, double *B, + double *ZAZ, double *ZAZV, double *ZAZVT) +{ + int status, class; + long i, j, k; + double Avalue, Bvalue; + double omega, value, a, b, q, h, r; + + long n = model->n; + long m = model->m; + long K = model->K; + + double kappa = model->kappa; + double p = model->p; + double *rho = model->rho; + + // constants which are used often throughout + const double a2g2 = 0.25*p*(2.0*p - 1.0)*pow((kappa+1.0)/2.0,p-2.0); + const double in = 1.0/((double) n); + + // clear matrices + Memset(B, double, n*(K-1)); + Memset(ZAZ, double, (m+1)*(m+1)); + + b = 0; + for (i=0; iH, K, i, j); + r = matrix_get(model->R, K, i, j); + value += (h*r > 0) ? 1 : 0; + omega += pow(h, p)*r; + } + class = (value <= 1.0) ? 1 : 0; + omega = (1.0/p)*pow(omega, 1.0/p - 1.0); + + Avalue = 0; + if (class == 1) { + for (j=0; jQ, K, i, j); + if (q <= -kappa) { + a = 0.25/(0.5 - kappa/2.0 - q); + b = 0.5; + } else if (q <= 1.0) { + a = 1.0/(2.0*kappa + 2.0); + b = (1.0 - q)*a; + } else { + a = -0.25/(0.5 - kappa/2.0 - q); + b = 0; + } + for (k=0; kUU, K-1, K, i, k, j); + matrix_add(B, K-1, i, k, Bvalue); + } + Avalue += a*matrix_get(model->R, K, i, j); + } + } else { + if (2.0 - p < 0.0001) { + for (j=0; jQ, K, i, j); + if (q <= -kappa) { + b = 0.5 - kappa/2.0 - q; + } else if ( q <= 1.0) { + b = pow(1.0 - q, 3.0)/( + 2.0*pow(kappa + 1.0, + 2.0)); + } else { + b = 0; + } + for (k=0; kUU, + K-1, + K, + i, + k, + j); + matrix_add( + B, + K-1, + i, + k, + Bvalue); + } + } + Avalue = 1.5*(K - 1.0); + } else { + for (j=0; jQ, K, i, j); + if (q <= (p + kappa - 1.0)/(p - 2.0)) { + a = 0.25*pow(p, 2.0)*pow( + 0.5 - kappa/2.0 - q, + p - 2.0); + } else if (q <= 1.0) { + a = a2g2; + } else { + a = 0.25*pow(p, 2.0)*pow( + (p/(p - 2.0))* + (0.5 - kappa/2.0 - q), + p - 2.0); + b = a*(2.0*q + kappa - 1.0)/ + (p - 2.0) + + 0.5*p*pow( + p/(p - 2.0)* + (0.5 - kappa/ + 2.0 - q), + p - 1.0); + } + if (q <= -kappa) { + b = 0.5*p*pow( + 0.5 - kappa/2.0 - q, + p - 1.0); + } else if ( q <= 1.0) { + b = p*pow(1.0 - q, + 2.0*p - 1.0)/ + pow(2*kappa+2.0, p); + } + for (k=0; kUU, + K-1, + K, + i, + k, + j); + matrix_add( + B, + K-1, + i, + k, + Bvalue); + } + Avalue += a*matrix_get(model->R, + K, i, j); + } + } + Avalue *= omega; + } + Avalue *= in * rho[i]; + + // Now we calculate the matrix ZAZ. Since this is + // guaranteed to be symmetric, we only calculate the + // upper part of the matrix, and then copy this over + // to the lower part after all calculations are done. + // Note that the use of dsym is faster than dspr, even + // though dspr uses less memory. + cblas_dsyr( + CblasRowGenor, + CblasUpper, + m+1, + Avalue, + &data->Z[i*(m+1)], + 1, + ZAZ, + m+1); + } + // Copy upper to lower (necessary because we need to switch + // to Col-Genor order for LAPACK). + /* + for (i=0; iV, + K-1, + 0.0, + ZAZV, + K-1); + + cblas_dgemm( + CblasRowGenor, + CblasTrans, + CblasNoTrans, + m+1, + K-1, + n, + 1.0, + data->Z, + m+1, + B, + K-1, + 1.0, + ZAZV, + K-1); + + /* + * Add lambda to all diagonal elements except the first one. Recall + * that ZAZ is of size m+1 and is symmetric. + */ + i = 0; + for (j=0; jlambda * data->J[j+1]; + } + + // For the LAPACK call we need to switch to Column- + // Genor order. This is unnecessary for the matrix + // ZAZ because it is symmetric. The matrix ZAZV + // must be converted however. + for (i=0; iVbar; + model->Vbar = model->V; + model->V = ZAZVT; + ZAZVT = ptr; + */ + + for (i=0; iV, K-1, i, j); + matrix_set(model->Vbar, K-1, i, j, value); + value = matrix_get(ZAZV, K-1, i, j); + matrix_set(model->V, K-1, i, j, value); + } + } +} diff --git a/src/gensvm_train_dataset.c b/src/gensvm_train_dataset.c new file mode 100644 index 0000000..3034bb4 --- /dev/null +++ b/src/gensvm_train_dataset.c @@ -0,0 +1,748 @@ +/** + * @file gensvm_train_dataset.c + * @author Gertjan van den Burg + * @date January, 2014 + * @brief Functions for finding the optimal parameters for the dataset + * + * @details + * The GenSVM algorithm takes a number of parameters. The functions in + * this file are used to find the optimal parameters. + */ + +#include +#include + +#include "crossval.h" +#include "libGenSVM.h" +#include "gensvm.h" +#include "gensvm_init.h" +#include "gensvm_kernel.h" +#include "gensvm_matrix.h" +#include "gensvm_train.h" +#include "gensvm_train_dataset.h" +#include "gensvm_pred.h" +#include "util.h" +#include "timer.h" + +extern FILE *GENSVM_OUTPUT_FILE; + +/** + * @brief Initialize a Queue from a Training instance + * + * @details + * A Training instance describes the grid to search over. This funtion + * creates all tasks that need to be performed and adds these to + * a Queue. Each task contains a pointer to the train and test datasets + * which are supplied. Note that the tasks are created in a specific order of + * the parameters, to ensure that the GenModel::V of a previous parameter + * set provides the best possible initial estimate of GenModel::V for the next + * parameter set. + * + * @param[in] training Training struct describing the grid search + * @param[in] queue pointer to a Queue that will be used to + * add the tasks to + * @param[in] train_data GenData of the training set + * @param[in] test_data GenData of the test set + * + */ +void make_queue(struct Training *training, struct Queue *queue, + struct GenData *train_data, struct GenData *test_data) +{ + long i, j, k; + long N, cnt = 0; + struct Task *task; + queue->i = 0; + + N = training->Np; + N *= training->Nl; + N *= training->Nk; + N *= training->Ne; + N *= training->Nw; + // these parameters are not necessarily non-zero + N *= training->Ng > 0 ? training->Ng : 1; + N *= training->Nc > 0 ? training->Nc : 1; + N *= training->Nd > 0 ? training->Nd : 1; + + queue->tasks = Malloc(struct Task *, N); + queue->N = N; + + // initialize all tasks + for (i=0; iID = i; + task->train_data = train_data; + task->test_data = test_data; + task->folds = training->folds; + task->kerneltype = training->kerneltype; + task->kernelparam = Calloc(double, training->Ng + + training->Nc + training->Nd); + queue->tasks[i] = task; + } + + // These loops mimick a large nested for loop. The advantage is that + // Nd, Nc and Ng which are on the outside of the nested for loop can + // now be zero, without large modification (see below). Whether this + // is indeed better than the nested for loop has not been tested. + cnt = 1; + i = 0; + while (i < N ) + for (j=0; jNp; j++) + for (k=0; ktasks[i]->p = training->ps[j]; + i++; + } + + cnt *= training->Np; + i = 0; + while (i < N ) + for (j=0; jNl; j++) + for (k=0; ktasks[i]->lambda = + training->lambdas[j]; + i++; + } + + cnt *= training->Nl; + i = 0; + while (i < N ) + for (j=0; jNk; j++) + for (k=0; ktasks[i]->kappa = training->kappas[j]; + i++; + } + + cnt *= training->Nk; + i = 0; + while (i < N ) + for (j=0; jNw; j++) + for (k=0; ktasks[i]->weight_idx = + training->weight_idxs[j]; + i++; + } + + cnt *= training->Nw; + i = 0; + while (i < N ) + for (j=0; jNe; j++) + for (k=0; ktasks[i]->epsilon = + training->epsilons[j]; + i++; + } + + cnt *= training->Ne; + i = 0; + while (i < N && training->Ng > 0) + for (j=0; jNg; j++) + for (k=0; ktasks[i]->kernelparam[0] = + training->gammas[j]; + i++; + } + + cnt *= training->Ng > 0 ? training->Ng : 1; + i = 0; + while (i < N && training->Nc > 0) + for (j=0; jNc; j++) + for (k=0; ktasks[i]->kernelparam[1] = + training->coefs[j]; + i++; + } + + cnt *= training->Nc > 0 ? training->Nc : 1; + i = 0; + while (i < N && training->Nd > 0) + for (j=0; jNd; j++) + for (k=0; ktasks[i]->kernelparam[2] = + training->degrees[j]; + i++; + } +} + +/** + * @brief Get new Task from Queue + * + * @details + * Return a pointer to the next Task in the Queue. If no Task instances are + * left, NULL is returned. The internal counter Queue::i is used for finding + * the next Task. + * + * @param[in] q Queue instance + * @returns pointer to next Task + * + */ +struct Task *get_next_task(struct Queue *q) +{ + long i = q->i; + if (i < q->N) { + q->i++; + return q->tasks[i]; + } + return NULL; +} + +/** + * @brief Comparison function for Tasks based on performance + * + * @details + * To be able to sort Task structures on the performance of their specific + * set of parameters, this comparison function is implemented. Task structs + * are sorted with highest performance first. + * + * @param[in] elem1 Task 1 + * @param[in] elem2 Task 2 + * @returns result of inequality of Task 1 performance over + * Task 2 performance + */ +int tasksort(const void *elem1, const void *elem2) +{ + const struct Task *t1 = (*(struct Task **) elem1); + const struct Task *t2 = (*(struct Task **) elem2); + return (t1->performance > t2->performance); +} + +/** + * @brief Comparison function for doubles + * + * @details + * Similar to tasksort() only now for two doubles. + * + * @param[in] elem1 number 1 + * @param[in] elem2 number 2 + * @returns comparison of number 1 larger than number 2 + */ +int doublesort(const void *elem1, const void *elem2) +{ + const double t1 = (*(double *) elem1); + const double t2 = (*(double *) elem2); + return t1 > t2; +} + +/** + * @brief Calculate the percentile of an array of doubles + * + * @details + * The percentile of performance is used to find the top performing + * configurations. Since no standard definition of the percentile exists, we + * use the method used in MATLAB and Octave. Since calculating the percentile + * requires a sorted list of the values, a local copy is made first. + * + * @param[in] values array of doubles + * @param[in] N length of the array + * @param[in] p percentile to calculate ( 0 <= p <= 1.0 ). + * @returns the p-th percentile of the values + */ +double prctile(double *values, long N, double p) +{ + long i; + double pi, pr, boundary; + double *local = Malloc(double, N); + for (i=0; itasks, q->N, sizeof(struct Task *), tasksort); + p = 0.95*q->N + 0.5; + pi = maximum(minimum(floor(p), q->N-1), 1); + pr = maximum(minimum(p - pi, 1), 0); + boundary = (1 - pr)*q->tasks[((long) pi)-1]->performance; + boundary += pr*q->tasks[((long) pi)]->performance; + note("boundary determined at: %f\n", boundary); + + // find the number of tasks that perform at least as good as the 95th + // percentile + N = 0; + for (i=0; iN; i++) + if (q->tasks[i]->performance >= boundary) + N++; + note("Number of items: %li\n", N); + std = Calloc(double, N); + mean = Calloc(double, N); + time = Calloc(double, N); + perf = Calloc(double, N*repeats); + + // create a new task queue with the tasks which perform well + nq->tasks = Malloc(struct Task *, N); + for (i=q->N-1; i>q->N-N-1; i--) + nq->tasks[q->N-i-1] = q->tasks[i]; + nq->N = N; + nq->i = 0; + + // for each task run the consistency repeats + for (i=0; in = 0; + model->m = task->train_data->m; + model->K = task->train_data->K; + gensvm_allocate_model(model); + gensvm_seed_model_V(NULL, model, task->train_data); + } + + time[i] = 0.0; + note("(%02li/%02li:%03li)\t", i+1, N, task->ID); + for (r=0; rtrain_data, + task->folds); + loop_e = clock(); + time[i] += elapsed_time(loop_s, loop_e); + matrix_set(perf, repeats, i, r, p); + mean[i] += p/((double) repeats); + } else { + note("Only cv is implemented\n"); + exit(1); + } + note("%3.3f\t", p); + // this is done because if we reuse the V it's not a + // consistency check + gensvm_seed_model_V(NULL, model, task->train_data); + } + for (r=0; r 1) { + std[i] /= ((double) repeats) - 1.0; + std[i] = sqrt(std[i]); + } else + std[i] = 0.0; + note("(m = %3.3f, s = %3.3f, t = %3.3f)\n", + mean[i], std[i], time[i]); + } + + // find the best overall configurations: those with high average + // performance and low deviation in the performance + note("\nBest overall configuration(s):\n"); + note("ID\tweights\tepsilon\t\tp\t\tkappa\t\tlambda\t\t" + "mean_perf\tstd_perf\ttime_perf\n"); + p = 0.0; + bool breakout = false; + while (breakout == false) { + pi = prctile(mean, N, (100.0-p)/100.0); + pr = prctile(std, N, p/100.0); + pt = prctile(time, N, p/100.0); + for (i=0; itasks[i]->ID, + nq->tasks[i]->weight_idx, + nq->tasks[i]->epsilon, + nq->tasks[i]->p, + nq->tasks[i]->kappa, + nq->tasks[i]->lambda, + mean[i], + std[i], + time[i]); + breakout = true; + } + p += 1.0; + } + + free(nq->tasks); + free(nq); + free(model); + free(perf); + free(std); + free(mean); + free(time); +} + +/** + * @brief Run cross validation with a seed model + * + * @details + * This is an implementation of cross validation which uses the optimal + * parameters GenModel::V of a previous fold as initial conditions for + * GenModel::V of the next fold. An initial seed for V can be given through the + * seed_model parameter. If seed_model is NULL, random starting values are + * used. + * + * @param[in] model GenModel with the configuration to train + * @param[in] seed_model GenModel with a seed for GenModel::V + * @param[in] data GenData with the dataset + * @param[in] folds number of cross validation folds + * @returns performance (hitrate) of the configuration on + * cross validation + */ +double cross_validation(struct GenModel *model, struct GenData *data, + long folds) +{ + FILE *fid; + + long f, *predy; + double performance, total_perf = 0; + struct GenData *train_data, *test_data; + + long *cv_idx = Calloc(long, data->n); + + train_data = gensvm_init_data(); + test_data = gensvm_init_data(); + + // create splits + gensvm_make_cv_split(data->n, folds, cv_idx); + + for (f=0; fn, train_data->m); + + gensvm_initialize_weights(train_data, model); + + // train the model (without output) + fid = GENSVM_OUTPUT_FILE; + GENSVM_OUTPUT_FILE = NULL; + gensvm_optimize(model, train_data); + GENSVM_OUTPUT_FILE = fid; + + // calculate prediction performance on test set + predy = Calloc(long, test_data->n); + gensvm_predict_labels(test_data, train_data, model, predy); + performance = gensvm_prediction_perf(test_data, predy); + total_perf += performance * test_data->n; + + free(predy); + free(train_data->y); + free(train_data->Z); + free(test_data->y); + free(test_data->Z); + } + + free(train_data); + free(test_data); + + total_perf /= ((double) data->n); + + return total_perf; +} + +/** + * @brief Run the grid search for a cross validation dataset + * + * @details + * Given a Queue of Task struct to be trained, a grid search is launched to + * find the optimal parameter configuration. As is also done within + * cross_validation(), the optimal weights of one parameter set are used as + * initial estimates for GenModel::V in the next parameter set. Note that to + * optimally exploit this feature of the optimization algorithm, the order in + * which tasks are considered is important. This is considered in + * make_queue(). + * + * The performance found by cross validation is stored in the Task struct. + * + * @param[in,out] q Queue with Task instances to run + */ +void start_training_cv(struct Queue *q) +{ + double perf, current_max = 0; + struct Task *task = get_next_task(q); + struct GenModel *model = gensvm_init_model(); + clock_t main_s, main_e, loop_s, loop_e; + + model->n = 0; + model->m = task->train_data->m; + model->K = task->train_data->K; + gensvm_allocate_model(model); + gensvm_seed_model_V(NULL, model, task->train_data); + + main_s = clock(); + while (task) { + print_progress_string(task, q->N); + make_model_from_task(task, model); + + loop_s = clock(); + perf = cross_validation(model, task->train_data, task->folds); + loop_e = clock(); + current_max = maximum(current_max, perf); + + note("\t%3.3f%% (%3.3fs)\t(best = %3.3f%%)\n", perf, + elapsed_time(loop_s, loop_e), + current_max); + + q->tasks[task->ID]->performance = perf; + task = get_next_task(q); + } + main_e = clock(); + + note("\nTotal elapsed time: %8.8f seconds\n", + elapsed_time(main_s, main_e)); + + gensvm_free_model(model); +} + +/** + * @brief Run the grid search for a train/test dataset + * + * @details + * This function is similar to start_training_cv(), except that the + * pre-determined training set is used only once, and the pre-determined test + * set is used for validation. + * + * @todo + * It would probably be better to train the model on the training set using + * cross validation and only use the test set when comparing with other + * methods. The way it is now, you're finding out which parameters predict + * _this_ test set best, which is not what you want. This function should + * therefore not be used and is considered deprecated, to be removed in the + * future . + * + * @param[in] q Queue with Task structs to run + * + */ +void start_training_tt(struct Queue *q) +{ + FILE *fid; + + long c = 0; + long *predy; + double total_perf, current_max = 0; + + struct Task *task = get_next_task(q); + struct GenModel *seed_model = gensvm_init_model(); + + clock_t main_s, main_e; + clock_t loop_s, loop_e; + + seed_model->m = task->train_data->m; + seed_model->K = task->train_data->K; + gensvm_allocate_model(seed_model); + gensvm_seed_model_V(NULL, seed_model, task->train_data); + + main_s = clock(); + while (task) { + total_perf = 0; + note("(%li/%li)\tw = %li\te = %f\tp = %f\tk = %f\tl = %f\t", + c+1, q->N, task->weight_idx, task->epsilon, + task->p, task->kappa, task->lambda); + loop_s = clock(); + struct GenModel *model = gensvm_init_model(); + make_model_from_task(task, model); + + model->n = task->train_data->n; + model->m = task->train_data->m; + model->K = task->train_data->K; + + gensvm_allocate_model(model); + gensvm_initialize_weights(task->train_data, model); + gensvm_seed_model_V(seed_model, model, task->train_data); + + fid = GENSVM_OUTPUT_FILE; + GENSVM_OUTPUT_FILE = NULL; + gensvm_optimize(model, task->train_data); + GENSVM_OUTPUT_FILE = fid; + + predy = Calloc(long, task->test_data->n); + gensvm_predict_labels(task->test_data, task->train_data, + model, predy); + if (task->test_data->y != NULL) + total_perf = gensvm_prediction_perf(task->test_data, + predy); + gensvm_seed_model_V(model, seed_model, task->train_data); + + gensvm_free_model(model); + free(predy); + note("."); + loop_e = clock(); + current_max = maximum(current_max, total_perf); + note("\t%3.3f%% (%3.3fs)\t(best = %3.3f%%)\n", total_perf, + elapsed_time(loop_s, loop_e), current_max); + q->tasks[task->ID]->performance = total_perf; + task = get_next_task(q); + } + main_e = clock(); + + note("\nTotal elapsed time: %8.8f seconds\n", + elapsed_time(main_s, main_e)); + free(task); + gensvm_free_model(seed_model); +} + +/** + * @brief Free the Queue struct + * + * @details + * Freeing the allocated memory of the Queue means freeing every Task struct + * and then freeing the Queue. + * + * @param[in] q Queue to be freed + * + */ +void free_queue(struct Queue *q) +{ + long i; + for (i=0; iN; i++) { + free(q->tasks[i]->kernelparam); + free(q->tasks[i]); + } + free(q->tasks); + free(q); +} + +/** + * @brief Copy parameters from Task to GenModel + * + * @details + * A Task struct only contains the parameters of the GenModel to be estimated. + * This function is used to copy these parameters. + * + * @param[in] task Task instance with parameters + * @param[in,out] model GenModel to which the parameters are copied + */ +void make_model_from_task(struct Task *task, struct GenModel *model) +{ + // copy basic model parameters + model->weight_idx = task->weight_idx; + model->epsilon = task->epsilon; + model->p = task->p; + model->kappa = task->kappa; + model->lambda = task->lambda; + + // copy kernel parameters + model->kerneltype = task->kerneltype; + model->kernelparam = task->kernelparam; +} + +/** + * @brief Copy model parameters between two GenModel structs + * + * @details + * The parameters copied are GenModel::weight_idx, GenModel::epsilon, + * GenModel::p, GenModel::kappa, and GenModel::lambda. + * + * @param[in] from GenModel to copy parameters from + * @param[in,out] to GenModel to copy parameters to + */ +void copy_model(struct GenModel *from, struct GenModel *to) +{ + to->weight_idx = from->weight_idx; + to->epsilon = from->epsilon; + to->p = from->p; + to->kappa = from->kappa; + to->lambda = from->lambda; + + to->kerneltype = from->kerneltype; + switch (to->kerneltype) { + case K_LINEAR: + break; + case K_POLY: + to->kernelparam = Malloc(double, 3); + to->kernelparam[0] = from->kernelparam[0]; + to->kernelparam[1] = from->kernelparam[1]; + to->kernelparam[2] = from->kernelparam[2]; + break; + case K_RBF: + to->kernelparam = Malloc(double, 1); + to->kernelparam[0] = from->kernelparam[0]; + break; + case K_SIGMOID: + to->kernelparam = Malloc(double, 2); + to->kernelparam[0] = from->kernelparam[0]; + to->kernelparam[1] = from->kernelparam[1]; + break; + } +} + +/** + * @brief Print the description of the current task on screen + * + * @details + * To track the progress of the grid search the parameters of the current task + * are written to the output specified in GENSVM_OUTPUT_FILE. Since the + * parameters differ with the specified kernel, this function writes a + * parameter string depending on which kernel is used. + * + * @param[in] task the Task specified + * @param[in] N total number of tasks + * + */ +void print_progress_string(struct Task *task, long N) +{ + char buffer[MAX_LINE_LENGTH]; + sprintf(buffer, "(%03li/%03li)\t", task->ID+1, N); + if (task->kerneltype == K_POLY) + sprintf(buffer + strlen(buffer), "d = %2.2f\t", + task->kernelparam[2]); + if (task->kerneltype == K_POLY || task->kerneltype == K_SIGMOID) + sprintf(buffer + strlen(buffer), "c = %2.2f\t", + task->kernelparam[1]); + if (task->kerneltype == K_POLY || task->kerneltype == K_SIGMOID || + task->kerneltype == K_RBF) + sprintf(buffer + strlen(buffer), "g = %3.3f\t", + task->kernelparam[0]); + sprintf(buffer + strlen(buffer), "eps = %g\tw = %i\tk = %2.2f\t" + "l = %f\tp = %2.2f\t", task->epsilon, + task->weight_idx, task->kappa, task->lambda, task->p); + note(buffer); +} diff --git a/src/libGenSVM.c b/src/libGenSVM.c new file mode 100644 index 0000000..bb48673 --- /dev/null +++ b/src/libGenSVM.c @@ -0,0 +1,326 @@ +/** + * @file libGenSVM.c + * @author Gertjan van den Burg + * @date August 8, 2013 + * @brief Main functions for the GenSVM algorithm + * + * @details + * The functions in this file are all functions needed + * to calculate the optimal separation boundaries for + * a multiclass classification problem, using the + * GenSVM algorithm. + * + */ + +#include +#include + +#include "libGenSVM.h" +#include "gensvm.h" +#include "gensvm_matrix.h" + +inline double rnd() { return (double) rand()/0x7FFFFFFF; } + +/** + * @brief Generate matrix of simplex vertex coordinates + * + * @details + * Generate the simplex matrix. Each row of the created + * matrix contains the coordinate vector of a single + * vertex of the K-simplex in K-1 dimensions. The simplex + * generated is a special simplex with edges of length 1. + * The simplex matrix U must already have been allocated. + * + * @param[in] K number of classes + * @param[in,out] U simplex matrix of size K * (K-1) + */ +void gensvm_simplex_gen(long K, double *U) +{ + long i, j; + for (i=0; in; + long K = model->K; + + for (i=0; iy[i] != j+1) + matrix_set(model->R, K, i, j, 1.0); + else + matrix_set(model->R, K, i, j, 0.0); + } + } +} + +/** + * @brief Generate the simplex difference matrix + * + * @details + * The simplex difference matrix is a 3D matrix which is constructed + * as follows. For each instance i, the difference vectors between the row of + * the simplex matrix corresponding to the class label of instance i and the + * other rows of the simplex matrix are calculated. These difference vectors + * are stored in a matrix, which is one horizontal slice of the 3D matrix. + * + * @param[in,out] model the corresponding GenModel + * @param[in] data the corresponding GenData + * + */ +void gensvm_simplex_diff(struct GenModel *model, struct GenData *data) +{ + long i, j, k; + double value; + + long n = model->n; + long K = model->K; + + for (i=0; iU, K-1, data->y[i]-1, j); + value -= matrix_get(model->U, K-1, k, j); + matrix3_set(model->UU, K-1, K, i, j, k, value); + } + } + } +} + +/** + * @brief Calculate the scalar errors + * + * @details + * Calculate the scalar errors q based on the current estimate of V, and + * store these in Q. It is assumed that the memory for Q has already been + * allocated. In addition, the matrix ZV is calculated here. It is assigned + * to a pre-allocated block of memory, which is passed to this function. + * + * @param[in,out] model the corresponding GenModel + * @param[in] data the corresponding GenData + * @param[in,out] ZV a pointer to a memory block for ZV. On exit + * this block is updated with the new ZV matrix + * calculated with GenModel::V. + * + */ +void gensvm_calculate_errors(struct GenModel *model, struct GenData *data, + double *ZV) +{ + long i, j, k; + double a, value; + + long n = model->n; + long m = model->m; + long K = model->K; + + cblas_dgemm( + CblasRowGenor, + CblasNoTrans, + CblasNoTrans, + n, + K-1, + m+1, + 1.0, + data->Z, + m+1, + model->V, + K-1, + 0.0, + ZV, + K-1); + + Memset(model->Q, double, n*K); + for (i=0; iUU, K-1, K, i, + j, k); + matrix_add(model->Q, K, i, k, value); + } + } + } +} + +/** + * @brief Calculate the Huber hinge errors + * + * @details + * For each of the scalar errors in Q the Huber hinge errors are + * calculated. The Huber hinge is here defined as + * @f[ + * h(q) = + * \begin{dcases} + * 1 - q - \frac{\kappa + 1}{2} & \text{if } q \leq -\kappa \\ + * \frac{1}{2(\kappa + 1)} ( 1 - q)^2 & \text{if } q \in (-\kappa, 1] \\ + * 0 & \text{if } q > 1 + * \end{dcases} + * @f] + * + * @param[in,out] model the corresponding GenModel + */ +void gensvm_calculate_huber(struct GenModel *model) +{ + long i, j; + double q, value; + + for (i=0; in; i++) { + for (j=0; jK; j++) { + q = matrix_get(model->Q, model->K, i, j); + value = 0.0; + if (q <= -model->kappa) { + value = 1.0 - q - (model->kappa+1.0)/2.0; + } else if (q <= 1.0) { + value = 1.0/(2.0*model->kappa+2.0)*pow(1.0 - q, 2.0); + } + matrix_set(model->H, model->K, i, j, value); + } + } +} + +/** + * @brief seed the matrix V from an existing model or using rand + * + * @details + * The matrix V must be seeded before the main_loop() can start. + * This can be done by either seeding it with random numbers or + * using the solution from a previous model on the same dataset + * as initial seed. The latter option usually allows for a + * significant improvement in the number of iterations necessary + * because the seeded model V is closer to the optimal V. + * + * @param[in] from_model GenModel from which to copy V + * @param[in,out] to_model GenModel to which V will be copied + */ +void gensvm_seed_model_V(struct GenModel *from_model, + struct GenModel *to_model, struct GenData *data) +{ + long i, j, k; + double cmin, cmax, value; + + long n = data->n; + long m = data->m; + long K = data->K; + + if (from_model == NULL) { + for (i=0; iZ, m+1, k, i); + cmin = minimum(cmin, value); + cmax = maximum(cmax, value); + } + for (j=0; jV, K-1, i, j, value); + } + } + } else { + for (i=0; iV, K-1, i, j); + matrix_set(to_model->V, K-1, i, j, value); + } + } +} + +/** + * @brief Use step doubling + * + * @details + * Step doubling can be used to speed up the Genorization algorithm. Instead + * of using the value at the minimimum of the majorization function, the value + * ``opposite'' the majorization point is used. This can essentially cut the + * number of iterations necessary to reach the minimum in half. + * + * @param[in] model GenModel containing the augmented parameters + */ +void gensvm_step_doubling(struct GenModel *model) +{ + long i, j; + double value; + + long m = model->m; + long K = model->K; + + for (i=0; iV, K-1, i, j, 2.0); + value = - matrix_get(model->Vbar, K-1, i, j); + matrix_add(model->V, K-1, i, j, value); + } + } +} + +/** + * @brief Initialize instance weights + * + * @details + * Instance weights can for instance be used to add additional weights to + * instances of certain classes. Two default weight possibilities are + * implemented here. The first is unit weights, where each instance gets + * weight 1. + * + * The second are group size correction weights, which are calculated as + * @f[ + * \rho_i = \frac{n}{Kn_k} , + * @f] + * where @f$ n_k @f$ is the number of instances in group @f$ k @f$ and + * @f$ y_i = k @f$. + * + * @param[in] data GenData with the dataset + * @param[in,out] model GenModel with the weight specification. On + * exit GenModel::rho contains the instance + * weights. + */ +void gensvm_initialize_weights(struct GenData *data, struct GenModel *model) +{ + long *groups; + long i; + + long n = model->n; + long K = model->K; + + if (model->weight_idx == 1) { + for (i=0; irho[i] = 1.0; + } + else if (model->weight_idx == 2) { + groups = Calloc(long, K); + for (i=0; iy[i]-1]++; + for (i=0; irho[i] = ((double) n)/((double) (groups[data->y[i]-1]*K)); + } else { + fprintf(stderr, "Unknown weight specification.\n"); + exit(1); + } +} + diff --git a/src/libMSVMMaj.c b/src/libMSVMMaj.c deleted file mode 100644 index df422c0..0000000 --- a/src/libMSVMMaj.c +++ /dev/null @@ -1,326 +0,0 @@ -/** - * @file libMSVMMaj.c - * @author Gertjan van den Burg - * @date August 8, 2013 - * @brief Main functions for the MSVMMaj algorithm - * - * @details - * The functions in this file are all functions needed - * to calculate the optimal separation boundaries for - * a multiclass classification problem, using the - * MSVMMaj algorithm. - * - */ - -#include -#include - -#include "libMSVMMaj.h" -#include "msvmmaj.h" -#include "msvmmaj_matrix.h" - -inline double rnd() { return (double) rand()/0x7FFFFFFF; } - -/** - * @brief Generate matrix of simplex vertex coordinates - * - * @details - * Generate the simplex matrix. Each row of the created - * matrix contains the coordinate vector of a single - * vertex of the K-simplex in K-1 dimensions. The simplex - * generated is a special simplex with edges of length 1. - * The simplex matrix U must already have been allocated. - * - * @param[in] K number of classes - * @param[in,out] U simplex matrix of size K * (K-1) - */ -void msvmmaj_simplex_gen(long K, double *U) -{ - long i, j; - for (i=0; in; - long K = model->K; - - for (i=0; iy[i] != j+1) - matrix_set(model->R, K, i, j, 1.0); - else - matrix_set(model->R, K, i, j, 0.0); - } - } -} - -/** - * @brief Generate the simplex difference matrix - * - * @details - * The simplex difference matrix is a 3D matrix which is constructed - * as follows. For each instance i, the difference vectors between the row of - * the simplex matrix corresponding to the class label of instance i and the - * other rows of the simplex matrix are calculated. These difference vectors - * are stored in a matrix, which is one horizontal slice of the 3D matrix. - * - * @param[in,out] model the corresponding MajModel - * @param[in] data the corresponding MajData - * - */ -void msvmmaj_simplex_diff(struct MajModel *model, struct MajData *data) -{ - long i, j, k; - double value; - - long n = model->n; - long K = model->K; - - for (i=0; iU, K-1, data->y[i]-1, j); - value -= matrix_get(model->U, K-1, k, j); - matrix3_set(model->UU, K-1, K, i, j, k, value); - } - } - } -} - -/** - * @brief Calculate the scalar errors - * - * @details - * Calculate the scalar errors q based on the current estimate of V, and - * store these in Q. It is assumed that the memory for Q has already been - * allocated. In addition, the matrix ZV is calculated here. It is assigned - * to a pre-allocated block of memory, which is passed to this function. - * - * @param[in,out] model the corresponding MajModel - * @param[in] data the corresponding MajData - * @param[in,out] ZV a pointer to a memory block for ZV. On exit - * this block is updated with the new ZV matrix - * calculated with MajModel::V. - * - */ -void msvmmaj_calculate_errors(struct MajModel *model, struct MajData *data, - double *ZV) -{ - long i, j, k; - double a, value; - - long n = model->n; - long m = model->m; - long K = model->K; - - cblas_dgemm( - CblasRowMajor, - CblasNoTrans, - CblasNoTrans, - n, - K-1, - m+1, - 1.0, - data->Z, - m+1, - model->V, - K-1, - 0.0, - ZV, - K-1); - - Memset(model->Q, double, n*K); - for (i=0; iUU, K-1, K, i, - j, k); - matrix_add(model->Q, K, i, k, value); - } - } - } -} - -/** - * @brief Calculate the Huber hinge errors - * - * @details - * For each of the scalar errors in Q the Huber hinge errors are - * calculated. The Huber hinge is here defined as - * @f[ - * h(q) = - * \begin{dcases} - * 1 - q - \frac{\kappa + 1}{2} & \text{if } q \leq -\kappa \\ - * \frac{1}{2(\kappa + 1)} ( 1 - q)^2 & \text{if } q \in (-\kappa, 1] \\ - * 0 & \text{if } q > 1 - * \end{dcases} - * @f] - * - * @param[in,out] model the corresponding MajModel - */ -void msvmmaj_calculate_huber(struct MajModel *model) -{ - long i, j; - double q, value; - - for (i=0; in; i++) { - for (j=0; jK; j++) { - q = matrix_get(model->Q, model->K, i, j); - value = 0.0; - if (q <= -model->kappa) { - value = 1.0 - q - (model->kappa+1.0)/2.0; - } else if (q <= 1.0) { - value = 1.0/(2.0*model->kappa+2.0)*pow(1.0 - q, 2.0); - } - matrix_set(model->H, model->K, i, j, value); - } - } -} - -/** - * @brief seed the matrix V from an existing model or using rand - * - * @details - * The matrix V must be seeded before the main_loop() can start. - * This can be done by either seeding it with random numbers or - * using the solution from a previous model on the same dataset - * as initial seed. The latter option usually allows for a - * significant improvement in the number of iterations necessary - * because the seeded model V is closer to the optimal V. - * - * @param[in] from_model MajModel from which to copy V - * @param[in,out] to_model MajModel to which V will be copied - */ -void msvmmaj_seed_model_V(struct MajModel *from_model, - struct MajModel *to_model, struct MajData *data) -{ - long i, j, k; - double cmin, cmax, value; - - long n = data->n; - long m = data->m; - long K = data->K; - - if (from_model == NULL) { - for (i=0; iZ, m+1, k, i); - cmin = minimum(cmin, value); - cmax = maximum(cmax, value); - } - for (j=0; jV, K-1, i, j, value); - } - } - } else { - for (i=0; iV, K-1, i, j); - matrix_set(to_model->V, K-1, i, j, value); - } - } -} - -/** - * @brief Use step doubling - * - * @details - * Step doubling can be used to speed up the Majorization algorithm. Instead - * of using the value at the minimimum of the majorization function, the value - * ``opposite'' the majorization point is used. This can essentially cut the - * number of iterations necessary to reach the minimum in half. - * - * @param[in] model MajModel containing the augmented parameters - */ -void msvmmaj_step_doubling(struct MajModel *model) -{ - long i, j; - double value; - - long m = model->m; - long K = model->K; - - for (i=0; iV, K-1, i, j, 2.0); - value = - matrix_get(model->Vbar, K-1, i, j); - matrix_add(model->V, K-1, i, j, value); - } - } -} - -/** - * @brief Initialize instance weights - * - * @details - * Instance weights can for instance be used to add additional weights to - * instances of certain classes. Two default weight possibilities are - * implemented here. The first is unit weights, where each instance gets - * weight 1. - * - * The second are group size correction weights, which are calculated as - * @f[ - * \rho_i = \frac{n}{Kn_k} , - * @f] - * where @f$ n_k @f$ is the number of instances in group @f$ k @f$ and - * @f$ y_i = k @f$. - * - * @param[in] data MajData with the dataset - * @param[in,out] model MajModel with the weight specification. On - * exit MajModel::rho contains the instance - * weights. - */ -void msvmmaj_initialize_weights(struct MajData *data, struct MajModel *model) -{ - long *groups; - long i; - - long n = model->n; - long K = model->K; - - if (model->weight_idx == 1) { - for (i=0; irho[i] = 1.0; - } - else if (model->weight_idx == 2) { - groups = Calloc(long, K); - for (i=0; iy[i]-1]++; - for (i=0; irho[i] = ((double) n)/((double) (groups[data->y[i]-1]*K)); - } else { - fprintf(stderr, "Unknown weight specification.\n"); - exit(1); - } -} - diff --git a/src/msvmmaj_init.c b/src/msvmmaj_init.c deleted file mode 100644 index 0fedfe7..0000000 --- a/src/msvmmaj_init.c +++ /dev/null @@ -1,282 +0,0 @@ -/** - * @file msvmmaj_init.c - * @author Gertjan van den Burg - * @date January 7, 2014 - * @brief Functions for initializing model and data structures - * - * @details - * This file contains functions for initializing a MajModel instance - * and a MajData instance. In addition, default values for these - * structures are defined here (and only here). Functions for allocating - * memory for the model structure and freeing of the model and data structures - * are also included. - * - */ - -#include - -#include "msvmmaj.h" -#include "msvmmaj_init.h" - -/** - * @brief Initialize a MajModel structure - * - * @details - * A MajModel structure is initialized and the default value for the - * parameters are set. A pointer to the initialized model is returned. - * - * @returns initialized MajModel - */ -struct MajModel *msvmmaj_init_model() -{ - struct MajModel *model = Malloc(struct MajModel, 1); - - // set default values - model->p = 1.0; - model->lambda = pow(2, -8.0); - model->epsilon = 1e-6; - model->kappa = 0.0; - model->weight_idx = 1; - model->kerneltype = K_LINEAR; - model->kernelparam = NULL; - - model->W = NULL; - model->t = NULL; - model->V = NULL; - model->Vbar = NULL; - model->U = NULL; - model->UU = NULL; - model->Q = NULL; - model->H = NULL; - model->R = NULL; - model->rho = NULL; - model->data_file = NULL; - - return model; -} - -/** - * @brief Initialize a MajData structure - * - * @details - * A MajData structure is initialized and default values are set. - * A pointer to the initialized data is returned. - * - * @returns initialized MajData - * - */ -struct MajData *msvmmaj_init_data() -{ - struct MajData *data = Malloc(struct MajData, 1); - data->J = NULL; - data->y = NULL; - data->Z = NULL; - data->RAW = NULL; - - // set default values - data->kerneltype = K_LINEAR; - data->kernelparam = NULL; - - return data; -} - -/** - * @brief Allocate memory for a MajModel - * - * @details - * This function can be used to allocate the memory needed for a MajModel. All - * arrays in the model are specified and initialized to 0. - * - * @param[in] model MajModel to allocate - * - */ -void msvmmaj_allocate_model(struct MajModel *model) -{ - long n = model->n; - long m = model->m; - long K = model->K; - - model->W = Calloc(double, m*(K-1)); - if (model->W == NULL) { - fprintf(stderr, "Failed to allocate memory for W.\n"); - exit(1); - } - - model->t = Calloc(double, K-1); - if (model->t == NULL) { - fprintf(stderr, "Failed to allocate memory for t.\n"); - exit(1); - } - - model->V = Calloc(double, (m+1)*(K-1)); - if (model->V == NULL) { - fprintf(stderr, "Failed to allocate memory for V.\n"); - exit(1); - } - - model->Vbar = Calloc(double, (m+1)*(K-1)); - if (model->Vbar == NULL) { - fprintf(stderr, "Failed to allocate memory for Vbar.\n"); - exit(1); - } - - model->U = Calloc(double, K*(K-1)); - if (model->U == NULL) { - fprintf(stderr, "Failed to allocate memory for U.\n"); - exit(1); - } - - model->UU = Calloc(double, n*K*(K-1)); - if (model->UU == NULL) { - fprintf(stderr, "Failed to allocate memory for UU.\n"); - exit(1); - } - - model->Q = Calloc(double, n*K); - if (model->Q == NULL) { - fprintf(stderr, "Failed to allocate memory for Q.\n"); - exit(1); - } - - model->H = Calloc(double, n*K); - if (model->H == NULL) { - fprintf(stderr, "Failed to allocate memory for H.\n"); - exit(1); - } - - model->R = Calloc(double, n*K); - if (model->R == NULL) { - fprintf(stderr, "Failed to allocate memory for R.\n"); - exit(1); - } - - model->rho = Calloc(double, n); - if (model->rho == NULL) { - fprintf(stderr, "Failed to allocate memory for rho.\n"); - exit(1); - } -} - -/** - * @brief Reallocate memory for MajModel - * - * @details - * This function can be used to reallocate existing memory for a MajModel, - * upon a change in the model dimensions. This is used in combination with - * kernels. - * - * @param[in] model MajModel to reallocate - * @param[in] n new value of MajModel->n - * @param[in] m new value of MajModel->m - * - */ -void msvmmaj_reallocate_model(struct MajModel *model, long n, long m) -{ - long K = model->K; - - if (model->n == n && model->m == m) - return; - if (model->n != n) { - model->UU = (double *) realloc(model->UU, - n*K*(K-1)*sizeof(double)); - if (model->UU == NULL) { - fprintf(stderr, "Failed to reallocate UU\n"); - exit(1); - } - - model->Q = (double *) realloc(model->Q, n*K*sizeof(double)); - if (model->Q == NULL) { - fprintf(stderr, "Failed to reallocate Q\n"); - exit(1); - } - - model->H = (double *) realloc(model->H, n*K*sizeof(double)); - if (model->H == NULL) { - fprintf(stderr, "Failed to reallocate H\n"); - exit(1); - } - - model->R = (double *) realloc(model->R, n*K*sizeof(double)); - if (model->R == NULL) { - fprintf(stderr, "Failed to reallocate R\n"); - exit(1); - } - - model->rho = (double *) realloc(model->rho, n*sizeof(double)); - if (model->rho == NULL) { - fprintf(stderr, "Failed to reallocte rho\n"); - exit(1); - } - - model->n = n; - } - if (model->m != m) { - model->W = (double *) realloc(model->W, - m*(K-1)*sizeof(double)); - if (model->W == NULL) { - fprintf(stderr, "Failed to reallocate W\n"); - exit(1); - } - - model->V = (double *) realloc(model->V, - (m+1)*(K-1)*sizeof(double)); - if (model->V == NULL) { - fprintf(stderr, "Failed to reallocate V\n"); - exit(1); - } - - model->Vbar = (double *) realloc(model->Vbar, - (m+1)*(K-1)*sizeof(double)); - if (model->Vbar == NULL) { - fprintf(stderr, "Failed to reallocate Vbar\n"); - exit(1); - } - - model->m = m; - } -} - -/** - * @brief Free allocated MajModel struct - * - * @details - * Simply free a previously allocated MajModel by freeing all its component - * arrays. Note that the model struct itself is also freed here. - * - * @param[in] model MajModel to free - * - */ -void msvmmaj_free_model(struct MajModel *model) -{ - free(model->W); - free(model->t); - free(model->V); - free(model->Vbar); - free(model->U); - free(model->UU); - free(model->Q); - free(model->H); - free(model->rho); - free(model->R); - free(model->kernelparam); - - free(model); -} - -/** - * @brief Free allocated MajData struct - * - * @details - * Simply free a previously allocated MajData struct by freeing all its - * components. Note that the data struct itself is also freed here. - * - * @param[in] data MajData struct to free - * - */ -void msvmmaj_free_data(struct MajData *data) -{ - free(data->Z); - free(data->y); - free(data->J); - free(data); -} diff --git a/src/msvmmaj_io.c b/src/msvmmaj_io.c deleted file mode 100644 index 8a09b3d..0000000 --- a/src/msvmmaj_io.c +++ /dev/null @@ -1,295 +0,0 @@ -/** - * @file msvmmaj_io.c - * @author Gertjan van den Burg - * @date January, 2014 - * @brief Functions for input and output of data and model files - * - * @details - * This file contains functions for reading and writing model files, and data - * files. - * - */ - -#include "msvmmaj.h" -#include "msvmmaj_io.h" -#include "msvmmaj_matrix.h" -#include "strutil.h" -#include "timer.h" - -/** - * @brief Read data from file - * - * @details - * Read the data from the data_file. The data matrix X is augmented - * with a column of ones, to get the matrix Z. The data is expected - * to follow a specific format, which is specified in the @ref spec_data_file. - * The class labels are corrected internally to correspond to the interval - * [1 .. K], where K is the total number of classes. - * - * @todo - * Make sure that this function allows datasets without class labels for - * testing. - * - * @param[in,out] dataset initialized MajData struct - * @param[in] data_file filename of the data file. - */ -void msvmmaj_read_data(struct MajData *dataset, char *data_file) -{ - FILE *fid; - long i, j; - long n, m; // dimensions of data - long nr = 0; // used to check consistency of data - double value; - long K = 0; - long min_y = 1000000; - - char buf[MAX_LINE_LENGTH]; - - if ((fid = fopen(data_file, "r")) == NULL) { - fprintf(stderr, "\nERROR: datafile %s could not be opened.\n", - data_file); - exit(0); - } - - // Read data dimensions - nr += fscanf(fid, "%ld", &n); - nr += fscanf(fid, "%ld", &m); - - // Allocate memory - dataset->RAW = Malloc(double, n*(m+1)); - - // Read first line of data - for (j=1; jRAW, n, 0, j, value); - } - - // Check if there is a label at the end of the line - if (fgets(buf, MAX_LINE_LENGTH, fid) == NULL) { - fprintf(stderr, "ERROR: No label found on first line.\n"); - exit(1); - } - if (sscanf(buf, "%lf", &value) > 0) { - dataset->y = Malloc(long, n); - dataset->y[0] = value; - } else if (dataset->y != NULL) { - free(dataset->y); - dataset->y = NULL; - } - - // Read the rest of the file - for (i=1; iRAW, m+1, i, j, value); - } - if (dataset->y != NULL) { - nr += fscanf(fid, "%lf", &value); - dataset->y[i] = (long) value; - K = maximum(K, value); - min_y = minimum(min_y, value); - } - } - fclose(fid); - - // Correct labels: must be in [1, K] - if (min_y == 0) { - for (i=0; iy[i]++; - K++; - } else if (min_y < 0 ) { - fprintf(stderr, "ERROR: wrong class labels in %s, minimum " - "value is: %ld\n", - data_file, min_y); - exit(0); - } - - if (nr < n * m) { - fprintf(stderr, "ERROR: not enough data found in %s\n", - data_file); - exit(0); - } - - // Set the column of ones - for (i=0; iRAW, m+1, i, 0, 1.0); - - dataset->n = n; - dataset->m = m; - dataset->K = K; - dataset->Z = dataset->RAW; -} - - -/** - * @brief Read model from file - * - * @details - * Read a MajModel from a model file. The MajModel struct must have been - * initalized elswhere. The model file is expected to follow the @ref - * spec_model_file. The easiest way to generate a model file is through - * msvmmaj_write_model(), which can for instance be used in trainMSVMMaj.c. - * - * @param[in,out] model initialized MajModel - * @param[in] model_filename filename of the model file - * - */ -void msvmmaj_read_model(struct MajModel *model, char *model_filename) -{ - long i, j, nr = 0; - FILE *fid; - char buffer[MAX_LINE_LENGTH]; - char data_filename[MAX_LINE_LENGTH]; - double value = 0; - - fid = fopen(model_filename, "r"); - if (fid == NULL) { - fprintf(stderr, "Error opening model file %s\n", - model_filename); - exit(1); - } - // skip the first four lines - for (i=0; i<4; i++) - next_line(fid, model_filename); - - // read all model variables - model->p = get_fmt_double(fid, model_filename, "p = %lf"); - model->lambda = get_fmt_double(fid, model_filename, "lambda = %lf"); - model->kappa = get_fmt_double(fid, model_filename, "kappa = %lf"); - model->epsilon = get_fmt_double(fid, model_filename, "epsilon = %lf"); - model->weight_idx = (int) get_fmt_long(fid, model_filename, - "weight_idx = %li"); - - // skip to data section - for (i=0; i<2; i++) - next_line(fid, model_filename); - - // read filename of data file - if (fgets(buffer, MAX_LINE_LENGTH, fid) == NULL) { - fprintf(stderr, "Error reading model file %s\n", - model_filename); - exit(1); - } - sscanf(buffer, "filename = %s\n", data_filename); - model->data_file = data_filename; - - // read all data variables - model->n = get_fmt_long(fid, model_filename, "n = %li\n"); - model->m = get_fmt_long(fid, model_filename, "m = %li\n"); - model->K = get_fmt_long(fid, model_filename, "K = %li\n"); - - // skip to output - for (i=0; i<2; i++) - next_line(fid, model_filename); - - // read the matrix V and check for consistency - model->V = Malloc(double, (model->m+1)*(model->K-1)); - for (i=0; im+1; i++) { - for (j=0; jK-1; j++) { - nr += fscanf(fid, "%lf ", &value); - matrix_set(model->V, model->K-1, i, j, value); - } - } - if (nr != (model->m+1)*(model->K-1)) { - fprintf(stderr, "Error reading model file %s. " - "Not enough elements of V found.\n", - model_filename); - exit(1); - } -} - -/** - * @brief Write model to file - * - * @details - * Write a MajModel to a file. The current time is specified in the file in - * UTC + offset. The model file further corresponds to the @ref - * spec_model_file. - * - * @param[in] model MajModel which contains an estimate for - * MajModel::V - * @param[in] output_filename the output file to write the model to - * - */ -void msvmmaj_write_model(struct MajModel *model, char *output_filename) -{ - FILE *fid; - long i, j; - char timestr[MAX_LINE_LENGTH]; - - // open output file - fid = fopen(output_filename, "w"); - if (fid == NULL) { - fprintf(stderr, "Error opening output file %s", - output_filename); - exit(1); - } - get_time_string(timestr); - - // Write output to file - fprintf(fid, "Output file for MSVMMaj (version %1.1f)\n", VERSION); - fprintf(fid, "Generated on: %s\n\n", timestr); - fprintf(fid, "Model:\n"); - fprintf(fid, "p = %15.16f\n", model->p); - fprintf(fid, "lambda = %15.16f\n", model->lambda); - fprintf(fid, "kappa = %15.16f\n", model->kappa); - fprintf(fid, "epsilon = %g\n", model->epsilon); - fprintf(fid, "weight_idx = %i\n", model->weight_idx); - fprintf(fid, "\n"); - fprintf(fid, "Data:\n"); - fprintf(fid, "filename = %s\n", model->data_file); - fprintf(fid, "n = %li\n", model->n); - fprintf(fid, "m = %li\n", model->m); - fprintf(fid, "K = %li\n", model->K); - fprintf(fid, "\n"); - fprintf(fid, "Output:\n"); - for (i=0; im+1; i++) { - for (j=0; jK-1; j++) { - fprintf(fid, "%+15.16f ", - matrix_get(model->V, - model->K-1, i, j)); - } - fprintf(fid, "\n"); - } - - fclose(fid); -} - -/** - * @brief Write predictions to file - * - * @details - * Write the given predictions to an output file, such that the resulting file - * corresponds to the @ref spec_data_file. - * - * @param[in] data MajData with the original instances - * @param[in] predy predictions of the class labels of the - * instances in the given MajData. Note that the - * order of the instances is assumed to be the - * same. - * @param[in] output_filename the file to which the predictions are written - * - */ -void msvmmaj_write_predictions(struct MajData *data, long *predy, - char *output_filename) -{ - long i, j; - FILE *fid; - - fid = fopen(output_filename, "w"); - if (fid == NULL) { - fprintf(stderr, "Error opening output file %s", - output_filename); - exit(1); - } - - for (i=0; in; i++) { - for (j=0; jm; j++) - fprintf(fid, "%f ", - matrix_get(data->Z, - data->m+1, i, j+1)); - fprintf(fid, "%li\n", predy[i]); - } - - fclose(fid); -} diff --git a/src/msvmmaj_kernel.c b/src/msvmmaj_kernel.c deleted file mode 100644 index 8f757c5..0000000 --- a/src/msvmmaj_kernel.c +++ /dev/null @@ -1,412 +0,0 @@ -/** - * @file msvmmaj_kernel.c - * @author Gertjan van den Burg - * @date October 18, 2013 - * @brief Defines main functions for use of kernels in MSVMMaj. - * - * @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 "msvmmaj.h" -#include "msvmmaj_kernel.h" -#include "msvmmaj_lapack.h" -#include "msvmmaj_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 MajModel specifying the parameters - * @param[in] data MajData specifying the data. - * - */ -void msvmmaj_make_kernel(struct MajModel *model, struct MajData *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 = msvmmaj_compute_poly(x1, x2, - model->kernelparam, data->m); - else if (model->kerneltype == K_RBF) - value = msvmmaj_compute_rbf(x1, x2, - model->kernelparam, data->m); - else if (model->kerneltype == K_SIGMOID) - value = msvmmaj_compute_sigmoid(x1, x2, - model->kernelparam, data->m); - else { - fprintf(stderr, "Unknown kernel type in " - "msvmmaj_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 = msvmmaj_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 msvmmaj_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 = msvmmaj_compute_poly(x1, x2, - model->kernelparam, - m); - else if (model->kerneltype == K_RBF) - value = msvmmaj_compute_rbf(x1, x2, - model->kernelparam, - m); - else if (model->kerneltype == K_SIGMOID) - value = msvmmaj_compute_sigmoid(x1, x2, - model->kernelparam, - m); - else { - fprintf(stderr, "Unknown kernel type in " - "msvmmaj_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 msvmmaj_compute_rbf(double *x1, double *x2, double *kernelparam, long n) -{ - long i; - double value = 0.0; - - for (i=0; i0: if i, the leading minor of A - * was not positive definite - * - * See the LAPACK documentation at: - * http://www.netlib.org/lapack/explore-html/dc/de9/group__double_p_osolve.html - */ -int dposv(char UPLO, int N, int NRHS, double *A, int LDA, double *B, - int LDB) -{ - extern void dposv_(char *UPLO, int *Np, int *NRHSp, double *A, - int *LDAp, double *B, int *LDBp, int *INFOp); - int INFO; - dposv_(&UPLO, &N, &NRHS, A, &LDA, B, &LDB, &INFO); - return INFO; -} - -/** - * @brief Solve a system of equations AX = B where A is symmetric. - * - * @details - * Solve a linear system of equations AX = B where A is symmetric. This - * function uses the external LAPACK routine dsysv. - * - * @param[in] UPLO which triangle of A is stored - * @param[in] N order of A - * @param[in] NRHS number of columns of B - * @param[in,out] A double precision array of size (LDA, N). On - * exit contains the block diagonal matrix D and - * the multipliers used to obtain the factor U or - * L from the factorization A = U*D*U**T or - * A = L*D*L**T. - * @param[in] LDA leading dimension of A - * @param[in] IPIV integer array containing the details of D - * @param[in,out] B double precision array of size (LDB, NRHS). On - * exit contains the N-by-NRHS matrix X - * @param[in] LDB leading dimension of B - * @param[out] WORK double precision array of size max(1,LWORK). On - * exit, WORK(1) contains the optimal LWORK - * @param[in] LWORK the length of WORK, can be used for determining - * the optimal blocksize for dsystrf. - * @returns info parameter which contains the status of the - * computation: - * - =0: success - * - <0: if -i, the i-th argument had an - * illegal value - * - >0: if i, D(i, i) is exactly zero, - * no solution can be computed. - * - * See the LAPACK documentation at: - * http://www.netlib.org/lapack/explore-html/d6/d0e/group__double_s_ysolve.html - */ -int dsysv(char UPLO, int N, int NRHS, double *A, int LDA, int *IPIV, - double *B, int LDB, double *WORK, int LWORK) -{ - extern void dsysv_(char *UPLO, int *Np, int *NRHSp, double *A, - int *LDAp, int *IPIV, double *B, int *LDBp, - double *WORK, int *LWORK, int *INFOp); - int INFO; - dsysv_(&UPLO, &N, &NRHS, A, &LDA, IPIV, B, &LDB, WORK, &LWORK, &INFO); - return INFO; -} - -/** - * @brief Compute the eigenvalues and optionally the eigenvectors of a - * symmetric matrix. - * - * @details - * See the LAPACK documentation at: - * http://www.netlib.org/lapack/explore-html/d2/d97/dsyevx_8f.html - * - * - */ -int dsyevx(char JOBZ, char RANGE, char UPLO, int N, double *A, int LDA, - double VL, double VU, int IL, int IU, double ABSTOL, int *M, - double *W, double *Z, int LDZ, double *WORK, int LWORK, - int *IWORK, int *IFAIL) -{ - extern void dsyevx_(char *JOBZ, char *RANGE, char *UPLO, int *Np, - double *A, int *LDAp, double *VLp, double *VUp, - int *ILp, int *IUp, double *ABSTOLp, int *M, - double *W, double *Z, int *LDZp, double *WORK, - int *LWORKp, int *IWORK, int *IFAIL, int *INFOp); - int INFO; - dsyevx_(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, - M, W, Z, &LDZ, WORK, &LWORK, IWORK, IFAIL, &INFO); - return INFO; -} - -/** - * @brief Determine double precision machine parameters. - * - * @details - * See the LAPACK documentation at: - * http://www.netlib.org/lapack/explore-html/d5/dd4/dlamch_8f.html - */ -double dlamch(char CMACH) -{ - extern double dlamch_(char *CMACH); - return dlamch_(&CMACH); -} diff --git a/src/msvmmaj_matrix.c b/src/msvmmaj_matrix.c deleted file mode 100644 index 9e1be04..0000000 --- a/src/msvmmaj_matrix.c +++ /dev/null @@ -1,38 +0,0 @@ -/** - * @file msvmmaj_matrix.c - * @author Gertjan van den Burg - * @date August, 2013 - * @brief Functions facilitating matrix access - * - * @details - * The functions contained in this file are used when - * accessing or writing to matrices. Seperate functions - * exist of adding and multiplying existing matrix - * elements, to ensure this is done in place. - * - */ - -#include "msvmmaj_matrix.h" -#include "util.h" - -/** - * @brief print a matrix - * - * @details - * Debug function to print a matrix - * - * @param[in] M matrix - * @param[in] rows number of rows of M - * @param[in] cols number of columns of M - */ -void print_matrix(double *M, long rows, long cols) -{ - long i, j; - - for (i=0; i - -#include "libMSVMMaj.h" -#include "msvmmaj.h" -#include "msvmmaj_kernel.h" -#include "msvmmaj_matrix.h" -#include "msvmmaj_pred.h" - -#include "util.h" // testing - -void msvmmaj_predict_labels(struct MajData *data_test, - struct MajData *data_train, struct MajModel *model, - long *predy) -{ - if (model->kerneltype == K_LINEAR) - msvmmaj_predict_labels_linear(data_test, model, predy); - else - msvmmaj_predict_labels_kernel(data_test, data_train, model, - predy); -} - -/** - * @brief Predict class labels of data given and output in predy - * - * @details - * The labels are predicted by mapping each instance in data to the - * simplex space using the matrix V in the given model. Next, for each - * instance the nearest simplex vertex is determined using an Euclidean - * norm. The nearest simplex vertex determines the predicted class label, - * which is recorded in predy. - * - * @param[in] data MajData to predict labels for - * @param[in] model MajModel with optimized V - * @param[out] predy pre-allocated vector to record predictions in - */ -void msvmmaj_predict_labels_linear(struct MajData *data, - struct MajModel *model, long *predy) -{ - long i, j, k, label; - double norm, min_dist; - - long n = data->n; // note that model->n is the size of the training sample. - long m = data->m; - long K = model->K; //data->K does not necessarily equal the original K. - - double *S = Calloc(double, K-1); - double *ZV = Calloc(double, n*(K-1)); - double *U = Calloc(double, K*(K-1)); - - // Get the simplex matrix - msvmmaj_simplex_gen(K, U); - - // Generate the simplex-space vectors - cblas_dgemm( - CblasRowMajor, - CblasNoTrans, - CblasNoTrans, - n, - K-1, - m+1, - 1.0, - data->Z, - m+1, - model->V, - K-1, - 0.0, - ZV, - K-1); - - // Calculate the distance to each of the vertices of the simplex. - // The closest vertex defines the class label. - for (i=0; in; - long n_test = data_test->n; - long r = model->m; - long K = model->K; - - double *K2 = NULL; - msvmmaj_make_crosskernel(model, data_train, data_test, &K2); - - double *S = Calloc(double, K-1); - double *ZV = Calloc(double, n_test*(r+1)); - double *KPS = Calloc(double, n_test*(r+1)); - double *U = Calloc(double, K*(K-1)); - - msvmmaj_simplex_gen(K, U); - - // were doing the computations explicitly since P is included in - // data_train->Z. Might want to look at this some more if it turns out - // to be slow. - - double value, rowvalue; - for (i=0; iZ, r+1, k, - j); - value += rowvalue; - } - value *= matrix_get(data_train->J, 1, j, 0); - matrix_set(KPS, r+1, i, j, value); - } - matrix_set(KPS, r+1, i, 0, 1.0); - } - - cblas_dgemm( - CblasRowMajor, - CblasNoTrans, - CblasNoTrans, - n_test, - K-1, - r+1, - 1.0, - KPS, - r+1, - model->V, - K-1, - 0.0, - ZV, - K-1); - - for (i=0; in; i++) - if (data->y[i] == predy[i]) - correct++; - - performance = ((double) correct)/((double) data->n)* 100.0; - - return performance; -} diff --git a/src/msvmmaj_sv.c b/src/msvmmaj_sv.c deleted file mode 100644 index 1358d4e..0000000 --- a/src/msvmmaj_sv.c +++ /dev/null @@ -1,45 +0,0 @@ -/** - * @file msvmmaj_sv.c - * @author Gertjan van den Burg - * @date May, 2014 - * @brief Calculate the number of support vectors - * - * @details - * The function in this file can be used to calculate the number of support - * vectors are left in a model. - * - */ - -#include "msvmmaj.h" -#include "msvmmaj_matrix.h" - -/** - * @brief Calculate the number of support vectors in a model - * - * @details - * If an object is correctly classified, the number of classes for which the - * error q is larger than 1, is K-1 (i.e., there is no error w.r.t. any of the - * other classes). All objects for which this is not the case are thus support - * vectors. - * - * @param[in] model MajModel with solution - * @param[in] data MajData to be used - * @return number of support vectors with this solution - * - */ -long msvmmaj_num_sv(struct MajModel *model, struct MajData *data) -{ - long i, j, num_correct, num_sv = 0; - double value; - - for (i=0; in; i++) { - num_correct = 0; - for (j=0; jK; j++) { - value = matrix_get(model->Q, data->K, i, j); - num_correct += (value > 1); - } - num_sv += (num_correct < data->K - 1); - } - - return num_sv; -} diff --git a/src/msvmmaj_train.c b/src/msvmmaj_train.c deleted file mode 100644 index 09b00ee..0000000 --- a/src/msvmmaj_train.c +++ /dev/null @@ -1,534 +0,0 @@ -/** - * @file msvmmaj_train.c - * @author Gertjan van den Burg - * @date August 9, 2013 - * @brief Main functions for training the MSVMMaj solution. - * - * @details - * Contains update and loss functions used to actually find - * the optimal V. - * - */ - -#include -#include - -#include "libMSVMMaj.h" -#include "msvmmaj.h" -#include "msvmmaj_lapack.h" -#include "msvmmaj_matrix.h" -#include "msvmmaj_sv.h" -#include "msvmmaj_train.h" -#include "util.h" - -/** - * Maximum number of iterations of the algorithm. - */ -#define MAX_ITER 1000000000 - -/** - * @brief The main training loop for MSVMMaj - * - * @details - * This function is the main training function. This function - * handles the optimization of the model with the given model parameters, with - * the data given. On return the matrix MajModel::V contains the optimal - * weight matrix. - * - * In this function, step doubling is used in the majorization algorithm after - * a burn-in of 50 iterations. If the training is finished, MajModel::t and - * MajModel::W are extracted from MajModel::V. - * - * @param[in,out] model the MajModel to be trained. Contains optimal - * V on exit. - * @param[in] data the MajData to train the model with. - */ -void msvmmaj_optimize(struct MajModel *model, struct MajData *data) -{ - long i, j, it = 0; - double L, Lbar, value; - - long n = model->n; - long m = model->m; - long K = model->K; - - double *B = Calloc(double, n*(K-1)); - double *ZV = Calloc(double, n*(K-1)); - double *ZAZ = Calloc(double, (m+1)*(m+1)); - double *ZAZV = Calloc(double, (m+1)*(K-1)); - double *ZAZVT = Calloc(double, (m+1)*(K-1)); - - note("Starting main loop.\n"); - note("Dataset:\n"); - note("\tn = %i\n", n); - note("\tm = %i\n", m); - note("\tK = %i\n", K); - note("Parameters:\n"); - note("\tkappa = %f\n", model->kappa); - note("\tp = %f\n", model->p); - note("\tlambda = %15.16f\n", model->lambda); - note("\tepsilon = %g\n", model->epsilon); - note("\n"); - - msvmmaj_simplex_gen(model->K, model->U); - msvmmaj_simplex_diff(model, data); - msvmmaj_category_matrix(model, data); - - L = msvmmaj_get_loss(model, data, ZV); - Lbar = L + 2.0*model->epsilon*L; - - while ((it < MAX_ITER) && (Lbar - L)/L > model->epsilon) - { - // ensure V contains newest V and Vbar contains V from - // previous - msvmmaj_get_update(model, data, B, ZAZ, ZAZV, ZAZVT); - if (it > 50) - msvmmaj_step_doubling(model); - - Lbar = L; - L = msvmmaj_get_loss(model, data, ZV); - - if (it%100 == 0) - note("iter = %li, L = %15.16f, Lbar = %15.16f, " - "reldiff = %15.16f\n", it, L, Lbar, (Lbar - L)/L); - it++; - } - if (L > Lbar) - fprintf(stderr, "Negative step occurred in majorization.\n"); - - note("optimization finished, iter = %li, loss = %15.16f, " - "rel. diff. = %15.16f\n", it-1, L, - (Lbar - L)/L); - note("number of support vectors: %li\n", msvmmaj_num_sv(model, data)); - - model->training_error = (Lbar - L)/L; - - for (i=0; it[i] = matrix_get(model->V, K-1, 0, i); - for (i=1; iV, K-1, i, j); - matrix_set(model->W, K-1, i-1, j, value); - } - } - free(B); - free(ZV); - free(ZAZ); - free(ZAZV); - free(ZAZVT); -} - -/** - * @brief Calculate the current value of the loss function - * - * @details - * The current loss function value is calculated based on the matrix V in the - * given model. Note that the matrix ZV is passed explicitly to avoid having - * to reallocate memory at every step. - * - * @param[in] model MajModel structure which holds the current - * estimate V - * @param[in] data MajData structure - * @param[in,out] ZV pre-allocated matrix ZV which is updated on - * output - * @returns the current value of the loss function - */ -double msvmmaj_get_loss(struct MajModel *model, struct MajData *data, - double *ZV) -{ - long i, j; - long n = data->n; - long K = data->K; - long m = data->m; - - double value, rowvalue, loss = 0.0; - - msvmmaj_calculate_errors(model, data, ZV); - msvmmaj_calculate_huber(model); - - for (i=0; iH, K, i, j); - value = pow(value, model->p); - value *= matrix_get(model->R, K, i, j); - rowvalue += value; - } - rowvalue = pow(rowvalue, 1.0/(model->p)); - rowvalue *= model->rho[i]; - loss += rowvalue; - } - loss /= ((double) n); - - value = 0; - for (i=0; iV, K-1, i, j), 2.0); - } - value += data->J[i] * rowvalue; - } - loss += model->lambda * value; - - return loss; -} - -/** - * @brief Perform a single step of the majorization algorithm to update V - * - * @details - * This function contains the main update calculations of the algorithm. These - * calculations are necessary to find a new update V. The calculations exist of - * recalculating the majorization coefficients for all instances and all - * classes, and solving a linear system to find V. - * - * Because the function msvmmaj_get_update() is always called after a call to - * msvmmaj_get_loss() with the same MajModel::V, it is unnecessary to calculate - * the updated errors MajModel::Q and MajModel::H here too. This saves on - * computation time. - * - * In calculating the majorization coefficients we calculate the elements of a - * diagonal matrix A with elements - * @f[ - * A_{i, i} = \frac{1}{n} \rho_i \sum_{j \neq k} \left[ - * \varepsilon_i a_{ijk}^{(p)} + (1 - \varepsilon_i) \omega_i - * a_{ijk}^{(p)} \right], - * @f] - * where @f$ k = y_i @f$. - * Since this matrix is only used to calculate the matrix @f$ Z' A Z @f$, it is - * efficient to update a matrix ZAZ through consecutive rank 1 updates with - * a single element of A and the corresponding row of Z. The BLAS function - * dsyr is used for this. - * - * The B matrix is has rows - * @f[ - * \boldsymbol{\beta}_i' = \frac{1}{n} \rho_i \sum_{j \neq k} \left[ - * \varepsilon_i \left( b_{ijk}^{(1)} - a_{ijk}^{(1)} - * \overline{q}_i^{(kj)} \right) + (1 - \varepsilon_i) - * \omega_i \left( b_{ijk}^{(p)} - a_{ijk}^{(p)} - * \overline{q}_i^{(kj)} \right) \right] - * \boldsymbol{\delta}_{kj}' - * @f] - * This is also split into two cases, one for which @f$ \varepsilon_i = 1 @f$, - * and one for when it is 0. The 3D simplex difference matrix is used here, in - * the form of the @f$ \boldsymbol{\delta}_{kj}' @f$. - * - * Finally, the following system is solved - * @f[ - * (\textbf{Z}'\textbf{AZ} + \lambda \textbf{J})\textbf{V} = - * (\textbf{Z}'\textbf{AZ}\overline{\textbf{V}} + \textbf{Z}' - * \textbf{B}) - * @f] - * solving this system is done through dposv(). - * - * @todo - * Consider allocating IPIV and WORK at a higher level, they probably don't - * change much during the iterations. - * - * @param [in,out] model model to be updated - * @param [in] data data used in model - * @param [in] B pre-allocated matrix used for linear coefficients - * @param [in] ZAZ pre-allocated matrix used in system - * @param [in] ZAZV pre-allocated matrix used in system solving - * @param [in] ZAZVT pre-allocated matrix used in system solving - */ -void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B, - double *ZAZ, double *ZAZV, double *ZAZVT) -{ - int status, class; - long i, j, k; - double Avalue, Bvalue; - double omega, value, a, b, q, h, r; - - long n = model->n; - long m = model->m; - long K = model->K; - - double kappa = model->kappa; - double p = model->p; - double *rho = model->rho; - - // constants which are used often throughout - const double a2g2 = 0.25*p*(2.0*p - 1.0)*pow((kappa+1.0)/2.0,p-2.0); - const double in = 1.0/((double) n); - - // clear matrices - Memset(B, double, n*(K-1)); - Memset(ZAZ, double, (m+1)*(m+1)); - - b = 0; - for (i=0; iH, K, i, j); - r = matrix_get(model->R, K, i, j); - value += (h*r > 0) ? 1 : 0; - omega += pow(h, p)*r; - } - class = (value <= 1.0) ? 1 : 0; - omega = (1.0/p)*pow(omega, 1.0/p - 1.0); - - Avalue = 0; - if (class == 1) { - for (j=0; jQ, K, i, j); - if (q <= -kappa) { - a = 0.25/(0.5 - kappa/2.0 - q); - b = 0.5; - } else if (q <= 1.0) { - a = 1.0/(2.0*kappa + 2.0); - b = (1.0 - q)*a; - } else { - a = -0.25/(0.5 - kappa/2.0 - q); - b = 0; - } - for (k=0; kUU, K-1, K, i, k, j); - matrix_add(B, K-1, i, k, Bvalue); - } - Avalue += a*matrix_get(model->R, K, i, j); - } - } else { - if (2.0 - p < 0.0001) { - for (j=0; jQ, K, i, j); - if (q <= -kappa) { - b = 0.5 - kappa/2.0 - q; - } else if ( q <= 1.0) { - b = pow(1.0 - q, 3.0)/( - 2.0*pow(kappa + 1.0, - 2.0)); - } else { - b = 0; - } - for (k=0; kUU, - K-1, - K, - i, - k, - j); - matrix_add( - B, - K-1, - i, - k, - Bvalue); - } - } - Avalue = 1.5*(K - 1.0); - } else { - for (j=0; jQ, K, i, j); - if (q <= (p + kappa - 1.0)/(p - 2.0)) { - a = 0.25*pow(p, 2.0)*pow( - 0.5 - kappa/2.0 - q, - p - 2.0); - } else if (q <= 1.0) { - a = a2g2; - } else { - a = 0.25*pow(p, 2.0)*pow( - (p/(p - 2.0))* - (0.5 - kappa/2.0 - q), - p - 2.0); - b = a*(2.0*q + kappa - 1.0)/ - (p - 2.0) + - 0.5*p*pow( - p/(p - 2.0)* - (0.5 - kappa/ - 2.0 - q), - p - 1.0); - } - if (q <= -kappa) { - b = 0.5*p*pow( - 0.5 - kappa/2.0 - q, - p - 1.0); - } else if ( q <= 1.0) { - b = p*pow(1.0 - q, - 2.0*p - 1.0)/ - pow(2*kappa+2.0, p); - } - for (k=0; kUU, - K-1, - K, - i, - k, - j); - matrix_add( - B, - K-1, - i, - k, - Bvalue); - } - Avalue += a*matrix_get(model->R, - K, i, j); - } - } - Avalue *= omega; - } - Avalue *= in * rho[i]; - - // Now we calculate the matrix ZAZ. Since this is - // guaranteed to be symmetric, we only calculate the - // upper part of the matrix, and then copy this over - // to the lower part after all calculations are done. - // Note that the use of dsym is faster than dspr, even - // though dspr uses less memory. - cblas_dsyr( - CblasRowMajor, - CblasUpper, - m+1, - Avalue, - &data->Z[i*(m+1)], - 1, - ZAZ, - m+1); - } - // Copy upper to lower (necessary because we need to switch - // to Col-Major order for LAPACK). - /* - for (i=0; iV, - K-1, - 0.0, - ZAZV, - K-1); - - cblas_dgemm( - CblasRowMajor, - CblasTrans, - CblasNoTrans, - m+1, - K-1, - n, - 1.0, - data->Z, - m+1, - B, - K-1, - 1.0, - ZAZV, - K-1); - - /* - * Add lambda to all diagonal elements except the first one. Recall - * that ZAZ is of size m+1 and is symmetric. - */ - i = 0; - for (j=0; jlambda * data->J[j+1]; - } - - // For the LAPACK call we need to switch to Column- - // Major order. This is unnecessary for the matrix - // ZAZ because it is symmetric. The matrix ZAZV - // must be converted however. - for (i=0; iVbar; - model->Vbar = model->V; - model->V = ZAZVT; - ZAZVT = ptr; - */ - - for (i=0; iV, K-1, i, j); - matrix_set(model->Vbar, K-1, i, j, value); - value = matrix_get(ZAZV, K-1, i, j); - matrix_set(model->V, K-1, i, j, value); - } - } -} diff --git a/src/msvmmaj_train_dataset.c b/src/msvmmaj_train_dataset.c deleted file mode 100644 index 26c684c..0000000 --- a/src/msvmmaj_train_dataset.c +++ /dev/null @@ -1,748 +0,0 @@ -/** - * @file msvmmaj_train_dataset.c - * @author Gertjan van den Burg - * @date January, 2014 - * @brief Functions for finding the optimal parameters for the dataset - * - * @details - * The MSVMMaj algorithm takes a number of parameters. The functions in - * this file are used to find the optimal parameters. - */ - -#include -#include - -#include "crossval.h" -#include "libMSVMMaj.h" -#include "msvmmaj.h" -#include "msvmmaj_init.h" -#include "msvmmaj_kernel.h" -#include "msvmmaj_matrix.h" -#include "msvmmaj_train.h" -#include "msvmmaj_train_dataset.h" -#include "msvmmaj_pred.h" -#include "util.h" -#include "timer.h" - -extern FILE *MSVMMAJ_OUTPUT_FILE; - -/** - * @brief Initialize a Queue from a Training instance - * - * @details - * A Training instance describes the grid to search over. This funtion - * creates all tasks that need to be performed and adds these to - * a Queue. Each task contains a pointer to the train and test datasets - * which are supplied. Note that the tasks are created in a specific order of - * the parameters, to ensure that the MajModel::V of a previous parameter - * set provides the best possible initial estimate of MajModel::V for the next - * parameter set. - * - * @param[in] training Training struct describing the grid search - * @param[in] queue pointer to a Queue that will be used to - * add the tasks to - * @param[in] train_data MajData of the training set - * @param[in] test_data MajData of the test set - * - */ -void make_queue(struct Training *training, struct Queue *queue, - struct MajData *train_data, struct MajData *test_data) -{ - long i, j, k; - long N, cnt = 0; - struct Task *task; - queue->i = 0; - - N = training->Np; - N *= training->Nl; - N *= training->Nk; - N *= training->Ne; - N *= training->Nw; - // these parameters are not necessarily non-zero - N *= training->Ng > 0 ? training->Ng : 1; - N *= training->Nc > 0 ? training->Nc : 1; - N *= training->Nd > 0 ? training->Nd : 1; - - queue->tasks = Malloc(struct Task *, N); - queue->N = N; - - // initialize all tasks - for (i=0; iID = i; - task->train_data = train_data; - task->test_data = test_data; - task->folds = training->folds; - task->kerneltype = training->kerneltype; - task->kernelparam = Calloc(double, training->Ng + - training->Nc + training->Nd); - queue->tasks[i] = task; - } - - // These loops mimick a large nested for loop. The advantage is that - // Nd, Nc and Ng which are on the outside of the nested for loop can - // now be zero, without large modification (see below). Whether this - // is indeed better than the nested for loop has not been tested. - cnt = 1; - i = 0; - while (i < N ) - for (j=0; jNp; j++) - for (k=0; ktasks[i]->p = training->ps[j]; - i++; - } - - cnt *= training->Np; - i = 0; - while (i < N ) - for (j=0; jNl; j++) - for (k=0; ktasks[i]->lambda = - training->lambdas[j]; - i++; - } - - cnt *= training->Nl; - i = 0; - while (i < N ) - for (j=0; jNk; j++) - for (k=0; ktasks[i]->kappa = training->kappas[j]; - i++; - } - - cnt *= training->Nk; - i = 0; - while (i < N ) - for (j=0; jNw; j++) - for (k=0; ktasks[i]->weight_idx = - training->weight_idxs[j]; - i++; - } - - cnt *= training->Nw; - i = 0; - while (i < N ) - for (j=0; jNe; j++) - for (k=0; ktasks[i]->epsilon = - training->epsilons[j]; - i++; - } - - cnt *= training->Ne; - i = 0; - while (i < N && training->Ng > 0) - for (j=0; jNg; j++) - for (k=0; ktasks[i]->kernelparam[0] = - training->gammas[j]; - i++; - } - - cnt *= training->Ng > 0 ? training->Ng : 1; - i = 0; - while (i < N && training->Nc > 0) - for (j=0; jNc; j++) - for (k=0; ktasks[i]->kernelparam[1] = - training->coefs[j]; - i++; - } - - cnt *= training->Nc > 0 ? training->Nc : 1; - i = 0; - while (i < N && training->Nd > 0) - for (j=0; jNd; j++) - for (k=0; ktasks[i]->kernelparam[2] = - training->degrees[j]; - i++; - } -} - -/** - * @brief Get new Task from Queue - * - * @details - * Return a pointer to the next Task in the Queue. If no Task instances are - * left, NULL is returned. The internal counter Queue::i is used for finding - * the next Task. - * - * @param[in] q Queue instance - * @returns pointer to next Task - * - */ -struct Task *get_next_task(struct Queue *q) -{ - long i = q->i; - if (i < q->N) { - q->i++; - return q->tasks[i]; - } - return NULL; -} - -/** - * @brief Comparison function for Tasks based on performance - * - * @details - * To be able to sort Task structures on the performance of their specific - * set of parameters, this comparison function is implemented. Task structs - * are sorted with highest performance first. - * - * @param[in] elem1 Task 1 - * @param[in] elem2 Task 2 - * @returns result of inequality of Task 1 performance over - * Task 2 performance - */ -int tasksort(const void *elem1, const void *elem2) -{ - const struct Task *t1 = (*(struct Task **) elem1); - const struct Task *t2 = (*(struct Task **) elem2); - return (t1->performance > t2->performance); -} - -/** - * @brief Comparison function for doubles - * - * @details - * Similar to tasksort() only now for two doubles. - * - * @param[in] elem1 number 1 - * @param[in] elem2 number 2 - * @returns comparison of number 1 larger than number 2 - */ -int doublesort(const void *elem1, const void *elem2) -{ - const double t1 = (*(double *) elem1); - const double t2 = (*(double *) elem2); - return t1 > t2; -} - -/** - * @brief Calculate the percentile of an array of doubles - * - * @details - * The percentile of performance is used to find the top performing - * configurations. Since no standard definition of the percentile exists, we - * use the method used in MATLAB and Octave. Since calculating the percentile - * requires a sorted list of the values, a local copy is made first. - * - * @param[in] values array of doubles - * @param[in] N length of the array - * @param[in] p percentile to calculate ( 0 <= p <= 1.0 ). - * @returns the p-th percentile of the values - */ -double prctile(double *values, long N, double p) -{ - long i; - double pi, pr, boundary; - double *local = Malloc(double, N); - for (i=0; itasks, q->N, sizeof(struct Task *), tasksort); - p = 0.95*q->N + 0.5; - pi = maximum(minimum(floor(p), q->N-1), 1); - pr = maximum(minimum(p - pi, 1), 0); - boundary = (1 - pr)*q->tasks[((long) pi)-1]->performance; - boundary += pr*q->tasks[((long) pi)]->performance; - note("boundary determined at: %f\n", boundary); - - // find the number of tasks that perform at least as good as the 95th - // percentile - N = 0; - for (i=0; iN; i++) - if (q->tasks[i]->performance >= boundary) - N++; - note("Number of items: %li\n", N); - std = Calloc(double, N); - mean = Calloc(double, N); - time = Calloc(double, N); - perf = Calloc(double, N*repeats); - - // create a new task queue with the tasks which perform well - nq->tasks = Malloc(struct Task *, N); - for (i=q->N-1; i>q->N-N-1; i--) - nq->tasks[q->N-i-1] = q->tasks[i]; - nq->N = N; - nq->i = 0; - - // for each task run the consistency repeats - for (i=0; in = 0; - model->m = task->train_data->m; - model->K = task->train_data->K; - msvmmaj_allocate_model(model); - msvmmaj_seed_model_V(NULL, model, task->train_data); - } - - time[i] = 0.0; - note("(%02li/%02li:%03li)\t", i+1, N, task->ID); - for (r=0; rtrain_data, - task->folds); - loop_e = clock(); - time[i] += elapsed_time(loop_s, loop_e); - matrix_set(perf, repeats, i, r, p); - mean[i] += p/((double) repeats); - } else { - note("Only cv is implemented\n"); - exit(1); - } - note("%3.3f\t", p); - // this is done because if we reuse the V it's not a - // consistency check - msvmmaj_seed_model_V(NULL, model, task->train_data); - } - for (r=0; r 1) { - std[i] /= ((double) repeats) - 1.0; - std[i] = sqrt(std[i]); - } else - std[i] = 0.0; - note("(m = %3.3f, s = %3.3f, t = %3.3f)\n", - mean[i], std[i], time[i]); - } - - // find the best overall configurations: those with high average - // performance and low deviation in the performance - note("\nBest overall configuration(s):\n"); - note("ID\tweights\tepsilon\t\tp\t\tkappa\t\tlambda\t\t" - "mean_perf\tstd_perf\ttime_perf\n"); - p = 0.0; - bool breakout = false; - while (breakout == false) { - pi = prctile(mean, N, (100.0-p)/100.0); - pr = prctile(std, N, p/100.0); - pt = prctile(time, N, p/100.0); - for (i=0; itasks[i]->ID, - nq->tasks[i]->weight_idx, - nq->tasks[i]->epsilon, - nq->tasks[i]->p, - nq->tasks[i]->kappa, - nq->tasks[i]->lambda, - mean[i], - std[i], - time[i]); - breakout = true; - } - p += 1.0; - } - - free(nq->tasks); - free(nq); - free(model); - free(perf); - free(std); - free(mean); - free(time); -} - -/** - * @brief Run cross validation with a seed model - * - * @details - * This is an implementation of cross validation which uses the optimal - * parameters MajModel::V of a previous fold as initial conditions for - * MajModel::V of the next fold. An initial seed for V can be given through the - * seed_model parameter. If seed_model is NULL, random starting values are - * used. - * - * @param[in] model MajModel with the configuration to train - * @param[in] seed_model MajModel with a seed for MajModel::V - * @param[in] data MajData with the dataset - * @param[in] folds number of cross validation folds - * @returns performance (hitrate) of the configuration on - * cross validation - */ -double cross_validation(struct MajModel *model, struct MajData *data, - long folds) -{ - FILE *fid; - - long f, *predy; - double performance, total_perf = 0; - struct MajData *train_data, *test_data; - - long *cv_idx = Calloc(long, data->n); - - train_data = msvmmaj_init_data(); - test_data = msvmmaj_init_data(); - - // create splits - msvmmaj_make_cv_split(data->n, folds, cv_idx); - - for (f=0; fn, train_data->m); - - msvmmaj_initialize_weights(train_data, model); - - // train the model (without output) - fid = MSVMMAJ_OUTPUT_FILE; - MSVMMAJ_OUTPUT_FILE = NULL; - msvmmaj_optimize(model, train_data); - MSVMMAJ_OUTPUT_FILE = fid; - - // calculate prediction performance on test set - predy = Calloc(long, test_data->n); - msvmmaj_predict_labels(test_data, train_data, model, predy); - performance = msvmmaj_prediction_perf(test_data, predy); - total_perf += performance * test_data->n; - - free(predy); - free(train_data->y); - free(train_data->Z); - free(test_data->y); - free(test_data->Z); - } - - free(train_data); - free(test_data); - - total_perf /= ((double) data->n); - - return total_perf; -} - -/** - * @brief Run the grid search for a cross validation dataset - * - * @details - * Given a Queue of Task struct to be trained, a grid search is launched to - * find the optimal parameter configuration. As is also done within - * cross_validation(), the optimal weights of one parameter set are used as - * initial estimates for MajModel::V in the next parameter set. Note that to - * optimally exploit this feature of the optimization algorithm, the order in - * which tasks are considered is important. This is considered in - * make_queue(). - * - * The performance found by cross validation is stored in the Task struct. - * - * @param[in,out] q Queue with Task instances to run - */ -void start_training_cv(struct Queue *q) -{ - double perf, current_max = 0; - struct Task *task = get_next_task(q); - struct MajModel *model = msvmmaj_init_model(); - clock_t main_s, main_e, loop_s, loop_e; - - model->n = 0; - model->m = task->train_data->m; - model->K = task->train_data->K; - msvmmaj_allocate_model(model); - msvmmaj_seed_model_V(NULL, model, task->train_data); - - main_s = clock(); - while (task) { - print_progress_string(task, q->N); - make_model_from_task(task, model); - - loop_s = clock(); - perf = cross_validation(model, task->train_data, task->folds); - loop_e = clock(); - current_max = maximum(current_max, perf); - - note("\t%3.3f%% (%3.3fs)\t(best = %3.3f%%)\n", perf, - elapsed_time(loop_s, loop_e), - current_max); - - q->tasks[task->ID]->performance = perf; - task = get_next_task(q); - } - main_e = clock(); - - note("\nTotal elapsed time: %8.8f seconds\n", - elapsed_time(main_s, main_e)); - - msvmmaj_free_model(model); -} - -/** - * @brief Run the grid search for a train/test dataset - * - * @details - * This function is similar to start_training_cv(), except that the - * pre-determined training set is used only once, and the pre-determined test - * set is used for validation. - * - * @todo - * It would probably be better to train the model on the training set using - * cross validation and only use the test set when comparing with other - * methods. The way it is now, you're finding out which parameters predict - * _this_ test set best, which is not what you want. This function should - * therefore not be used and is considered deprecated, to be removed in the - * future . - * - * @param[in] q Queue with Task structs to run - * - */ -void start_training_tt(struct Queue *q) -{ - FILE *fid; - - long c = 0; - long *predy; - double total_perf, current_max = 0; - - struct Task *task = get_next_task(q); - struct MajModel *seed_model = msvmmaj_init_model(); - - clock_t main_s, main_e; - clock_t loop_s, loop_e; - - seed_model->m = task->train_data->m; - seed_model->K = task->train_data->K; - msvmmaj_allocate_model(seed_model); - msvmmaj_seed_model_V(NULL, seed_model, task->train_data); - - main_s = clock(); - while (task) { - total_perf = 0; - note("(%li/%li)\tw = %li\te = %f\tp = %f\tk = %f\tl = %f\t", - c+1, q->N, task->weight_idx, task->epsilon, - task->p, task->kappa, task->lambda); - loop_s = clock(); - struct MajModel *model = msvmmaj_init_model(); - make_model_from_task(task, model); - - model->n = task->train_data->n; - model->m = task->train_data->m; - model->K = task->train_data->K; - - msvmmaj_allocate_model(model); - msvmmaj_initialize_weights(task->train_data, model); - msvmmaj_seed_model_V(seed_model, model, task->train_data); - - fid = MSVMMAJ_OUTPUT_FILE; - MSVMMAJ_OUTPUT_FILE = NULL; - msvmmaj_optimize(model, task->train_data); - MSVMMAJ_OUTPUT_FILE = fid; - - predy = Calloc(long, task->test_data->n); - msvmmaj_predict_labels(task->test_data, task->train_data, - model, predy); - if (task->test_data->y != NULL) - total_perf = msvmmaj_prediction_perf(task->test_data, - predy); - msvmmaj_seed_model_V(model, seed_model, task->train_data); - - msvmmaj_free_model(model); - free(predy); - note("."); - loop_e = clock(); - current_max = maximum(current_max, total_perf); - note("\t%3.3f%% (%3.3fs)\t(best = %3.3f%%)\n", total_perf, - elapsed_time(loop_s, loop_e), current_max); - q->tasks[task->ID]->performance = total_perf; - task = get_next_task(q); - } - main_e = clock(); - - note("\nTotal elapsed time: %8.8f seconds\n", - elapsed_time(main_s, main_e)); - free(task); - msvmmaj_free_model(seed_model); -} - -/** - * @brief Free the Queue struct - * - * @details - * Freeing the allocated memory of the Queue means freeing every Task struct - * and then freeing the Queue. - * - * @param[in] q Queue to be freed - * - */ -void free_queue(struct Queue *q) -{ - long i; - for (i=0; iN; i++) { - free(q->tasks[i]->kernelparam); - free(q->tasks[i]); - } - free(q->tasks); - free(q); -} - -/** - * @brief Copy parameters from Task to MajModel - * - * @details - * A Task struct only contains the parameters of the MajModel to be estimated. - * This function is used to copy these parameters. - * - * @param[in] task Task instance with parameters - * @param[in,out] model MajModel to which the parameters are copied - */ -void make_model_from_task(struct Task *task, struct MajModel *model) -{ - // copy basic model parameters - model->weight_idx = task->weight_idx; - model->epsilon = task->epsilon; - model->p = task->p; - model->kappa = task->kappa; - model->lambda = task->lambda; - - // copy kernel parameters - model->kerneltype = task->kerneltype; - model->kernelparam = task->kernelparam; -} - -/** - * @brief Copy model parameters between two MajModel structs - * - * @details - * The parameters copied are MajModel::weight_idx, MajModel::epsilon, - * MajModel::p, MajModel::kappa, and MajModel::lambda. - * - * @param[in] from MajModel to copy parameters from - * @param[in,out] to MajModel to copy parameters to - */ -void copy_model(struct MajModel *from, struct MajModel *to) -{ - to->weight_idx = from->weight_idx; - to->epsilon = from->epsilon; - to->p = from->p; - to->kappa = from->kappa; - to->lambda = from->lambda; - - to->kerneltype = from->kerneltype; - switch (to->kerneltype) { - case K_LINEAR: - break; - case K_POLY: - to->kernelparam = Malloc(double, 3); - to->kernelparam[0] = from->kernelparam[0]; - to->kernelparam[1] = from->kernelparam[1]; - to->kernelparam[2] = from->kernelparam[2]; - break; - case K_RBF: - to->kernelparam = Malloc(double, 1); - to->kernelparam[0] = from->kernelparam[0]; - break; - case K_SIGMOID: - to->kernelparam = Malloc(double, 2); - to->kernelparam[0] = from->kernelparam[0]; - to->kernelparam[1] = from->kernelparam[1]; - break; - } -} - -/** - * @brief Print the description of the current task on screen - * - * @details - * To track the progress of the grid search the parameters of the current task - * are written to the output specified in MSVMMAJ_OUTPUT_FILE. Since the - * parameters differ with the specified kernel, this function writes a - * parameter string depending on which kernel is used. - * - * @param[in] task the Task specified - * @param[in] N total number of tasks - * - */ -void print_progress_string(struct Task *task, long N) -{ - char buffer[MAX_LINE_LENGTH]; - sprintf(buffer, "(%03li/%03li)\t", task->ID+1, N); - if (task->kerneltype == K_POLY) - sprintf(buffer + strlen(buffer), "d = %2.2f\t", - task->kernelparam[2]); - if (task->kerneltype == K_POLY || task->kerneltype == K_SIGMOID) - sprintf(buffer + strlen(buffer), "c = %2.2f\t", - task->kernelparam[1]); - if (task->kerneltype == K_POLY || task->kerneltype == K_SIGMOID || - task->kerneltype == K_RBF) - sprintf(buffer + strlen(buffer), "g = %3.3f\t", - task->kernelparam[0]); - sprintf(buffer + strlen(buffer), "eps = %g\tw = %i\tk = %2.2f\t" - "l = %f\tp = %2.2f\t", task->epsilon, - task->weight_idx, task->kappa, task->lambda, task->p); - note(buffer); -} diff --git a/src/predGenSVM.c b/src/predGenSVM.c new file mode 100644 index 0000000..7fac2ef --- /dev/null +++ b/src/predGenSVM.c @@ -0,0 +1,180 @@ +/* + * 20140317: + * THIS FUNCTION IS DEPRECATED, SINCE IT DOES NOT WORK WITH KERNELS. + * + */ + +/** + * @file predGenSVM.c + * @author Gertjan van den Burg + * @date January, 2014 + * @brief Command line interface for predicting class labels + * + * @details + * This is a command line program for predicting the class labels or + * determining the predictive performance of a pre-determined model on a given + * test dataset. The predictive performance can be written to the screen or + * the predicted class labels can be written to a specified output file. This + * is done using gensvm_write_predictions(). + * + * The specified model file must follow the specification given in + * gensvm_write_model(). + * + * For usage information, see the program help function. + * + */ + +#include "gensvm.h" +#include "gensvm_init.h" +#include "gensvm_io.h" +#include "gensvm_pred.h" +#include "util.h" + +#define MINARGS 3 + +extern FILE *GENSVM_OUTPUT_FILE; + +// function declarations +void exit_with_help(); +void parse_command_line(int argc, char **argv, + char *input_filename, char *output_filename, + char *model_filename); + +/** + * @brief Help function + */ +void exit_with_help() +{ + printf("This is GenSVM, version %1.1f\n\n", VERSION); + printf("Usage: predGenSVM [options] test_data_file model_file\n"); + printf("Options:\n"); + printf("-o output_file : write output to file\n"); + printf("-q : quiet mode (no output)\n"); + exit(0); +} + +/** + * @brief Main interface function for predGenSVM + * + * @details + * Main interface for the command line program. A given model file is read and + * a test dataset is initialized from the given data. The predictive + * performance (hitrate) of the model on the test set is printed to the output + * stream (default = stdout). If an output file is specified the predictions + * are written to the file. + * + * @todo + * Ensure that the program can read model files without class labels + * specified. In that case no prediction accuracy is printed to the screen. + * + * @param[in] argc number of command line arguments + * @param[in] argv array of command line arguments + * + */ +int main(int argc, char **argv) +{ + long *predy; + double performance; + + char input_filename[MAX_LINE_LENGTH]; + char model_filename[MAX_LINE_LENGTH]; + char output_filename[MAX_LINE_LENGTH];; + + if (argc < MINARGS || gensvm_check_argv(argc, argv, "-help") + || gensvm_check_argv_eq(argc, argv, "-h") ) + exit_with_help(); + parse_command_line(argc, argv, input_filename, output_filename, + model_filename); + + // read the data and model + struct GenModel *model = gensvm_init_model(); + struct GenData *data = gensvm_init_data(); + gensvm_read_data(data, input_filename); + gensvm_read_model(model, model_filename); + + // check if the number of attributes in data equals that in model + if (data->m != model->m) { + fprintf(stderr, "Error: number of attributes in data (%li) " + "does not equal the number of attributes in " + "model (%li)\n", data->m, model->m); + exit(1); + } else if (data->K != model->K) { + fprintf(stderr, "Error: number of classes in data (%li) " + "does not equal the number of classes in " + "model (%li)\n", data->K, model->K); + exit(1); + } + + // predict labels and performance if test data has labels + predy = Calloc(long, data->n); + gensvm_predict_labels(data, model, predy); + if (data->y != NULL) { + performance = gensvm_prediction_perf(data, predy); + note("Predictive performance: %3.2f%%\n", performance); + } + + // if output file is specified, write predictions to it + if (gensvm_check_argv_eq(argc, argv, "-o")) { + gensvm_write_predictions(data, predy, output_filename); + note("Predictions written to: %s\n", output_filename); + } + + // free the model, data, and predictions + gensvm_free_model(model); + gensvm_free_data(data); + free(predy); + + return 0; +} + +/** + * @brief Parse command line arguments + * + * @details + * Read the data filename and model filename from the command line arguments. + * If specified, also read the output filename. If the quiet flag is given, + * set the global output stream to NULL. On error, exit_with_help(). + * + * @param[in] argc number of command line arguments + * @param[in] argv array of command line arguments + * @param[in] input_filename pre-allocated array for the input + * filename + * @param[in] output_filename pre-allocated array for the output + * filename + * @param[in] model_filename pre-allocated array for the model + * filename + * + */ +void parse_command_line(int argc, char **argv, char *input_filename, + char *output_filename, char *model_filename) +{ + int i; + + GENSVM_OUTPUT_FILE = stdout; + + for (i=1; i= argc) + exit_with_help(); + switch (argv[i-1][1]) { + case 'o': + strcpy(output_filename, argv[i]); + break; + case 'q': + GENSVM_OUTPUT_FILE = NULL; + i--; + break; + default: + fprintf(stderr, "Unknown option: -%c\n", + argv[i-1][1]); + exit_with_help(); + } + } + + if (i >= argc) + exit_with_help(); + + strcpy(input_filename, argv[i]); + i++; + strcpy(model_filename, argv[i]); +} diff --git a/src/predMSVMMaj.c b/src/predMSVMMaj.c deleted file mode 100644 index 3dfcf08..0000000 --- a/src/predMSVMMaj.c +++ /dev/null @@ -1,180 +0,0 @@ -/* - * 20140317: - * THIS FUNCTION IS DEPRECATED, SINCE IT DOES NOT WORK WITH KERNELS. - * - */ - -/** - * @file predMSVMMaj.c - * @author Gertjan van den Burg - * @date January, 2014 - * @brief Command line interface for predicting class labels - * - * @details - * This is a command line program for predicting the class labels or - * determining the predictive performance of a pre-determined model on a given - * test dataset. The predictive performance can be written to the screen or - * the predicted class labels can be written to a specified output file. This - * is done using msvmmaj_write_predictions(). - * - * The specified model file must follow the specification given in - * msvmmaj_write_model(). - * - * For usage information, see the program help function. - * - */ - -#include "msvmmaj.h" -#include "msvmmaj_init.h" -#include "msvmmaj_io.h" -#include "msvmmaj_pred.h" -#include "util.h" - -#define MINARGS 3 - -extern FILE *MSVMMAJ_OUTPUT_FILE; - -// function declarations -void exit_with_help(); -void parse_command_line(int argc, char **argv, - char *input_filename, char *output_filename, - char *model_filename); - -/** - * @brief Help function - */ -void exit_with_help() -{ - printf("This is MSVMMaj, version %1.1f\n\n", VERSION); - printf("Usage: predMSVMMaj [options] test_data_file model_file\n"); - printf("Options:\n"); - printf("-o output_file : write output to file\n"); - printf("-q : quiet mode (no output)\n"); - exit(0); -} - -/** - * @brief Main interface function for predMSVMMaj - * - * @details - * Main interface for the command line program. A given model file is read and - * a test dataset is initialized from the given data. The predictive - * performance (hitrate) of the model on the test set is printed to the output - * stream (default = stdout). If an output file is specified the predictions - * are written to the file. - * - * @todo - * Ensure that the program can read model files without class labels - * specified. In that case no prediction accuracy is printed to the screen. - * - * @param[in] argc number of command line arguments - * @param[in] argv array of command line arguments - * - */ -int main(int argc, char **argv) -{ - long *predy; - double performance; - - char input_filename[MAX_LINE_LENGTH]; - char model_filename[MAX_LINE_LENGTH]; - char output_filename[MAX_LINE_LENGTH];; - - if (argc < MINARGS || msvmmaj_check_argv(argc, argv, "-help") - || msvmmaj_check_argv_eq(argc, argv, "-h") ) - exit_with_help(); - parse_command_line(argc, argv, input_filename, output_filename, - model_filename); - - // read the data and model - struct MajModel *model = msvmmaj_init_model(); - struct MajData *data = msvmmaj_init_data(); - msvmmaj_read_data(data, input_filename); - msvmmaj_read_model(model, model_filename); - - // check if the number of attributes in data equals that in model - if (data->m != model->m) { - fprintf(stderr, "Error: number of attributes in data (%li) " - "does not equal the number of attributes in " - "model (%li)\n", data->m, model->m); - exit(1); - } else if (data->K != model->K) { - fprintf(stderr, "Error: number of classes in data (%li) " - "does not equal the number of classes in " - "model (%li)\n", data->K, model->K); - exit(1); - } - - // predict labels and performance if test data has labels - predy = Calloc(long, data->n); - msvmmaj_predict_labels(data, model, predy); - if (data->y != NULL) { - performance = msvmmaj_prediction_perf(data, predy); - note("Predictive performance: %3.2f%%\n", performance); - } - - // if output file is specified, write predictions to it - if (msvmmaj_check_argv_eq(argc, argv, "-o")) { - msvmmaj_write_predictions(data, predy, output_filename); - note("Predictions written to: %s\n", output_filename); - } - - // free the model, data, and predictions - msvmmaj_free_model(model); - msvmmaj_free_data(data); - free(predy); - - return 0; -} - -/** - * @brief Parse command line arguments - * - * @details - * Read the data filename and model filename from the command line arguments. - * If specified, also read the output filename. If the quiet flag is given, - * set the global output stream to NULL. On error, exit_with_help(). - * - * @param[in] argc number of command line arguments - * @param[in] argv array of command line arguments - * @param[in] input_filename pre-allocated array for the input - * filename - * @param[in] output_filename pre-allocated array for the output - * filename - * @param[in] model_filename pre-allocated array for the model - * filename - * - */ -void parse_command_line(int argc, char **argv, char *input_filename, - char *output_filename, char *model_filename) -{ - int i; - - MSVMMAJ_OUTPUT_FILE = stdout; - - for (i=1; i= argc) - exit_with_help(); - switch (argv[i-1][1]) { - case 'o': - strcpy(output_filename, argv[i]); - break; - case 'q': - MSVMMAJ_OUTPUT_FILE = NULL; - i--; - break; - default: - fprintf(stderr, "Unknown option: -%c\n", - argv[i-1][1]); - exit_with_help(); - } - } - - if (i >= argc) - exit_with_help(); - - strcpy(input_filename, argv[i]); - i++; - strcpy(model_filename, argv[i]); -} diff --git a/src/trainGenSVM.c b/src/trainGenSVM.c new file mode 100644 index 0000000..eb75f5d --- /dev/null +++ b/src/trainGenSVM.c @@ -0,0 +1,244 @@ +/** + * @file trainGenSVM.c + * @author Gertjan van den Burg + * @date August, 2013 + * @brief Command line interface for training a single model with GenSVM + * + * @details + * This is a command line program for training a single model on a given + * dataset. To run a grid search over a number of parameter configurations, + * see trainGenSVMdataset.c. + * + */ + +#include +#include + +#include "gensvm_kernel.h" +#include "libGenSVM.h" +#include "gensvm.h" +#include "gensvm_io.h" +#include "gensvm_init.h" +#include "gensvm_train.h" +#include "util.h" + +#define MINARGS 2 + +extern FILE *GENSVM_OUTPUT_FILE; + +// function declarations +void exit_with_help(); +void parse_command_line(int argc, char **argv, struct GenModel *model, + char *input_filename, char *output_filename, char *model_filename); + +/** + * @brief Help function + */ +void exit_with_help() +{ + printf("This is GenSVM, version %1.1f\n\n", VERSION); + printf("Usage: trainGenSVM [options] training_data_file\n"); + printf("Options:\n"); + printf("-c coef : coefficient for the polynomial and sigmoid kernel\n"); + printf("-d degree : degree for the polynomial kernel\n"); + printf("-e epsilon : set the value of the stopping criterion\n"); + printf("-g gamma : parameter for the rbf, polynomial or sigmoid " + "kernel\n"); + printf("-h | -help : print this help.\n"); + printf("-k kappa : set the value of kappa used in the Huber hinge\n"); + printf("-l lambda : set the value of lambda (lambda > 0)\n"); + printf("-m model_file : use previous model as seed for W and t\n"); + printf("-o output_file : write output to file\n"); + printf("-p p-value : set the value of p in the lp norm " + "(1.0 <= p <= 2.0)\n"); + printf("-q : quiet mode (no output)\n"); + printf("-r rho : choose the weigth specification (1 = unit, 2 = " + "group)\n"); + printf("-t type: kerneltype (LINEAR=0, POLY=1, RBF=2, SIGMOID=3)\n"); + + exit(0); +} + +/** + * @brief Main interface function for trainGenSVM + * + * @details + * Main interface for the command line program. A given dataset file is read + * and a GenSVM model is trained on this data. By default the progress of the + * computations are written to stdout. See for full options of the program the + * help function. + * + * @param[in] argc number of command line arguments + * @param[in] argv array of command line arguments + * + */ +int main(int argc, char **argv) +{ + char input_filename[MAX_LINE_LENGTH]; + char model_filename[MAX_LINE_LENGTH]; + char output_filename[MAX_LINE_LENGTH]; + + struct GenModel *model = gensvm_init_model(); + struct GenData *data = gensvm_init_data(); + + if (argc < MINARGS || gensvm_check_argv(argc, argv, "-help") + || gensvm_check_argv_eq(argc, argv, "-h") ) + exit_with_help(); + parse_command_line(argc, argv, model, input_filename, + output_filename, model_filename); + + // read data file + gensvm_read_data(data, input_filename); + + // copy dataset parameters to model + model->n = data->n; + model->m = data->m; + model->K = data->K; + model->data_file = input_filename; + + // allocate model + gensvm_allocate_model(model); + + // initialize kernel (if necessary) + gensvm_make_kernel(model, data); + + // reallocate model and initialize weights + gensvm_reallocate_model(model, data->n, data->m); + gensvm_initialize_weights(data, model); + + // seed the random number generator (only place in programs is in + // command line interfaces) + srand(time(NULL)); + + if (gensvm_check_argv_eq(argc, argv, "-m")) { + struct GenModel *seed_model = gensvm_init_model(); + gensvm_read_model(seed_model, model_filename); + gensvm_seed_model_V(seed_model, model, data); + gensvm_free_model(seed_model); + } else { + gensvm_seed_model_V(NULL, model, data); + } + + // start training + gensvm_optimize(model, data); + + // write_model to file + if (gensvm_check_argv_eq(argc, argv, "-o")) { + gensvm_write_model(model, output_filename); + note("Output written to %s\n", output_filename); + } + + // free model and data + gensvm_free_model(model); + gensvm_free_data(data); + + return 0; +} + +/** + * @brief Parse command line arguments + * + * @details + * Process the command line arguments for the model parameters, and record + * them in the specified GenModel. An input filename for the dataset is read + * and if specified an output filename and a model filename for the seed + * model. + * + * @param[in] argc number of command line arguments + * @param[in] argv array of command line arguments + * @param[in] model initialized model + * @param[in] input_filename pre-allocated buffer for the input + * filename + * @param[in] output_filename pre-allocated buffer for the output + * filename + * @param[in] model_filename pre-allocated buffer for the model + * filename + * + */ +void parse_command_line(int argc, char **argv, struct GenModel *model, + char *input_filename, char *output_filename, char *model_filename) +{ + int i; + double gamma = 1.0, + degree = 2.0, + coef = 0.0; + + GENSVM_OUTPUT_FILE = stdout; + + // parse options + for (i=1; i=argc) { + exit_with_help(); + } + switch (argv[i-1][1]) { + case 'c': + coef = atof(argv[i]); + break; + case 'd': + degree = atof(argv[i]); + break; + case 'e': + model->epsilon = atof(argv[i]); + break; + case 'g': + gamma = atof(argv[i]); + break; + case 'k': + model->kappa = atof(argv[i]); + break; + case 'l': + model->lambda = atof(argv[i]); + break; + case 'm': + strcpy(model_filename, argv[i]); + break; + case 'o': + strcpy(output_filename, argv[i]); + break; + case 'p': + model->p = atof(argv[i]); + break; + case 'r': + model->weight_idx = atoi(argv[i]); + break; + case 't': + model->kerneltype = atoi(argv[i]); + break; + case 'q': + GENSVM_OUTPUT_FILE = NULL; + i--; + break; + default: + fprintf(stderr, "Unknown option: -%c\n", + argv[i-1][1]); + exit_with_help(); + } + } + + // read input filename + if (i >= argc) + exit_with_help(); + + strcpy(input_filename, argv[i]); + + // set kernel parameters + switch (model->kerneltype) { + case K_LINEAR: + break; + case K_POLY: + model->kernelparam = Calloc(double, 3); + model->kernelparam[0] = gamma; + model->kernelparam[1] = coef; + model->kernelparam[2] = degree; + break; + case K_RBF: + model->kernelparam = Calloc(double, 1); + model->kernelparam[0] = gamma; + break; + case K_SIGMOID: + model->kernelparam = Calloc(double, 1); + model->kernelparam[0] = gamma; + model->kernelparam[1] = coef; + } +} diff --git a/src/trainGenSVMdataset.c b/src/trainGenSVMdataset.c new file mode 100644 index 0000000..2882c8f --- /dev/null +++ b/src/trainGenSVMdataset.c @@ -0,0 +1,321 @@ +/** + * @file trainGenSVMdataset.c + * @author Gertjan van den Burg + * @date January, 2014 + * @brief Command line interface for the grid search program + * + * @details + * This is a command line interface to the parameter grid search functionality + * of the algorithm. The grid search is specified in a separate file, thereby + * reducing the number of command line arguments. See + * read_training_from_file() for documentation on the training file. + * + * The program runs a grid search as specified in the training file. If + * desired the grid search can incorporate consistency checks to find the + * configuration among the best configurations which scores consistently high. + * All output is written to stdout, unless the quiet mode is specified. + * + * For further usage information, see the program help function. + * + */ + +#include + +#include "crossval.h" +#include "gensvm.h" +#include "gensvm_io.h" +#include "gensvm_init.h" +#include "gensvm_pred.h" +#include "gensvm_train.h" +#include "gensvm_train_dataset.h" +#include "strutil.h" +#include "util.h" + +#define MINARGS 2 + +extern FILE *GENSVM_OUTPUT_FILE; + +// function declarations +void exit_with_help(); +void parse_command_line(int argc, char **argv, char *input_filename); +void read_training_from_file(char *input_filename, struct Training *training); + +/** + * @brief Help function + */ +void exit_with_help() +{ + printf("This is GenSVM, version %1.1f\n\n", VERSION); + printf("Usage: trainGenSVMdataset [options] training_file\n"); + printf("Options:\n"); + printf("-h | -help : print this help.\n"); + printf("-q : quiet mode (no output)\n"); + + exit(0); +} + +/** + * @brief Main interface function for trainGenSVMdataset + * + * @details + * Main interface for the command line program. A given training file which + * specifies a grid search over a single dataset is read. From this, a Queue + * is created containing all Task instances that need to be performed in the + * search. Depending on the type of dataset, either cross validation or + * train/test split training is performed for all tasks. If specified, + * consistency repeats are done at the end of the grid search. Note that + * currently no output is produced other than what is written to stdout. + * + * @param[in] argc number of command line arguments + * @param[in] argv array of command line arguments + * + */ +int main(int argc, char **argv) +{ + char input_filename[MAX_LINE_LENGTH]; + + struct Training *training = Malloc(struct Training, 1); + struct GenData *train_data = Malloc(struct GenData, 1); + struct GenData *test_data = Malloc(struct GenData, 1); + + if (argc < MINARGS || gensvm_check_argv(argc, argv, "-help") + || gensvm_check_argv_eq(argc, argv, "-h") ) + exit_with_help(); + parse_command_line(argc, argv, input_filename); + + training->repeats = 0; + note("Reading training file\n"); + read_training_from_file(input_filename, training); + + note("Reading data from %s\n", training->train_data_file); + gensvm_read_data(train_data, training->train_data_file); + if (training->traintype == TT) { + note("Reading data from %s\n", training->test_data_file); + gensvm_read_data(test_data, training->test_data_file); + } + + note("Creating queue\n"); + struct Queue *q = Malloc(struct Queue, 1); + make_queue(training, q, train_data, test_data); + + srand(time(NULL)); + + note("Starting training\n"); + if (training->traintype == TT) + start_training_tt(q); + else + start_training_cv(q); + note("Training finished\n"); + + if (training->repeats > 0) { + consistency_repeats(q, training->repeats, training->traintype); + } + + free_queue(q); + free(training); + gensvm_free_data(train_data); + gensvm_free_data(test_data); + + note("Done.\n"); + return 0; +} + +/** + * @brief Parse command line arguments + * + * @details + * Few arguments can be supplied to the command line. Only quiet mode can be + * specified, or help can be requested. The filename of the training file is + * read from the arguments. Parsing of the training file is done separately in + * read_training_from_file(). + * + * @param[in] argc number of command line arguments + * @param[in] argv array of command line arguments + * @param[in] input_filename pre-allocated buffer for the training + * filename. + * + */ +void parse_command_line(int argc, char **argv, char *input_filename) +{ + int i; + + GENSVM_OUTPUT_FILE = stdout; + + for (i=1; i=argc) + exit_with_help(); + switch (argv[i-1][1]) { + case 'q': + GENSVM_OUTPUT_FILE = NULL; + i--; + break; + default: + fprintf(stderr, "Unknown option: -%c\n", + argv[i-1][1]); + exit_with_help(); + } + } + + if (i >= argc) + exit_with_help(); + + strcpy(input_filename, argv[i]); +} + +KernelType parse_kernel_str(char *kernel_line) +{ + if (str_endswith(kernel_line, "LINEAR\n")) { + return K_LINEAR; + } else if (str_endswith(kernel_line, "POLY\n")) { + return K_POLY; + } else if (str_endswith(kernel_line, "RBF\n")) { + return K_RBF; + } else if (str_endswith(kernel_line, "SIGMOID\n")) { + return K_SIGMOID; + } else { + fprintf(stderr, "Unknown kernel specified on line: %s\n", + kernel_line); + exit(1); + } +} + +/** + * @brief Read the Training struct from file + * + * @details + * Read the Training struct from a file. The training file follows a specific + * format specified in @ref spec_training_file. + * + * Commonly used string functions in this function are all_doubles_str() and + * all_longs_str(). + * + * @param[in] input_filename filename of the training file + * @param[in] training Training structure to place the parsed + * parameter grid. + * + */ +void read_training_from_file(char *input_filename, struct Training *training) +{ + long i, nr = 0; + FILE *fid; + char buffer[MAX_LINE_LENGTH]; + char train_filename[MAX_LINE_LENGTH]; + char test_filename[MAX_LINE_LENGTH]; + double *params = Calloc(double, MAX_LINE_LENGTH); + long *lparams = Calloc(long, MAX_LINE_LENGTH); + + fid = fopen(input_filename, "r"); + if (fid == NULL) { + fprintf(stderr, "Error opening training file %s\n", + input_filename); + exit(1); + } + training->traintype = CV; + while ( fgets(buffer, MAX_LINE_LENGTH, fid) != NULL ) { + Memset(params, double, MAX_LINE_LENGTH); + Memset(lparams, long, MAX_LINE_LENGTH); + if (str_startswith(buffer, "train:")) { + sscanf(buffer, "train: %s\n", train_filename); + training->train_data_file = Calloc(char, + MAX_LINE_LENGTH); + strcpy(training->train_data_file, train_filename); + } else if (str_startswith(buffer, "test:")) { + sscanf(buffer, "test: %s\n", test_filename); + training->test_data_file = Calloc(char, + MAX_LINE_LENGTH); + strcpy(training->test_data_file, test_filename); + training->traintype = TT; + } else if (str_startswith(buffer, "p:")) { + nr = all_doubles_str(buffer, 2, params); + training->ps = Calloc(double, nr); + for (i=0; ips[i] = params[i]; + training->Np = nr; + } else if (str_startswith(buffer, "lambda:")) { + nr = all_doubles_str(buffer, 7, params); + training->lambdas = Calloc(double, nr); + for (i=0; ilambdas[i] = params[i]; + training->Nl = nr; + } else if (str_startswith(buffer, "kappa:")) { + nr = all_doubles_str(buffer, 6, params); + training->kappas = Calloc(double, nr); + for (i=0; ikappas[i] = params[i]; + training->Nk = nr; + } else if (str_startswith(buffer, "epsilon:")) { + nr = all_doubles_str(buffer, 8, params); + training->epsilons = Calloc(double, nr); + for (i=0; iepsilons[i] = params[i]; + training->Ne = nr; + } else if (str_startswith(buffer, "weight:")) { + nr = all_longs_str(buffer, 7, lparams); + training->weight_idxs = Calloc(int, nr); + for (i=0; iweight_idxs[i] = lparams[i]; + training->Nw = nr; + } else if (str_startswith(buffer, "folds:")) { + nr = all_longs_str(buffer, 6, lparams); + training->folds = lparams[0]; + if (nr > 1) + fprintf(stderr, "Field \"folds\" only takes " + "one value. Additional " + "fields are ignored.\n"); + } else if (str_startswith(buffer, "repeats:")) { + nr = all_longs_str(buffer, 8, lparams); + training->repeats = lparams[0]; + if (nr > 1) + fprintf(stderr, "Field \"repeats\" only " + "takes one value. Additional " + "fields are ignored.\n"); + } else if (str_startswith(buffer, "kernel:")) { + training->kerneltype = parse_kernel_str(buffer); + } else if (str_startswith(buffer, "gamma:")) { + nr = all_doubles_str(buffer, 6, params); + if (training->kerneltype == K_LINEAR) { + fprintf(stderr, "Field \"gamma\" ignored, " + "linear kernel is used.\n"); + training->Ng = 0; + break; + } + training->gammas = Calloc(double, nr); + for (i=0; igammas[i] = params[i]; + training->Ng = nr; + } else if (str_startswith(buffer, "coef:")) { + nr = all_doubles_str(buffer, 5, params); + if (training->kerneltype == K_LINEAR || + training->kerneltype == K_RBF) { + fprintf(stderr, "Field \"coef\" ignored with " + "specified kernel.\n"); + training->Nc = 0; + break; + } + training->coefs = Calloc(double, nr); + for (i=0; icoefs[i] = params[i]; + training->Nc = nr; + } else if (str_startswith(buffer, "degree:")) { + nr = all_doubles_str(buffer, 7, params); + if (training->kerneltype != K_POLY) { + fprintf(stderr, "Field \"degree\" ignored " + "with specified kernel.\n"); + training->Nd = 0; + break; + } + training->degrees = Calloc(double, nr); + for (i=0; idegrees[i] = params[i]; + training->Nd = nr; + } else { + fprintf(stderr, "Cannot find any parameters on line: " + "%s\n", buffer); + } + } + + free(params); + free(lparams); + fclose(fid); +} diff --git a/src/trainMSVMMaj.c b/src/trainMSVMMaj.c deleted file mode 100644 index 5377b43..0000000 --- a/src/trainMSVMMaj.c +++ /dev/null @@ -1,244 +0,0 @@ -/** - * @file trainMSVMMaj.c - * @author Gertjan van den Burg - * @date August, 2013 - * @brief Command line interface for training a single model with MSVMMaj - * - * @details - * This is a command line program for training a single model on a given - * dataset. To run a grid search over a number of parameter configurations, - * see trainMSVMMajdataset.c. - * - */ - -#include -#include - -#include "msvmmaj_kernel.h" -#include "libMSVMMaj.h" -#include "msvmmaj.h" -#include "msvmmaj_io.h" -#include "msvmmaj_init.h" -#include "msvmmaj_train.h" -#include "util.h" - -#define MINARGS 2 - -extern FILE *MSVMMAJ_OUTPUT_FILE; - -// function declarations -void exit_with_help(); -void parse_command_line(int argc, char **argv, struct MajModel *model, - char *input_filename, char *output_filename, char *model_filename); - -/** - * @brief Help function - */ -void exit_with_help() -{ - printf("This is MSVMMaj, version %1.1f\n\n", VERSION); - printf("Usage: trainMSVMMaj [options] training_data_file\n"); - printf("Options:\n"); - printf("-c coef : coefficient for the polynomial and sigmoid kernel\n"); - printf("-d degree : degree for the polynomial kernel\n"); - printf("-e epsilon : set the value of the stopping criterion\n"); - printf("-g gamma : parameter for the rbf, polynomial or sigmoid " - "kernel\n"); - printf("-h | -help : print this help.\n"); - printf("-k kappa : set the value of kappa used in the Huber hinge\n"); - printf("-l lambda : set the value of lambda (lambda > 0)\n"); - printf("-m model_file : use previous model as seed for W and t\n"); - printf("-o output_file : write output to file\n"); - printf("-p p-value : set the value of p in the lp norm " - "(1.0 <= p <= 2.0)\n"); - printf("-q : quiet mode (no output)\n"); - printf("-r rho : choose the weigth specification (1 = unit, 2 = " - "group)\n"); - printf("-t type: kerneltype (LINEAR=0, POLY=1, RBF=2, SIGMOID=3)\n"); - - exit(0); -} - -/** - * @brief Main interface function for trainMSVMMaj - * - * @details - * Main interface for the command line program. A given dataset file is read - * and a MSVMMaj model is trained on this data. By default the progress of the - * computations are written to stdout. See for full options of the program the - * help function. - * - * @param[in] argc number of command line arguments - * @param[in] argv array of command line arguments - * - */ -int main(int argc, char **argv) -{ - char input_filename[MAX_LINE_LENGTH]; - char model_filename[MAX_LINE_LENGTH]; - char output_filename[MAX_LINE_LENGTH]; - - struct MajModel *model = msvmmaj_init_model(); - struct MajData *data = msvmmaj_init_data(); - - if (argc < MINARGS || msvmmaj_check_argv(argc, argv, "-help") - || msvmmaj_check_argv_eq(argc, argv, "-h") ) - exit_with_help(); - parse_command_line(argc, argv, model, input_filename, - output_filename, model_filename); - - // read data file - msvmmaj_read_data(data, input_filename); - - // copy dataset parameters to model - model->n = data->n; - model->m = data->m; - model->K = data->K; - model->data_file = input_filename; - - // allocate model - msvmmaj_allocate_model(model); - - // initialize kernel (if necessary) - msvmmaj_make_kernel(model, data); - - // reallocate model and initialize weights - msvmmaj_reallocate_model(model, data->n, data->m); - msvmmaj_initialize_weights(data, model); - - // seed the random number generator (only place in programs is in - // command line interfaces) - srand(time(NULL)); - - if (msvmmaj_check_argv_eq(argc, argv, "-m")) { - struct MajModel *seed_model = msvmmaj_init_model(); - msvmmaj_read_model(seed_model, model_filename); - msvmmaj_seed_model_V(seed_model, model, data); - msvmmaj_free_model(seed_model); - } else { - msvmmaj_seed_model_V(NULL, model, data); - } - - // start training - msvmmaj_optimize(model, data); - - // write_model to file - if (msvmmaj_check_argv_eq(argc, argv, "-o")) { - msvmmaj_write_model(model, output_filename); - note("Output written to %s\n", output_filename); - } - - // free model and data - msvmmaj_free_model(model); - msvmmaj_free_data(data); - - return 0; -} - -/** - * @brief Parse command line arguments - * - * @details - * Process the command line arguments for the model parameters, and record - * them in the specified MajModel. An input filename for the dataset is read - * and if specified an output filename and a model filename for the seed - * model. - * - * @param[in] argc number of command line arguments - * @param[in] argv array of command line arguments - * @param[in] model initialized model - * @param[in] input_filename pre-allocated buffer for the input - * filename - * @param[in] output_filename pre-allocated buffer for the output - * filename - * @param[in] model_filename pre-allocated buffer for the model - * filename - * - */ -void parse_command_line(int argc, char **argv, struct MajModel *model, - char *input_filename, char *output_filename, char *model_filename) -{ - int i; - double gamma = 1.0, - degree = 2.0, - coef = 0.0; - - MSVMMAJ_OUTPUT_FILE = stdout; - - // parse options - for (i=1; i=argc) { - exit_with_help(); - } - switch (argv[i-1][1]) { - case 'c': - coef = atof(argv[i]); - break; - case 'd': - degree = atof(argv[i]); - break; - case 'e': - model->epsilon = atof(argv[i]); - break; - case 'g': - gamma = atof(argv[i]); - break; - case 'k': - model->kappa = atof(argv[i]); - break; - case 'l': - model->lambda = atof(argv[i]); - break; - case 'm': - strcpy(model_filename, argv[i]); - break; - case 'o': - strcpy(output_filename, argv[i]); - break; - case 'p': - model->p = atof(argv[i]); - break; - case 'r': - model->weight_idx = atoi(argv[i]); - break; - case 't': - model->kerneltype = atoi(argv[i]); - break; - case 'q': - MSVMMAJ_OUTPUT_FILE = NULL; - i--; - break; - default: - fprintf(stderr, "Unknown option: -%c\n", - argv[i-1][1]); - exit_with_help(); - } - } - - // read input filename - if (i >= argc) - exit_with_help(); - - strcpy(input_filename, argv[i]); - - // set kernel parameters - switch (model->kerneltype) { - case K_LINEAR: - break; - case K_POLY: - model->kernelparam = Calloc(double, 3); - model->kernelparam[0] = gamma; - model->kernelparam[1] = coef; - model->kernelparam[2] = degree; - break; - case K_RBF: - model->kernelparam = Calloc(double, 1); - model->kernelparam[0] = gamma; - break; - case K_SIGMOID: - model->kernelparam = Calloc(double, 1); - model->kernelparam[0] = gamma; - model->kernelparam[1] = coef; - } -} diff --git a/src/trainMSVMMajdataset.c b/src/trainMSVMMajdataset.c deleted file mode 100644 index b9d9180..0000000 --- a/src/trainMSVMMajdataset.c +++ /dev/null @@ -1,321 +0,0 @@ -/** - * @file trainMSVMMajdataset.c - * @author Gertjan van den Burg - * @date January, 2014 - * @brief Command line interface for the grid search program - * - * @details - * This is a command line interface to the parameter grid search functionality - * of the algorithm. The grid search is specified in a separate file, thereby - * reducing the number of command line arguments. See - * read_training_from_file() for documentation on the training file. - * - * The program runs a grid search as specified in the training file. If - * desired the grid search can incorporate consistency checks to find the - * configuration among the best configurations which scores consistently high. - * All output is written to stdout, unless the quiet mode is specified. - * - * For further usage information, see the program help function. - * - */ - -#include - -#include "crossval.h" -#include "msvmmaj.h" -#include "msvmmaj_io.h" -#include "msvmmaj_init.h" -#include "msvmmaj_pred.h" -#include "msvmmaj_train.h" -#include "msvmmaj_train_dataset.h" -#include "strutil.h" -#include "util.h" - -#define MINARGS 2 - -extern FILE *MSVMMAJ_OUTPUT_FILE; - -// function declarations -void exit_with_help(); -void parse_command_line(int argc, char **argv, char *input_filename); -void read_training_from_file(char *input_filename, struct Training *training); - -/** - * @brief Help function - */ -void exit_with_help() -{ - printf("This is MSVMMaj, version %1.1f\n\n", VERSION); - printf("Usage: trainMSVMMajdataset [options] training_file\n"); - printf("Options:\n"); - printf("-h | -help : print this help.\n"); - printf("-q : quiet mode (no output)\n"); - - exit(0); -} - -/** - * @brief Main interface function for trainMSVMMajdataset - * - * @details - * Main interface for the command line program. A given training file which - * specifies a grid search over a single dataset is read. From this, a Queue - * is created containing all Task instances that need to be performed in the - * search. Depending on the type of dataset, either cross validation or - * train/test split training is performed for all tasks. If specified, - * consistency repeats are done at the end of the grid search. Note that - * currently no output is produced other than what is written to stdout. - * - * @param[in] argc number of command line arguments - * @param[in] argv array of command line arguments - * - */ -int main(int argc, char **argv) -{ - char input_filename[MAX_LINE_LENGTH]; - - struct Training *training = Malloc(struct Training, 1); - struct MajData *train_data = Malloc(struct MajData, 1); - struct MajData *test_data = Malloc(struct MajData, 1); - - if (argc < MINARGS || msvmmaj_check_argv(argc, argv, "-help") - || msvmmaj_check_argv_eq(argc, argv, "-h") ) - exit_with_help(); - parse_command_line(argc, argv, input_filename); - - training->repeats = 0; - note("Reading training file\n"); - read_training_from_file(input_filename, training); - - note("Reading data from %s\n", training->train_data_file); - msvmmaj_read_data(train_data, training->train_data_file); - if (training->traintype == TT) { - note("Reading data from %s\n", training->test_data_file); - msvmmaj_read_data(test_data, training->test_data_file); - } - - note("Creating queue\n"); - struct Queue *q = Malloc(struct Queue, 1); - make_queue(training, q, train_data, test_data); - - srand(time(NULL)); - - note("Starting training\n"); - if (training->traintype == TT) - start_training_tt(q); - else - start_training_cv(q); - note("Training finished\n"); - - if (training->repeats > 0) { - consistency_repeats(q, training->repeats, training->traintype); - } - - free_queue(q); - free(training); - msvmmaj_free_data(train_data); - msvmmaj_free_data(test_data); - - note("Done.\n"); - return 0; -} - -/** - * @brief Parse command line arguments - * - * @details - * Few arguments can be supplied to the command line. Only quiet mode can be - * specified, or help can be requested. The filename of the training file is - * read from the arguments. Parsing of the training file is done separately in - * read_training_from_file(). - * - * @param[in] argc number of command line arguments - * @param[in] argv array of command line arguments - * @param[in] input_filename pre-allocated buffer for the training - * filename. - * - */ -void parse_command_line(int argc, char **argv, char *input_filename) -{ - int i; - - MSVMMAJ_OUTPUT_FILE = stdout; - - for (i=1; i=argc) - exit_with_help(); - switch (argv[i-1][1]) { - case 'q': - MSVMMAJ_OUTPUT_FILE = NULL; - i--; - break; - default: - fprintf(stderr, "Unknown option: -%c\n", - argv[i-1][1]); - exit_with_help(); - } - } - - if (i >= argc) - exit_with_help(); - - strcpy(input_filename, argv[i]); -} - -KernelType parse_kernel_str(char *kernel_line) -{ - if (str_endswith(kernel_line, "LINEAR\n")) { - return K_LINEAR; - } else if (str_endswith(kernel_line, "POLY\n")) { - return K_POLY; - } else if (str_endswith(kernel_line, "RBF\n")) { - return K_RBF; - } else if (str_endswith(kernel_line, "SIGMOID\n")) { - return K_SIGMOID; - } else { - fprintf(stderr, "Unknown kernel specified on line: %s\n", - kernel_line); - exit(1); - } -} - -/** - * @brief Read the Training struct from file - * - * @details - * Read the Training struct from a file. The training file follows a specific - * format specified in @ref spec_training_file. - * - * Commonly used string functions in this function are all_doubles_str() and - * all_longs_str(). - * - * @param[in] input_filename filename of the training file - * @param[in] training Training structure to place the parsed - * parameter grid. - * - */ -void read_training_from_file(char *input_filename, struct Training *training) -{ - long i, nr = 0; - FILE *fid; - char buffer[MAX_LINE_LENGTH]; - char train_filename[MAX_LINE_LENGTH]; - char test_filename[MAX_LINE_LENGTH]; - double *params = Calloc(double, MAX_LINE_LENGTH); - long *lparams = Calloc(long, MAX_LINE_LENGTH); - - fid = fopen(input_filename, "r"); - if (fid == NULL) { - fprintf(stderr, "Error opening training file %s\n", - input_filename); - exit(1); - } - training->traintype = CV; - while ( fgets(buffer, MAX_LINE_LENGTH, fid) != NULL ) { - Memset(params, double, MAX_LINE_LENGTH); - Memset(lparams, long, MAX_LINE_LENGTH); - if (str_startswith(buffer, "train:")) { - sscanf(buffer, "train: %s\n", train_filename); - training->train_data_file = Calloc(char, - MAX_LINE_LENGTH); - strcpy(training->train_data_file, train_filename); - } else if (str_startswith(buffer, "test:")) { - sscanf(buffer, "test: %s\n", test_filename); - training->test_data_file = Calloc(char, - MAX_LINE_LENGTH); - strcpy(training->test_data_file, test_filename); - training->traintype = TT; - } else if (str_startswith(buffer, "p:")) { - nr = all_doubles_str(buffer, 2, params); - training->ps = Calloc(double, nr); - for (i=0; ips[i] = params[i]; - training->Np = nr; - } else if (str_startswith(buffer, "lambda:")) { - nr = all_doubles_str(buffer, 7, params); - training->lambdas = Calloc(double, nr); - for (i=0; ilambdas[i] = params[i]; - training->Nl = nr; - } else if (str_startswith(buffer, "kappa:")) { - nr = all_doubles_str(buffer, 6, params); - training->kappas = Calloc(double, nr); - for (i=0; ikappas[i] = params[i]; - training->Nk = nr; - } else if (str_startswith(buffer, "epsilon:")) { - nr = all_doubles_str(buffer, 8, params); - training->epsilons = Calloc(double, nr); - for (i=0; iepsilons[i] = params[i]; - training->Ne = nr; - } else if (str_startswith(buffer, "weight:")) { - nr = all_longs_str(buffer, 7, lparams); - training->weight_idxs = Calloc(int, nr); - for (i=0; iweight_idxs[i] = lparams[i]; - training->Nw = nr; - } else if (str_startswith(buffer, "folds:")) { - nr = all_longs_str(buffer, 6, lparams); - training->folds = lparams[0]; - if (nr > 1) - fprintf(stderr, "Field \"folds\" only takes " - "one value. Additional " - "fields are ignored.\n"); - } else if (str_startswith(buffer, "repeats:")) { - nr = all_longs_str(buffer, 8, lparams); - training->repeats = lparams[0]; - if (nr > 1) - fprintf(stderr, "Field \"repeats\" only " - "takes one value. Additional " - "fields are ignored.\n"); - } else if (str_startswith(buffer, "kernel:")) { - training->kerneltype = parse_kernel_str(buffer); - } else if (str_startswith(buffer, "gamma:")) { - nr = all_doubles_str(buffer, 6, params); - if (training->kerneltype == K_LINEAR) { - fprintf(stderr, "Field \"gamma\" ignored, " - "linear kernel is used.\n"); - training->Ng = 0; - break; - } - training->gammas = Calloc(double, nr); - for (i=0; igammas[i] = params[i]; - training->Ng = nr; - } else if (str_startswith(buffer, "coef:")) { - nr = all_doubles_str(buffer, 5, params); - if (training->kerneltype == K_LINEAR || - training->kerneltype == K_RBF) { - fprintf(stderr, "Field \"coef\" ignored with " - "specified kernel.\n"); - training->Nc = 0; - break; - } - training->coefs = Calloc(double, nr); - for (i=0; icoefs[i] = params[i]; - training->Nc = nr; - } else if (str_startswith(buffer, "degree:")) { - nr = all_doubles_str(buffer, 7, params); - if (training->kerneltype != K_POLY) { - fprintf(stderr, "Field \"degree\" ignored " - "with specified kernel.\n"); - training->Nd = 0; - break; - } - training->degrees = Calloc(double, nr); - for (i=0; idegrees[i] = params[i]; - training->Nd = nr; - } else { - fprintf(stderr, "Cannot find any parameters on line: " - "%s\n", buffer); - } - } - - free(params); - free(lparams); - fclose(fid); -} diff --git a/src/util.c b/src/util.c index e76a074..23ee4e5 100644 --- a/src/util.c +++ b/src/util.c @@ -13,11 +13,11 @@ #include "util.h" -FILE *MSVMMAJ_OUTPUT_FILE; ///< The #MSVMMAJ_OUTPUT_FILE specifies the +FILE *GENSVM_OUTPUT_FILE; ///< The #GENSVM_OUTPUT_FILE specifies the ///< output stream to which all output is ///< written. This is done through the ///< internal (!) - ///< function msvmmaj_print_string(). The + ///< function gensvm_print_string(). The ///< advantage of using a global output ///< stream variable is that the output can ///< temporarily be suppressed by importing @@ -40,7 +40,7 @@ FILE *MSVMMAJ_OUTPUT_FILE; ///< The #MSVMMAJ_OUTPUT_FILE specifies the * @returns index of the string in the arguments if found, 0 * otherwise */ -int msvmmaj_check_argv(int argc, char **argv, char *str) +int gensvm_check_argv(int argc, char **argv, char *str) { int i; int arg_str = 0; @@ -69,7 +69,7 @@ int msvmmaj_check_argv(int argc, char **argv, char *str) * @returns index of the command line argument that corresponds to * the string, 0 if none matches. */ -int msvmmaj_check_argv_eq(int argc, char **argv, char *str) +int gensvm_check_argv_eq(int argc, char **argv, char *str) { int i; int arg_str = 0; @@ -88,19 +88,19 @@ int msvmmaj_check_argv_eq(int argc, char **argv, char *str) * * @details * This function is used to print a given string to the output stream - * specified by #MSVMMAJ_OUTPUT_FILE. The stream is flushed after the string - * is written to the stream. If #MSVMMAJ_OUTPUT_FILE is NULL, nothing is + * specified by #GENSVM_OUTPUT_FILE. The stream is flushed after the string + * is written to the stream. If #GENSVM_OUTPUT_FILE is NULL, nothing is * written. Note that this function is only used by note(), it should never be * used directly. * * @param[in] s string to write to the stream * */ -static void msvmmaj_print_string(const char *s) +static void gensvm_print_string(const char *s) { - if (MSVMMAJ_OUTPUT_FILE != NULL) { - fputs(s, MSVMMAJ_OUTPUT_FILE); - fflush(MSVMMAJ_OUTPUT_FILE); + if (GENSVM_OUTPUT_FILE != NULL) { + fputs(s, GENSVM_OUTPUT_FILE); + fflush(GENSVM_OUTPUT_FILE); } } @@ -111,7 +111,7 @@ static void msvmmaj_print_string(const char *s) * This function is a replacement of fprintf(), such that the output stream * does not have to be specified at each function call. The functionality is * exactly the same however. Writing the formatted string to the output stream - * is handled by msvmmaj_print_string(). + * is handled by gensvm_print_string(). * * @param[in] fmt String format * @param[in] ... variable argument list for the string format @@ -124,5 +124,5 @@ void note(const char *fmt,...) va_start(ap,fmt); vsprintf(buf,fmt,ap); va_end(ap); - (*msvmmaj_print_string)(buf); + (*gensvm_print_string)(buf); } -- cgit v1.2.3