diff --git a/fastlib/u/garryb/svm/build.py b/fastlib/u/garryb/svm/build.py deleted file mode 100644 index 0071909434..0000000000 --- a/fastlib/u/garryb/svm/build.py +++ /dev/null @@ -1,7 +0,0 @@ - -binrule( - name = "svm_bin", - sources = ["svm.cc"], - headers = ["smo.h", "svm.h"], - linkables = ["fastlib:fastlib"]) - diff --git a/fastlib/u/garryb/svm/smo.h b/fastlib/u/garryb/svm/smo.h deleted file mode 100644 index f4c64118e4..0000000000 --- a/fastlib/u/garryb/svm/smo.h +++ /dev/null @@ -1,445 +0,0 @@ -#ifndef U_SVM_SMO_H -#define U_SVM_SMO_H - -#include "fastlib/fastlib.h" - -/* TODO: I don't actually want these to be public */ -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_COPY(SMO); - - public: - typedef TKernel Kernel; - - private: - Matrix kernel_cache_sign_; - Kernel kernel_; - const Dataset *dataset_; - index_t n_data_; - Matrix matrix_; - Vector alpha_; - Vector error_; - double thresh_; - double c_; - int budget_; - double sum_alpha_; - - public: - SMO() {} - ~SMO() {} - - /** - * Initializes an SMO problem. - * - * You must initialize separately the kernel. - */ - void Init(const Dataset* dataset_in, double c_in, int budget_in) { - c_ = c_in; - - 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 Kernel& kernel() const { - return kernel_; - } - - Kernel& kernel() { - return kernel_; - } - - double threshold() const { - return thresh_; - } - - 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); - - double FixAlpha_(double alpha) const { - if (alpha < SMO_ZERO) { - alpha = 0; - } else if (alpha > c_ - SMO_ZERO) { - alpha = c_; - } - return alpha; - } - - bool IsBound_(double alpha) const { - return alpha <= 0 || alpha >= c_; - } - - int GetLabelSign_(index_t i) const { - return matrix_.get(matrix_.n_rows()-1, i) != 0 ? 1 : -1; - } - - void GetVector_(index_t i, Vector *v) const { - matrix_.MakeColumnSubvector(i, 0, matrix_.n_rows()-1, v); - } - - double Error_(index_t i) const { - double val; - if (!IsBound_(alpha_[i])) { - val = error_[i]; - DEBUG_MSG(0, "error values %f and %f", error_[i], Evaluate_(i) - GetLabelSign_(i)); - } else { - val = CalculateError_(i); - } - return val; - } - - double CalculateError_(index_t i) const { - return Evaluate_(i) - GetLabelSign_(i); - } - - double Evaluate_(index_t i) const; - - double EvalKernel_(index_t i, index_t j) const { - return kernel_cache_sign_.get(i, j) * (GetLabelSign_(i) * GetLabelSign_(j)); - } - - void CalcKernels_() { - kernel_cache_sign_.Init(n_data_, n_data_); - fprintf(stderr, "Kernel Start\n"); - for (index_t i = 0; i < n_data_; i++) { - for (index_t j = 0; j < n_data_; j++) { - Vector v_i; - GetVector_(i, &v_i); - 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)); - } - } - fprintf(stderr, "Kernel Stop\n"); - } -}; - -template -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; - - CalcKernels_(); - - 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) { - examine_all = true; - } - - if (++n_iter == 200) { - fprintf(stderr, "Max iterations %f!!!!!!!!!!!!!!!!!!!!!!!!!!\n", - sum_alpha_); - break; - } - - if (n_iter % 100 == 0 && budget_ < n_data_) { - MinHeap alphas; - - alphas.Init(); - for (index_t i = 0; i < n_data_; i++) { - alphas.Put(-alpha_[i], i); - } - - for (index_t i = 0; i < budget_; i++) { - alphas.Pop(); - } - - while (alphas.size() > 0) { - alpha_[alphas.Pop()] = 0; - } - - sum_alpha_ = 0; - - for (index_t i = 0; i < n_data_; i++) { - if (!IsBound_(alpha_[i])) { - error_[i] = CalculateError_(i); - } else { - error_[i] = 0; - } - sum_alpha_ += alpha_[i]; - } - } - } -} - -template -index_t SMO::TrainIteration_(bool examine_all) { - index_t num_changed = 0; - - for (index_t i = 0; i < n_data_; i++) { - if ((examine_all || !IsBound_(alpha_[i])) && TryChange_(i)) { - num_changed++; - } - } - - return num_changed; -} - -template -bool SMO::TryChange_(index_t j) { - double error_j = Error_(j); - double rj = error_j * GetLabelSign_(j); - - DEBUG_GOT_HERE(0); - - if (!((rj < -SMO_TOLERANCE && alpha_[j] < c_) - || (rj > SMO_TOLERANCE && alpha_[j] > 0))) { - return false; // nothing to change - } - - // first try the one we suspect to have the largest yield - - if (error_j != 0) { - index_t i = -1; - - double diff_max = 0; - - for (index_t k = 0; k < n_data_; k++) { - if (!IsBound_(alpha_[k])) { - double error_k = error_[k]; - double diff_k = fabs(error_k - error_j); - if (unlikely(diff_k > diff_max)) { - diff_max = diff_k; - i = k; - } - } - } - if (i != -1 && TakeStep_(i, j, error_j)) { - return true; - } - } - - DEBUG_GOT_HERE(0); - // try searching through non-bound examples - index_t start_i = rand() % n_data_; - index_t i = start_i; - - do { - //if (alpha_[i] != 0 && TakeStep_(i, j, error_j)) { - if (!IsBound_(alpha_[i]) && TakeStep_(i, j, error_j)) { - return true; - } - i = (i + 1) % n_data_; - } while (i != start_i); - - DEBUG_GOT_HERE(0); - // try searching through all examples - start_i = rand() % n_data_; - i = start_i; - - do { - if (IsBound_(alpha_[i]) && TakeStep_(i, j, error_j)) { - return true; - } - i = (i + 1) % n_data_; - } while (i != start_i); - - return false; -} - -template -bool SMO::TakeStep_(index_t i, index_t j, double error_j) { - if (i == j) { - DEBUG_GOT_HERE(0); - return false; - } - - int yi = GetLabelSign_(i); - int yj = GetLabelSign_(j); - double alpha_i = alpha_[i]; - double alpha_j = alpha_[j]; - double delta_thresh; - double l; - double u; - int s = (yi == yj) ? 1 : -1; - double error_i = Error_(i); - double r; - double budget_upper_bound; - - if (s < 0) { - DEBUG_ASSERT(s == -1); - r = alpha_j - alpha_i; // target values are not equal - //budget_upper_bound = (budget_*c_ - sum_alpha_ + 2*alpha_j) / 2; - double gamma = alpha_i - alpha_j; - budget_upper_bound = (gamma - budget_*c_ + sum_alpha_ - alpha_i - alpha_j) / (-2); - } else { - r = alpha_j + alpha_i - c_; // target values are equal - budget_upper_bound = DBL_MAX; - } - - l = math::ClampNonNegative(r); - u = min(c_ + math::ClampNonPositive(r), budget_upper_bound); - - if (l >= u - SMO_TOLERANCE) { - // TODO: might put in some tolerance - DEBUG_MSG(0, "l=%f, u=%f, r=%f, c_=%f, s=%f", l, u, r, c_, s); - DEBUG_GOT_HERE(0); - return false; - } - - // cached kernel values - double kii = EvalKernel_(i, i); - double kij = EvalKernel_(i, j); - double kjj = EvalKernel_(j, j); - // second derivative of objective function - double eta = +2*kij - kii - kjj; - - DEBUG_MSG(0, "kij=%f, kii=%f, kjj=%f", kij, kii, kjj); - - if (likely(eta < 0)) { - 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 { - DEBUG_MSG(0, "Uncommon case"); - //abort(); - double c1 = eta/2; - double c2 = yj * (error_i - error_j) - eta * alpha_j; - double objlower = c1*l*l + c2*l; - double objupper = c1*u*u + c2*u; - - if (objlower > objupper + SMO_EPS) { - alpha_j = l; - } else if (objlower < objupper - SMO_EPS) { - alpha_j = u; - } else { - alpha_j = alpha_[j]; - } - DEBUG_ASSERT(alpha_j == FixAlpha_(alpha_j)); - } - - alpha_j = FixAlpha_(alpha_j); - - double delta_alpha_j = alpha_j - alpha_[j]; - - // check if there is progress - if (fabs(delta_alpha_j) < SMO_EPS*(alpha_j + alpha_[j] + SMO_EPS)) { - DEBUG_GOT_HERE(0); - return false; - } - - alpha_i = alpha_i - (s)*(delta_alpha_j); - if (alpha_i < SMO_ZERO) { - alpha_j += s * alpha_i; - alpha_i = 0; - } else if (alpha_i > c_ - SMO_ZERO) { - double t = alpha_i - c_; - alpha_j += s * t; - alpha_i = c_; - } - - double delta_alpha_i = alpha_i - alpha_[i]; - - // calculate threshold - double delta_thresh_i = error_i + yi*delta_alpha_i*kii + yj*delta_alpha_j*kij; - double delta_thresh_j = error_j + yi*delta_alpha_i*kij + yj*delta_alpha_j*kjj; - - if (!IsBound_(alpha_i)) { - delta_thresh = delta_thresh_i; - } else if (!IsBound_(alpha_j)) { - delta_thresh = delta_thresh_j; - } else { - delta_thresh = (delta_thresh_i + delta_thresh_j) / 2.0; - } - - Vector kernel_i; - Vector kernel_j; - - kernel_cache_sign_.MakeColumnVector(i, &kernel_i); - kernel_cache_sign_.MakeColumnVector(j, &kernel_j); - - 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; - } - } - - thresh_ += delta_thresh; - - alpha_[i] = alpha_i; - alpha_[j] = alpha_j; - sum_alpha_ += delta_alpha_i + delta_alpha_j; - - // this is only necessary when i or j are not bound, but there is nothing - // wrong with doing this all the time - error_[i] = 0; - error_[j] = 0; - - DEBUG_GOT_HERE(0); - return true; -} - -#endif diff --git a/fastlib/u/garryb/svm/svm.cc b/fastlib/u/garryb/svm/svm.cc deleted file mode 100644 index c109b254b4..0000000000 --- a/fastlib/u/garryb/svm/svm.cc +++ /dev/null @@ -1,130 +0,0 @@ -#include "svm.h" - -void DoSvmNormalize(Dataset* dataset) { - Matrix m; - Vector sums; - - m.Init(dataset->n_features()-1, dataset->n_points()); - sums.Init(dataset->n_features() - 1); - sums.SetZero(); - - for (index_t i = 0; i < dataset->n_points(); i++) { - Vector s; - Vector d; - dataset->matrix().MakeColumnSubvector(i, 0, dataset->n_features()-1, &s); - m.MakeColumnVector(i, &d); - d.CopyValues(s); - la::AddTo(s, &sums); - } - - la::Scale(-1.0 / dataset->n_points(), &sums); - for (index_t i = 0; i < dataset->n_points(); i++) { - Vector d; - m.MakeColumnVector(i, &d); - la::AddTo(sums, &d); - } - - Matrix cov; - - - la::MulTransBInit(m, m, &cov); - - Vector d; - Matrix u; // eigenvectors - Matrix ui; // the inverse of eigenvectors - - //cov.PrintDebug("cov"); - la::EigenvectorsInit(cov, &d, &u); - la::TransposeInit(u, &ui); - - for (index_t i = 0; i < d.length(); i++) { - d[i] = 1.0 / sqrt(d[i] / (dataset->n_points() - 1)); - } - - la::ScaleRows(d, &ui); - - Matrix cov_inv_half; - la::MulInit(u, ui, &cov_inv_half); - - Matrix final; - la::MulInit(cov_inv_half, m, &final); - - for (index_t i = 0; i < dataset->n_points(); i++) { - Vector s; - Vector d; - dataset->matrix().MakeColumnSubvector(i, 0, dataset->n_features()-1, &d); - final.MakeColumnVector(i, &s); - d.CopyValues(s); - } - - //dataset->matrix().PrintDebug("m"); -} - -int main(int argc, char *argv[]) { - fx_init(argc, argv); - - 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" - - 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); - } - - 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 deleted file mode 100644 index 8b63bc9573..0000000000 --- a/fastlib/u/garryb/svm/svm.h +++ /dev/null @@ -1,113 +0,0 @@ -#ifndef U_SVM_SVM_H -#define U_SVM_SVM_H - -#include "smo.h" - -#include "fastlib/fastlib.h" - -#include - -struct SVMLinearKernel { - void Init(datanode *node) {} - - void Copy(const SVMLinearKernel& other) {} - - double Eval(const Vector& a, const Vector& b) const { - return la::Dot(a, b); - } -}; - -class SVMRBFKernel { - private: - double sigma_; - double gamma_; - - public: - - void Init(datanode *node) { - sigma_ = fx_param_double_req(NULL, "sigma"); - gamma_ = -1.0 / (2 * math::Sqr(sigma_)); - } - - 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(gamma_ * distance_squared); - } - - double sigma() const { - return sigma_; - } -}; - -template -class SVM { - public: - typedef TKernel Kernel; - - 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); - - int Classify(const Vector& vector); -}; - -template -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_)); - - fx_format_result(module, "n_support", "%"LI"d", alpha_.length()); -} - -template -int SVM::Classify(const Vector& datum) { - 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; - } - - DEBUG_MSG(0, "summation=%f, thresh_=%f", summation, thresh_); - - return (summation - thresh_ > 0.0) ? 1 : 0; -} - -#endif