From 24fedb504c3e96ffc0900cf09d3337ca72e2aabe Mon Sep 17 00:00:00 2001 From: houyang Date: Fri, 18 Jan 2008 17:46:31 +0000 Subject: [PATCH] --- fastlib/u/garryb/svm/smo.h | 196 +++++++-------- fastlib/u/garryb/svm/svm.cc | 214 +++++----------- fastlib/u/garryb/svm/svm.h | 483 ++++++------------------------------ fastlib/u/houyang/svm/svm.h | 67 +++-- 4 files changed, 267 insertions(+), 693 deletions(-) diff --git a/fastlib/u/garryb/svm/smo.h b/fastlib/u/garryb/svm/smo.h index 117ef62272..f4c64118e4 100644 --- a/fastlib/u/garryb/svm/smo.h +++ b/fastlib/u/garryb/svm/smo.h @@ -4,14 +4,13 @@ #include "fastlib/fastlib.h" /* TODO: I don't actually want these to be public */ -/* but sometimes we should provide freedoms for our advanced users */ const double SMO_ZERO = 1.0e-8; const double SMO_EPS = 1.0e-4; const double SMO_TOLERANCE = 1.0e-4; template class SMO { - FORBID_ACCIDENTAL_COPIES(SMO); + FORBID_COPY(SMO); public: typedef TKernel Kernel; @@ -20,11 +19,9 @@ class SMO { Matrix kernel_cache_sign_; Kernel kernel_; const Dataset *dataset_; - index_t n_data_; // number of data samples - index_t n_classes_; // number of classes - Matrix matrix_; // alias for the data matrix - Vector alpha_; // the alphas, to be optimized - index_t n_sv_; // number of support vectors + index_t n_data_; + Matrix matrix_; + Vector alpha_; Vector error_; double thresh_; double c_; @@ -40,16 +37,26 @@ class SMO { * * You must initialize separately the kernel. */ - void Init(int n_classes_in, double c_in, int budget_in) { + void Init(const Dataset* dataset_in, double c_in, int budget_in) { c_ = c_in; - n_classes_ = n_classes_in; - budget_ = budget_in; - thresh_ = 0.0; - n_sv_ = 0; + dataset_ = dataset_in; + matrix_.Alias(dataset_->matrix()); + + n_data_ = matrix_.n_cols(); + budget_ = min(budget_in, n_data_); + + alpha_.Init(n_data_); + alpha_.SetZero(); + sum_alpha_ = 0; + + error_.Init(n_data_); + error_.SetZero(); + + thresh_ = 0; } - void Train(const Dataset* dataset_in); + void Train(); const Kernel& kernel() const { return kernel_; @@ -63,15 +70,11 @@ class SMO { return thresh_; } - //index_t num_sv() const { - // return n_sv_; - //} - - void GetSVM(ArrayList &dataset_index, ArrayList &coef, ArrayList &sv_indicator); + void GetSVM(Matrix *support_vectors, Vector *support_alpha) const; private: index_t TrainIteration_(bool examine_all); - + bool TryChange_(index_t j); bool TakeStep_(index_t i, index_t j, double error_j); @@ -89,7 +92,6 @@ class SMO { return alpha <= 0 || alpha >= c_; } - // labels: the last row of the data matrix, 0 or 1 int GetLabelSign_(index_t i) const { return matrix_.get(matrix_.n_rows()-1, i) != 0 ? 1 : -1; } @@ -102,7 +104,7 @@ class SMO { double val; if (!IsBound_(alpha_[i])) { val = error_[i]; - VERBOSE_MSG(0, "error values %f and %f", error_[i], Evaluate_(i) - GetLabelSign_(i)); + DEBUG_MSG(0, "error values %f and %f", error_[i], Evaluate_(i) - GetLabelSign_(i)); } else { val = CalculateError_(i); } @@ -129,6 +131,7 @@ class SMO { Vector v_j; GetVector_(j, &v_j); double k = kernel_.Eval(v_i, v_j); + kernel_cache_sign_.set(j, i, k * GetLabelSign_(i) * GetLabelSign_(j)); } } @@ -136,50 +139,72 @@ class SMO { } }; -// Budget SMO training for 2-classes template -void SMO::Train(const Dataset* dataset_in) { +void SMO::GetSVM(Matrix *support_vectors, Vector *support_alpha) const { + index_t n_support = 0; + index_t i_support = 0; + + for (index_t i = 0; i < n_data_; i++) { + if (unlikely(alpha_[i] != 0)) { + n_support++; + } + } + + support_vectors->Init(matrix_.n_rows() - 1, n_support); + support_alpha->Init(n_support); + + for (index_t i = 0; i < n_data_; i++) { + if (unlikely(alpha_[i] != 0)) { + Vector source; + Vector dest; + + GetVector_(i, &source); + support_vectors->MakeColumnVector(i_support, &dest); + dest.CopyValues(source); + + (*support_alpha)[i_support] = alpha_[i] * GetLabelSign_(i); + + i_support++; + } + } +} + +template +double SMO::Evaluate_(index_t i) const { + // TODO: This only handles linear + Vector kernel_values; + double summation = 0; + + kernel_cache_sign_.MakeColumnVector(i, &kernel_values); + summation = la::Dot(alpha_, kernel_values) * GetLabelSign_(i); + + return (summation - thresh_); +} + +template +void SMO::Train() { bool examine_all = true; index_t num_changed = 0; int n_iter = 0; - - // data-dependent initialization - dataset_ = dataset_in; - matrix_.Alias(dataset_->matrix()); - n_data_ = matrix_.n_cols(); - budget_ = min(budget_, n_data_); - alpha_.Init(n_data_); - - alpha_.SetZero(); - sum_alpha_ = 0; - - error_.Init(n_data_); - error_.SetZero(); - - // calculate kernel_cache_sign_: [k_ij* y_i* y_j] CalcKernels_(); - while (num_changed > 0 || examine_all) { // TODO: other stopping criteria - VERBOSE_GOT_HERE(0); - // SMO iterations + while ((num_changed > 0 || examine_all)) { + DEBUG_GOT_HERE(0); num_changed = TrainIteration_(examine_all); if (examine_all) { examine_all = false; - } else if (num_changed == 0) { + } else if (num_changed == 0) { examine_all = true; } - // if exceed the maximum number of iterations, finished - // 200...TODO if (++n_iter == 200) { fprintf(stderr, "Max iterations %f!!!!!!!!!!!!!!!!!!!!!!!!!!\n", sum_alpha_); break; } - // for every 100 iterations, do budget shrink if (n_iter % 100 == 0 && budget_ < n_data_) { MinHeap alphas; @@ -202,16 +227,14 @@ void SMO::Train(const Dataset* dataset_in) { if (!IsBound_(alpha_[i])) { error_[i] = CalculateError_(i); } else { - error_[i] = 0; + error_[i] = 0; } sum_alpha_ += alpha_[i]; } - } // if - - } // while + } + } } -// SMO training iterations template index_t SMO::TrainIteration_(bool examine_all) { index_t num_changed = 0; @@ -225,16 +248,15 @@ index_t SMO::TrainIteration_(bool examine_all) { return num_changed; } -// try to find working set (maximal violating pair) template bool SMO::TryChange_(index_t j) { - double error_j = Error_(j); // -y_j g_j^* -thresh + double error_j = Error_(j); double rj = error_j * GetLabelSign_(j); - VERBOSE_GOT_HERE(0); + DEBUG_GOT_HERE(0); - if (!( (rj < -SMO_TOLERANCE && alpha_[j] < c_) - ||(rj > SMO_TOLERANCE && alpha_[j] > 0) )) { + if (!((rj < -SMO_TOLERANCE && alpha_[j] < c_) + || (rj > SMO_TOLERANCE && alpha_[j] > 0))) { return false; // nothing to change } @@ -245,11 +267,10 @@ bool SMO::TryChange_(index_t j) { double diff_max = 0; - // find the max(abs(y_i g_i^*)) for (index_t k = 0; k < n_data_; k++) { - if (!IsBound_(alpha_[k])) { // if 0 < alpha_[k] < c_ + if (!IsBound_(alpha_[k])) { double error_k = error_[k]; - double diff_k = fabs(error_k - error_j); // abs(y_j g_j^* - y_k g_k^*) + double diff_k = fabs(error_k - error_j); if (unlikely(diff_k > diff_max)) { diff_max = diff_k; i = k; @@ -261,7 +282,7 @@ bool SMO::TryChange_(index_t j) { } } - VERBOSE_GOT_HERE(0); + DEBUG_GOT_HERE(0); // try searching through non-bound examples index_t start_i = rand() % n_data_; index_t i = start_i; @@ -274,7 +295,7 @@ bool SMO::TryChange_(index_t j) { i = (i + 1) % n_data_; } while (i != start_i); - VERBOSE_GOT_HERE(0); + DEBUG_GOT_HERE(0); // try searching through all examples start_i = rand() % n_data_; i = start_i; @@ -289,11 +310,10 @@ bool SMO::TryChange_(index_t j) { return false; } -// search direction, update gradient template bool SMO::TakeStep_(index_t i, index_t j, double error_j) { if (i == j) { - VERBOSE_GOT_HERE(0); + DEBUG_GOT_HERE(0); return false; } @@ -325,8 +345,8 @@ bool SMO::TakeStep_(index_t i, index_t j, double error_j) { if (l >= u - SMO_TOLERANCE) { // TODO: might put in some tolerance - VERBOSE_MSG(0, "l=%f, u=%f, r=%f, c_=%f, s=%f", l, u, r, c_, s); - VERBOSE_GOT_HERE(0); + DEBUG_MSG(0, "l=%f, u=%f, r=%f, c_=%f, s=%f", l, u, r, c_, s); + DEBUG_GOT_HERE(0); return false; } @@ -334,18 +354,17 @@ bool SMO::TakeStep_(index_t i, index_t j, double error_j) { double kii = EvalKernel_(i, i); double kij = EvalKernel_(i, j); double kjj = EvalKernel_(j, j); - // second derivative of the objective function + // second derivative of objective function double eta = +2*kij - kii - kjj; - VERBOSE_MSG(0, "kij=%f, kii=%f, kjj=%f", kij, kii, kjj); + DEBUG_MSG(0, "kij=%f, kii=%f, kjj=%f", kij, kii, kjj); - // update alpha_j if (likely(eta < 0)) { - VERBOSE_MSG(0, "Common case"); + DEBUG_MSG(0, "Common case"); alpha_j = alpha_[j] - yj * (error_i - error_j) / eta; alpha_j = FixAlpha_(math::ClampRange(alpha_j, l, u)); } else { - VERBOSE_MSG(0, "Uncommon case"); + DEBUG_MSG(0, "Uncommon case"); //abort(); double c1 = eta/2; double c2 = yj * (error_i - error_j) - eta * alpha_j; @@ -368,11 +387,10 @@ bool SMO::TakeStep_(index_t i, index_t j, double error_j) { // check if there is progress if (fabs(delta_alpha_j) < SMO_EPS*(alpha_j + alpha_[j] + SMO_EPS)) { - VERBOSE_GOT_HERE(0); + DEBUG_GOT_HERE(0); return false; } - // update alpha_i alpha_i = alpha_i - (s)*(delta_alpha_j); if (alpha_i < SMO_ZERO) { alpha_j += s * alpha_i; @@ -403,7 +421,6 @@ bool SMO::TakeStep_(index_t i, index_t j, double error_j) { kernel_cache_sign_.MakeColumnVector(i, &kernel_i); kernel_cache_sign_.MakeColumnVector(j, &kernel_j); - // update gradient for (index_t k = 0; k < n_data_; k++) { if (likely(k != i) && likely(k != j) && !IsBound_(alpha_[k])) { error_[k] += (delta_alpha_i*kernel_i[k] + delta_alpha_j*kernel_j[k]) * GetLabelSign_(k) - delta_thresh; @@ -421,41 +438,8 @@ bool SMO::TakeStep_(index_t i, index_t j, double error_j) { error_[i] = 0; error_[j] = 0; - VERBOSE_GOT_HERE(0); + DEBUG_GOT_HERE(0); return true; } - -template -double SMO::Evaluate_(index_t i) const { - // TODO: This only handles linear - Vector kernel_values; - double summation = 0; - - kernel_cache_sign_.MakeColumnVector(i, &kernel_values); - summation = la::Dot(alpha_, kernel_values) * GetLabelSign_(i); - - return (summation - thresh_); -} - -// Get SVM results:coefficients, number and indecies of SVs -template -void SMO::GetSVM(ArrayList &dataset_index, ArrayList &coef, ArrayList &sv_indicator) { - coef.Init(n_data_); - for (index_t i = 0; i < n_data_; i++) { - if (alpha_[i] != 0) { // support vectors - coef[i] = alpha_[i] * GetLabelSign_(i); - sv_indicator[dataset_index[i]] = true; - n_sv_++; - //sv_index[i_sv]= dataset_bi_index[i]; - //GetVector_(i, &source); - //support_vectors->MakeColumnVector(i_support, &dest); - //dest.CopyValues(source); - } - else { - coef[i] = 0; - } - } -} - #endif diff --git a/fastlib/u/garryb/svm/svm.cc b/fastlib/u/garryb/svm/svm.cc index 0975e7a323..c109b254b4 100644 --- a/fastlib/u/garryb/svm/svm.cc +++ b/fastlib/u/garryb/svm/svm.cc @@ -26,6 +26,7 @@ void DoSvmNormalize(Dataset* dataset) { Matrix cov; + la::MulTransBInit(m, m, &cov); Vector d; @@ -33,7 +34,7 @@ void DoSvmNormalize(Dataset* dataset) { Matrix ui; // the inverse of eigenvectors //cov.PrintDebug("cov"); - PASSED(la::EigenvectorsInit(cov, &d, &u)); + la::EigenvectorsInit(cov, &d, &u); la::TransposeInit(u, &ui); for (index_t i = 0; i < d.length(); i++) { @@ -57,170 +58,73 @@ void DoSvmNormalize(Dataset* dataset) { } //dataset->matrix().PrintDebug("m"); - - if (fx_param_bool(NULL, "save", 0)) { - fx_default_param(NULL, "kfold/save", "1"); - dataset->WriteCsv("m_normalized.csv"); - } -} - -void CreateArtificialDataset(Dataset* dataset){ - Matrix m; - index_t n = fx_param_int(NULL, "n", 30); - double offset = fx_param_double(NULL, "offset", 0.0); - double range = fx_param_double(NULL, "range", 1.0); - double slope = fx_param_double(NULL, "slope", 1.0); - double margin = fx_param_double(NULL, "margin", 1.0); - double var = fx_param_double(NULL, "var", 1.0); - double intercept = fx_param_double(NULL, "intercept", 0.0); - - // 2 dimensional dataset, size n, 3 classes - m.Init(3, n); - for (index_t i = 0; i < n; i += 3) { - double x; - double y; - - x = (rand() * range / RAND_MAX) + offset; - y = margin / 2 + (rand() * var / RAND_MAX); - m.set(0, i, x); - m.set(1, i, x*slope + y + intercept); - m.set(2, i, 0); // labels - - x = (rand() * range / RAND_MAX) + offset; - y = margin / 2 + (rand() * var / RAND_MAX); - m.set(0, i+1, 10*x); - m.set(1, i+1, x*slope + y + intercept); - m.set(2, i+1, 1); // labels - - x = (rand() * range / RAND_MAX) + offset; - y = margin / 2 + (rand() * var / RAND_MAX); - m.set(0, i+2, 20*x); - m.set(1, i+2, x*slope + y + intercept); - m.set(2, i+2, 2); // labels - } - - data::Save("m.csv", m); - dataset->OwnMatrix(&m); -} - -int LoadData(Dataset* dataset, String datafilename){ - if (fx_param_exists(NULL, datafilename)) { - // when a data file is specified, use it. - if ( !PASSED(dataset->InitFromFile(fx_param_str_req(NULL, datafilename))) ) { - fprintf(stderr, "Couldn't open the data file.\n"); - return 0; - } - } - else { - fprintf(stderr, "No data file exist. Generating artificial dataset.\n"); - // otherwise, create an artificial dataset and save it to "m.csv" - CreateArtificialDataset(dataset); - } - - if (fx_param_bool(NULL, "normalize", 1)) { - fprintf(stderr, "Normalizing\n"); - DoSvmNormalize(dataset); - } else { - fprintf(stderr, "Skipping normalize\n"); - } - return 1; } int main(int argc, char *argv[]) { fx_init(argc, argv); - srand(time(NULL)); - - String mode = fx_param_str_req(NULL, "mode"); - String kernel = fx_param_str_req(NULL, "kernel"); - // TODO: more kernels to be supported - - // Cross Validation Mode, need cross validation data - if(mode == "cv") { - fprintf(stderr, "SVM Cross Validation... \n"); - - // Load cross validation data - Dataset cvset; - if (LoadData(&cvset, "cv_data") == 0) - return 1; - - if (kernel == "linear") { - SimpleCrossValidator< SVM > cross_validator; - // Initialize n_folds_, confusion_matrix_; k_cv: number of cross-validation folds, need k_cv>1 - cross_validator.Init(&cvset,cvset.n_labels(),fx_param_int_req(NULL,"k_cv"), fx_root, "svm"); - // k_cv folds cross validation; (true): do training set permutation - cross_validator.Run(true); - cross_validator.confusion_matrix().PrintDebug("confusion matrix"); - } - else if (kernel == "gaussian") { - SimpleCrossValidator< SVM > cross_validator; - // Initialize n_folds_, confusion_matrix_; k_cv: number of cross-validation folds - cross_validator.Init(&cvset,cvset.n_labels(),fx_param_int_req(NULL,"k_cv"), fx_root, "svm"); - // k_cv folds cross validation; (true): do training set permutation - cross_validator.Run(true); - cross_validator.confusion_matrix().PrintDebug("confusion matrix"); - } - } - // Training Mode, need training data | Training + Testing(online) Mode, need training data + testing data - else if (mode=="train" || mode=="train_test"){ - fprintf(stderr, "SVM Training... \n"); - - // Load training data - Dataset trainset; - if (LoadData(&trainset, "train_data") == 0) // TODO:param_req + Dataset dataset; + + if (fx_param_exists(NULL, "data")) { + // if a data file is specified, use it. + if (!PASSED(dataset.InitFromFile(fx_param_str_req(NULL, "data")))) { + fprintf(stderr, "Couldn't open the data file.\n"); return 1; + } + } else { + // create an artificial dataset and save it to "m.csv" - // Begin SVM Training | Training and Testing - datanode *svm_module = fx_submodule(fx_root, NULL, "svm"); - - if (kernel == "linear") { - SVM svm; - svm.InitTrain(trainset, trainset.n_labels(), svm_module); - if (mode=="train_test"){ // training and testing, thus no need to load model from file - fprintf(stderr, "SVM Classifying... \n"); - // Load testing data - Dataset testset; - if (LoadData(&testset, "test_data") == 0) // TODO:param_req - return 1; - svm.BatchClassify(&testset, "test_labels"); - } - } - else if (kernel == "gaussian") { - SVM svm; - svm.InitTrain(trainset, trainset.n_labels(), svm_module); - if (mode=="train_test"){ // training and testing, thus no need to load model from file - fprintf(stderr, "SVM Classifying... \n"); - // Load testing data - Dataset testset; - if (LoadData(&testset, "test_data") == 0) // TODO:param_req - return 1; - svm.BatchClassify(&testset, "test_labels"); // TODO:param_req - } + Matrix m; + index_t n = fx_param_int(NULL, "n", 30); + double offset = fx_param_double(NULL, "offset", 0.0); + double range = fx_param_double(NULL, "range", 1.0); + double slope = fx_param_double(NULL, "slope", 1.0); + double margin = fx_param_double(NULL, "margin", 1.0); + double var = fx_param_double(NULL, "var", 1.0); + double intercept = fx_param_double(NULL, "intercept", 0.0); + + // 3 dimensional dataset, size n + m.Init(3, n); + + for (index_t i = 0; i < n; i += 2) { + double x; + double y; + + x = (rand() * range / RAND_MAX) + offset; + y = margin / 2 + (rand() * var / RAND_MAX); + m.set(0, i, x); + m.set(1, i, x*slope + y + intercept); + m.set(2, i, 0); + + x = (rand() * range / RAND_MAX) + offset; + y = margin / 2 + (rand() * var / RAND_MAX); + m.set(0, i+1, x); + m.set(1, i+1, x*slope - y + intercept); + m.set(2, i+1, 1); } + data::Save("m.csv", m); + dataset.OwnMatrix(&m); } - // Testing(offline) Mode, need loading model file and testing data - else if (mode=="test") { - fprintf(stderr, "SVM Classifying... \n"); - - // Load testing data - Dataset testset; - if (LoadData(&testset, "test_data") == 0) // TODO:param_req - return 1; - - // Begin Classification - datanode *svm_module = fx_submodule(fx_root, NULL, "svm"); - - if (kernel == "linear") { - SVM svm; - svm.Init(testset, testset.n_labels(), svm_module); // TODO:n_labels() -> num_classes_ - svm.LoadModelBatchClassify(&testset, "svm_model", "test_labels"); // TODO:param_req - } - else if (kernel == "gaussian") { - SVM svm; - svm.Init(testset, testset.n_labels(), svm_module); // TODO:n_labels() -> num_classes_ - svm.LoadModelBatchClassify(&testset, "svm_model", "test_labels"); // TODO:param_req - } + + if (fx_param_bool(NULL, "normalize", 1)) { + fprintf(stderr, "Normalizing\n"); + DoSvmNormalize(&dataset); + } else { + fprintf(stderr, "Skipping normalize\n"); } + + if (fx_param_bool(NULL, "save", 0)) { + fx_default_param(NULL, "kfold/save", "1"); + dataset.WriteCsv("normalized.csv"); + } + + SimpleCrossValidator< SVM > cross_validator; + // k_cv: number of cross-validation folds + cross_validator.Init(&dataset, 2,fx_param_int_req(NULL,"k_cv"), fx_root, "svm"); + cross_validator.Run(true); + + cross_validator.confusion_matrix().PrintDebug("confusion matrix"); + fx_done(); } diff --git a/fastlib/u/garryb/svm/svm.h b/fastlib/u/garryb/svm/svm.h index b80d7bad42..8b63bc9573 100644 --- a/fastlib/u/garryb/svm/svm.h +++ b/fastlib/u/garryb/svm/svm.h @@ -7,450 +7,107 @@ #include -#define ID_LINEAR 0 -#define ID_GAUSSIAN 1 - struct SVMLinearKernel { - ArrayList kpara_; // kernel parameters - void Init(datanode *node) { //TODO:node - kpara_.Init(); - } - void GetName(String* kname) { - kname->Copy("linear"); - } - int GetTypeId() { - return ID_LINEAR; - } + void Init(datanode *node) {} + + void Copy(const SVMLinearKernel& other) {} + double Eval(const Vector& a, const Vector& b) const { return la::Dot(a, b); } - void SaveParam(FILE* fp) { - } }; class SVMRBFKernel { + private: + double sigma_; + double gamma_; + public: - ArrayList kpara_; // kernel parameters - void Init(datanode *node) { //TODO: node.. - kpara_.Init(2); - kpara_[0] = fx_param_double_req(NULL, "sigma"); //sigma - kpara_[1] = -1.0 / (2 * math::Sqr(kpara_[0])); //gamma + + void Init(datanode *node) { + sigma_ = fx_param_double_req(NULL, "sigma"); + gamma_ = -1.0 / (2 * math::Sqr(sigma_)); } - void GetName(String* kname) { - kname->Copy("gaussian"); - } - int GetTypeId() { - return ID_GAUSSIAN; + + void Copy(const SVMRBFKernel& other) { + sigma_ = other.sigma_; + gamma_ = other.gamma_; } + double Eval(const Vector& a, const Vector& b) const { double distance_squared = la::DistanceSqEuclidean(a, b); - return exp(kpara_[1] * distance_squared); + return exp(gamma_ * distance_squared); } - void SaveParam(FILE* fp) { - fprintf(fp, "sigma %g\n", kpara_[0]); - fprintf(fp, "gamma %g\n", kpara_[1]); + + double sigma() const { + return sigma_; } }; template class SVM { - private: - // models for the binary classifiers - struct SVM_MODELS { - double thresh_; - ArrayList bi_coef_; // all coefficients of the binary dataset, not necessarily thoes of SVs - }; - ArrayList models_; - - ArrayList train_labels_list_; // set of labels, need to be integers - - // total set of support vectors and their coefficients - Matrix sv_; - Matrix sv_coef_; - - index_t total_num_sv_; - ArrayList sv_indicator_; - ArrayList sv_index_; - ArrayList sv_list_startpos_; - ArrayList sv_list_ct_; - - struct SVM_PARAMETERS { - TKernel kernel_; - String kernelname_; - int kerneltypeid_; - double c_; - int b_; - }; - SVM_PARAMETERS param_; // same for every binary model - - int num_classes_; - int num_models_; - int num_features_; - public: typedef TKernel Kernel; - - void Init(const Dataset& dataset, int n_classes, datanode *module); + + private: + TKernel kernel_; + double c_; + int b_; + Matrix support_vectors_; + double thresh_; + Vector alpha_; + + public: void InitTrain(const Dataset& dataset, int n_classes, datanode *module); - void SaveModel(String modelfilename); - void LoadModel(Dataset* testset, String modelfilename); + int Classify(const Vector& vector); - void BatchClassify(Dataset* testset, String testlabelfilename); - void LoadModelBatchClassify(Dataset* testset, String modelfilename, String testlabelfilename); }; template -void SVM::Init(const Dataset& dataset, int n_classes, datanode *module){ - models_.Init(); - sv_indicator_.Init(dataset.n_points()); - sv_index_.Init(); +void SVM::InitTrain( + const Dataset& dataset, int n_classes, datanode *module) { + DEBUG_ASSERT_MSG(n_classes == 2, "SVM is only a binary classifier"); + + fx_set_param(module, "kernel_type", typeid(TKernel).name()); + kernel_.Init(fx_submodule(module, "kernel", "kernel")); + c_ = fx_param_double(module, "c", 1.0); + b_ = fx_param_int(module, "b", dataset.n_points()); + + SMO smo; + smo.Init(&dataset, c_, b_); + smo.kernel().Copy(kernel_); + smo.Train(); + + thresh_ = smo.threshold(); + smo.GetSVM(&support_vectors_, &alpha_); + DEBUG_ASSERT(alpha_.length() != 0); + DEBUG_ASSERT(alpha_.length() == support_vectors_.n_cols()); + + DEBUG_ONLY(fprintf(stderr, "----------------------\n")); + DEBUG_ONLY(support_vectors_.PrintDebug("support vectors")); + DEBUG_ONLY(alpha_.PrintDebug("support vector weights")); + DEBUG_ONLY(fprintf(stderr, "-- THRESHOLD: %f\n", thresh_)); - param_.kernel_.Init(fx_submodule(module, "kernel", "kernel")); - param_.kernel_.GetName(¶m_.kernelname_); - param_.kerneltypeid_ = param_.kernel_.GetTypeId(); - - param_.c_ = fx_param_double_req(NULL, "c"); - // budget parameter, contorls # of support vectors; default: # of data samples (use all) - param_.b_ = fx_param_int(module, "b", dataset.n_points()); // TODO: param_req - - num_classes_ = n_classes; - num_models_ = num_classes_ * (num_classes_-1) / 2; - num_features_ = 0; - sv_list_startpos_.Init(n_classes); - sv_list_ct_.Init(n_classes); - - for (index_t i=0; i -void SVM::InitTrain(const Dataset& dataset, int n_classes, datanode *module) { - Init(dataset, n_classes, module); - num_features_ = dataset.n_features()-1; // last column in dataset is for labels - - // Group labels, split the training dataset for training bi-class SVM classifiers - ArrayList train_labels_index; // e.g. [c1[0,5,6,7,10,13,17],c2[1,2,4,8,9],c3[...]] - ArrayList train_labels_ct; // counter, e.g. [7,5,8] - ArrayList train_labels_startpos; // start position, e.g. [0,7,12] - dataset.GetLabels(train_labels_list_, train_labels_index, train_labels_ct, train_labels_startpos); - // Train n_classes*(n_classes-1)/2 bi-class(labels: 0, 1) models using SMO - index_t ct = 0; - index_t i; index_t j; - for (i = 0; i < n_classes; i++) { - for (j = i+1; j < n_classes; j++) { - models_.AddBack(); - - // Initialize parameters c_, budget_, alpha_, error_, thresh_ - SMO smo; - smo.Init(n_classes, param_.c_, param_.b_); - smo.kernel().Init(fx_submodule(module, "kernel", "kernel")); - // Construct dataset consists of two classes i and j (reassign labels 0 and 1) - Dataset dataset_bi; - dataset_bi.InitBlank(); - dataset_bi.info().Init(); - dataset_bi.matrix().Init(num_features_+1, train_labels_ct[i]+train_labels_ct[j]); - ArrayList dataset_bi_index; - dataset_bi_index.Init(train_labels_ct[i]+train_labels_ct[j]); - - for (index_t m = 0; m < train_labels_ct[i]; m++) { - Vector source, dest; - dataset_bi.matrix().MakeColumnVector(m, &dest); - dataset.matrix().MakeColumnVector(train_labels_index[train_labels_startpos[i]+m], &source); - dest.CopyValues(source); - dataset_bi.matrix().set(num_features_, m, 0); // last row for labels 0 - dataset_bi_index[m] = train_labels_index[train_labels_startpos[i]+m]; - } - for (index_t n = 0; n < train_labels_ct[j]; n++) { - Vector source, dest; - dataset_bi.matrix().MakeColumnVector(n+train_labels_ct[i], &dest); - dataset.matrix().MakeColumnVector(train_labels_index[train_labels_startpos[j]+n], &source); - dest.CopyValues(source); - dataset_bi.matrix().set(num_features_, n+train_labels_ct[i], 1); // last row for labels 1 - dataset_bi_index[n+train_labels_ct[i]] = train_labels_index[train_labels_startpos[j]+n]; - } - // 2-classes SVM-SMO training - smo.Train(&dataset_bi); - // Get the trained bi-class model - models_[ct].thresh_ = smo.threshold(); - smo.GetSVM(dataset_bi_index, models_[ct].bi_coef_, sv_indicator_); - //smo.GetSVM(dataset_bi_index, models_[ct].bi_coef_, sv_indicator_, &(models_[ct].num_sv_bi_)); - - //DEBUG_ASSERT(models_[ct].alpha_.length() != 0); - //DEBUG_ASSERT(models_[ct].alpha_.length() == models_[ct].support_vectors_.n_cols()); - //DEBUG_ONLY(fprintf(stderr, "----------------------\n")); - //DEBUG_ONLY(models_[ct].support_vectors_.PrintDebug("support vectors")); - //DEBUG_ONLY(models_[ct].alpha_.PrintDebug("support vector weights")); - //DEBUG_ONLY(fprintf(stderr, "-- THRESHOLD: %f\n", models_[ct].thresh_)); - //fx_format_result(module, "n_support", "%"LI"d", models_[ct].alpha_.length()); - //fx_format_result(module, "n_support", "%"LI"d", models_[ct].sv_bi_coef_.size()); - ct++; - } - } - // Get total set of SVs; - index_t k; - sv_list_startpos_[0] = 0; - total_num_sv_ = 0; - for (i = 0; i < n_classes; i++) { - ct = 0; - for (j = 0; j < train_labels_ct[i]; j++) { - if (sv_indicator_[train_labels_startpos[i] + j]) { - *sv_index_.AddBack() = train_labels_index[train_labels_startpos[i] + j]; - total_num_sv_++; - ct++; - } - } - sv_list_ct_[i] = ct; - if (i >= 1) - sv_list_startpos_[i] = sv_list_startpos_[i-1] + sv_list_ct_[i-1]; - } - sv_.Init(num_features_, total_num_sv_); - for (i = 0; i < total_num_sv_; i++) { - Vector source, dest; - sv_.MakeColumnVector(i, &dest); - dataset.matrix().MakeColumnSubvector(sv_index_[i], 0, num_features_, &source); // last row of dataset is for labels - dest.CopyValues(source); - } - // Get coefficients for the total set of SVs, i.e. models_[x].bi_coef_ -> sv_coef_ - index_t p; - index_t ct_model = 0; - sv_coef_.Init(n_classes-1, total_num_sv_); - sv_coef_.SetZero(); - for (i = 0; i < n_classes; i++) { - for (j = i+1; j < n_classes; j++) { - p = sv_list_startpos_[i]; - for (k = 0; k < train_labels_ct[i]; k++) { - if (sv_indicator_[train_labels_startpos[i]+k]) { - sv_coef_.set(j-1, p++, models_[ct_model].bi_coef_[k]); - } - } - p = sv_list_startpos_[j]; - for (k = 0; k < train_labels_ct[j]; k++) { - if (sv_indicator_[train_labels_startpos[j]+k]) { - sv_coef_.set(i, p++, models_[ct_model].bi_coef_[k]); - } - } - ct_model++; - } - } - // Save models - SaveModel("svm_model"); // TODO: param_req, and for CV mode -} - -// Save SVM model to file, TODO: use XML -template -void SVM::SaveModel(String modelfilename) { - FILE *fp = fopen(modelfilename,"w"); - if (fp==NULL) - fprintf(stderr, "Cannot save trained model to file!"); - index_t i; index_t j; - - fprintf(fp, "svm_type svm_c\n"); // TODO: svm-mu, svm-regression... - fprintf(fp, "num_classes %d\n", num_classes_); // TODO: only for svm_c - fprintf(fp, "kernel_name %s\n", param_.kernelname_.c_str()); - fprintf(fp, "kernel_typeid %d\n", param_.kerneltypeid_); - // save kernel parameters - param_.kernel_.SaveParam(fp); - fprintf(fp, "total_num_sv %d\n", total_num_sv_); - fprintf(fp, "labels "); - for (i = 0; i < num_classes_; i++) - fprintf(fp, "%d ", train_labels_list_[i]); - fprintf(fp, "\n"); - // save models - fprintf(fp, "thresholds "); - for (i = 0; i < num_models_; i++) - fprintf(fp, "%f ", models_[i].thresh_); - fprintf(fp, "\n"); - fprintf(fp, "sv_list_startpos "); - for (i =0; i < num_classes_; i++) - fprintf(fp, "%d ", sv_list_startpos_[i]); - fprintf(fp, "\n"); - fprintf(fp, "sv_list_ct "); - for (i =0; i < num_classes_; i++) - fprintf(fp, "%d ", sv_list_ct_[i]); - fprintf(fp, "\n"); - // save coefficients and support vectors - fprintf(fp, "SV_coefs\n"); - for (i = 0; i < total_num_sv_; i++) { - for (j = 0; j < num_classes_-1; j++) { - fprintf(fp, "%f ", sv_coef_.get(j,i)); - } - fprintf(fp, "\n"); - } - fprintf(fp, "SVs\n"); - for (i = 0; i < total_num_sv_; i++) { - for (j = 0; j < num_features_; j++) { // n_rows-1 - fprintf(fp, "%f ", sv_.get(j,i)); - } - fprintf(fp, "\n"); - } - fclose(fp); -} - -// Load SVM model file, TODO: use XML -template -void SVM::LoadModel(Dataset* testset, String modelfilename) { - //Init - train_labels_list_.Init(num_classes_); - num_features_ = testset->n_features(); - - //load model file - FILE *fp = fopen(modelfilename, "r"); - char cmd[80]; - int i, j; int temp_d; double temp_f; - for (i = 0; i < num_models_; i++) { - models_.AddBack(); - } - while (1) { - fscanf(fp,"%80s",cmd); - if(strcmp(cmd,"svm_type")==0) { - fscanf(fp,"%80s",cmd); - if(strcmp(cmd,"svm_c")==0) { - fprintf(stderr, "SVM_C\n"); - } - } - else if (strcmp(cmd, "num_classes")==0) { - fscanf(fp,"%d",&num_classes_); - } - else if (strcmp(cmd, "kernel_name")==0) { - fscanf(fp,"%80s",param_.kernelname_.c_str()); - } - else if (strcmp(cmd, "kernel_typeid")==0) { - fscanf(fp,"%d",¶m_.kerneltypeid_); - } - else if (strcmp(cmd, "sigma")==0) { - fscanf(fp,"%lf",¶m_.kernel_.kpara_[0]); // for gaussian kernels - } - else if (strcmp(cmd, "gamma")==0) { - fscanf(fp,"%lf",¶m_.kernel_.kpara_[1]); // for gaussian kernels - } - else if (strcmp(cmd, "total_num_sv")==0) { - fscanf(fp,"%d",&total_num_sv_); - } - else if (strcmp(cmd, "labels")==0) { - for (i=0; i int SVM::Classify(const Vector& datum) { - index_t i; index_t j; index_t k; - ArrayList keval; - keval.Init(total_num_sv_); - for (i = 0; i < total_num_sv_; i++) { - Vector support_vector_i; - sv_.MakeColumnVector(i, &support_vector_i); - keval[i] = param_.kernel_.Eval(datum, support_vector_i); + double summation = 0; + + for (index_t i = 0; i < alpha_.length(); i++) { + Vector support_vector; + support_vectors_.MakeColumnVector(i, &support_vector); + double term = alpha_[i] * kernel_.Eval(datum, support_vector); + + DEBUG_MSG(0, "alpha %f, term %f", alpha_[i], term); + + summation += term; } - ArrayList values; - values.Init(num_models_); - index_t ct = 0; - for (i = 0; i < num_classes_; i++) { - for (j = i+1; j < num_classes_; j++) { - double sum = 0; - for(k = 0; k < sv_list_ct_[i]; k++) { - sum += sv_coef_.get(j-1, sv_list_startpos_[i]+k) * keval[sv_list_startpos_[i]+k]; - } - for(k = 0; k < sv_list_ct_[j]; k++) - sum += sv_coef_.get(i, sv_list_startpos_[j]+k) * keval[sv_list_startpos_[j]+k]; - sum -= models_[ct].thresh_; - values[ct] = sum; - fprintf(stderr, "%f\n", values[ct]); - ct++; - } - } - - ArrayList vote; - vote.Init(num_classes_); - ct = 0; - for (i = 0; i < num_classes_; i++) { - for (j = i+1; j < num_classes_; j++) { - if(values[ct] > 0) - ++vote[i]; - else - ++vote[j]; - ct++; - } - } - index_t vote_max_idx = 0; - for (i = 1; i < num_classes_; i++) { - if (vote[i] > vote[vote_max_idx]) - vote_max_idx = i; - } - return train_labels_list_[vote_max_idx]; -} - -// Batch classification for multiple testing vectors, no need to load model file -// Note: for test set, if no true test labels provided, just put some dummy labels (e.g. all -1) in the last row of testset -template -void SVM::BatchClassify(Dataset* testset, String testlablefilename) { - FILE *fp = fopen(testlablefilename,"w"); - if (fp==NULL) - fprintf(stderr, "Cannot save test labels to file!"); - num_features_ = testset->n_features()-1; - for (index_t i=0; in_points(); i++) { - Vector testvec; - testset->matrix().MakeColumnSubvector(i, 0, num_features_, &testvec); - int testlabel = Classify(testvec); - fprintf(fp, "%d\n", testlabel); - } - fclose(fp); -} - -// Load models from a file, and perform batch classification for multiple testing vectors -template -void SVM::LoadModelBatchClassify(Dataset* testset, String modelfilename, String testlabelfilename) { - LoadModel(testset, modelfilename); - BatchClassify(testset, testlabelfilename); + + DEBUG_MSG(0, "summation=%f, thresh_=%f", summation, thresh_); + + return (summation - thresh_ > 0.0) ? 1 : 0; } #endif diff --git a/fastlib/u/houyang/svm/svm.h b/fastlib/u/houyang/svm/svm.h index a1f69d5870..b80d7bad42 100644 --- a/fastlib/u/houyang/svm/svm.h +++ b/fastlib/u/houyang/svm/svm.h @@ -108,7 +108,7 @@ void SVM::Init(const Dataset& dataset, int n_classes, datanode *module) param_.kernel_.Init(fx_submodule(module, "kernel", "kernel")); param_.kernel_.GetName(¶m_.kernelname_); param_.kerneltypeid_ = param_.kernel_.GetTypeId(); - // c; default:1 + param_.c_ = fx_param_double_req(NULL, "c"); // budget parameter, contorls # of support vectors; default: # of data samples (use all) param_.b_ = fx_param_int(module, "b", dataset.n_points()); // TODO: param_req @@ -255,25 +255,33 @@ void SVM::SaveModel(String modelfilename) { param_.kernel_.SaveParam(fp); fprintf(fp, "total_num_sv %d\n", total_num_sv_); fprintf(fp, "labels "); - for (i=0; i::LoadModel(Dataset* testset, String modelfilename) { //Init train_labels_list_.Init(num_classes_); num_features_ = testset->n_features(); + //load model file FILE *fp = fopen(modelfilename, "r"); char cmd[80]; int i, j; int temp_d; double temp_f; - for (i=0; i::LoadModel(Dataset* testset, String modelfilename) { if(strcmp(cmd,"svm_type")==0) { fscanf(fp,"%80s",cmd); if(strcmp(cmd,"svm_c")==0) { - fprintf(stderr, "SVM Classification"); + fprintf(stderr, "SVM_C\n"); } } else if (strcmp(cmd, "num_classes")==0) { @@ -327,29 +336,48 @@ void SVM::LoadModel(Dataset* testset, String modelfilename) { } } else if (strcmp(cmd, "thresholds")==0) { - for (i=0; i::Classify(const Vector& datum) { for (i = 0; i < num_classes_; i++) { for (j = i+1; j < num_classes_; j++) { double sum = 0; - for(k = 0; k < sv_list_ct_[i]; k++) + for(k = 0; k < sv_list_ct_[i]; k++) { sum += sv_coef_.get(j-1, sv_list_startpos_[i]+k) * keval[sv_list_startpos_[i]+k]; + } for(k = 0; k < sv_list_ct_[j]; k++) sum += sv_coef_.get(i, sv_list_startpos_[j]+k) * keval[sv_list_startpos_[j]+k]; sum -= models_[ct].thresh_;