This is an automated email from the git hooks/post-receive script. tille pushed a commit to branch master in repository psortb.
commit d133b9eb470ef3ff131ad2291578dfbe904f0128 Author: Carnë Draug <[email protected]> Date: Thu Apr 20 20:11:28 2017 +0100 Patch to remove bundled libsvm from bio-tools-psort-svmloc. --- debian/control | 1 + debian/patches/remove-bundled-party-libsvm.patch | 3147 ++++++++++++++++++++++ debian/patches/series | 1 + 3 files changed, 3149 insertions(+) diff --git a/debian/control b/debian/control index 3cd6911..f8c0468 100644 --- a/debian/control +++ b/debian/control @@ -4,6 +4,7 @@ Uploaders: Andreas Tille <[email protected]> Section: science Priority: optional Build-Depends: debhelper (>= 10), + libsvm-dev, bioperl, ncbi-blast+-legacy, libsquid-dev (>= 1.9g+cvs20050121-9~), diff --git a/debian/patches/remove-bundled-party-libsvm.patch b/debian/patches/remove-bundled-party-libsvm.patch new file mode 100644 index 0000000..960ba43 --- /dev/null +++ b/debian/patches/remove-bundled-party-libsvm.patch @@ -0,0 +1,3147 @@ +Description: Remove source from bundled libsvm and adjust build system. +Author: Carnë Draug <[email protected]> +Last-Update: 2017-04-20 +--- a/bio-tools-psort-svmloc/MANIFEST ++++ b/bio-tools-psort-svmloc/MANIFEST +@@ -8,5 +8,4 @@ + lib/Bio/Tools/PSort/SVMLoc.pm + sample.model + fre_patterns.txt +-libsvm.h + typemap +--- a/bio-tools-psort-svmloc/Makefile.PL ++++ b/bio-tools-psort-svmloc/Makefile.PL +@@ -9,7 +9,7 @@ + + @libs = qw/svmloc/; + %paths = (); +-$libpath = '-lm -lstdc++'; ++$libpath = '-lm -lsvm -lstdc++'; + + $defaults_file = '../psortb.defaults'; + if( -f $defaults_file ) { +--- a/bio-tools-psort-svmloc/bindings.h ++++ b/bio-tools-psort-svmloc/bindings.h +@@ -9,7 +9,9 @@ + #include <set> + #include <assert.h> + +-#include "libsvm.h" ++#include <stdlib.h> ++#include <string.h> ++#include <libsvm/svm.h> + + using namespace std; + +--- a/bio-tools-psort-svmloc/libsvm.cpp ++++ /dev/null +@@ -1,3036 +0,0 @@ +-#include <math.h> +-#include <stdio.h> +-#include <stdlib.h> +-#include <ctype.h> +-#include <float.h> +-#include <string.h> +-#include <stdarg.h> +-#include "libsvm.h" +-typedef float Qfloat; +-typedef signed char schar; +-#ifndef min +-template <class T> inline T min(T x,T y) { return (x<y)?x:y; } +-#endif +-#ifndef max +-template <class T> inline T max(T x,T y) { return (x>y)?x:y; } +-#endif +-template <class T> inline void swap(T& x, T& y) { T t=x; x=y; y=t; } +-template <class S, class T> inline void clone(T*& dst, S* src, int n) +-{ +- dst = new T[n]; +- memcpy((void *)dst,(void *)src,sizeof(T)*n); +-} +-#define INF HUGE_VAL +-#define TAU 1e-12 +-#define Malloc(type,n) (type *)malloc((n)*sizeof(type)) +-#if 1 +-void info(char *fmt,...) +-{ +- va_list ap; +- va_start(ap,fmt); +- vprintf(fmt,ap); +- va_end(ap); +-} +-void info_flush() +-{ +- fflush(stdout); +-} +-#else +-void info(char *fmt,...) {} +-void info_flush() {} +-#endif +- +-// +-// Kernel Cache +-// +-// l is the number of total data items +-// size is the cache size limit in bytes +-// +-class Cache +-{ +-public: +- Cache(int l,int size); +- ~Cache(); +- +- // request data [0,len) +- // return some position p where [p,len) need to be filled +- // (p >= len if nothing needs to be filled) +- int get_data(const int index, Qfloat **data, int len); +- void swap_index(int i, int j); // future_option +-private: +- int l; +- int size; +- struct head_t +- { +- head_t *prev, *next; // a cicular list +- Qfloat *data; +- int len; // data[0,len) is cached in this entry +- }; +- +- head_t *head; +- head_t lru_head; +- void lru_delete(head_t *h); +- void lru_insert(head_t *h); +-}; +- +-Cache::Cache(int l_,int size_):l(l_),size(size_) +-{ +- head = (head_t *)calloc(l,sizeof(head_t)); // initialized to 0 +- size /= sizeof(Qfloat); +- size -= l * sizeof(head_t) / sizeof(Qfloat); +- size = max(size, 2*l); // cache must be large enough for two columns +- lru_head.next = lru_head.prev = &lru_head; +-} +- +-Cache::~Cache() +-{ +- for(head_t *h = lru_head.next; h != &lru_head; h=h->next) +- free(h->data); +- free(head); +-} +- +-void Cache::lru_delete(head_t *h) +-{ +- // delete from current location +- h->prev->next = h->next; +- h->next->prev = h->prev; +-} +- +-void Cache::lru_insert(head_t *h) +-{ +- // insert to last position +- h->next = &lru_head; +- h->prev = lru_head.prev; +- h->prev->next = h; +- h->next->prev = h; +-} +- +-int Cache::get_data(const int index, Qfloat **data, int len) +-{ +- head_t *h = &head[index]; +- if(h->len) lru_delete(h); +- int more = len - h->len; +- +- if(more > 0) +- { +- // free old space +- while(size < more) +- { +- head_t *old = lru_head.next; +- lru_delete(old); +- free(old->data); +- size += old->len; +- old->data = 0; +- old->len = 0; +- } +- +- // allocate new space +- h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len); +- size -= more; +- swap(h->len,len); +- } +- +- lru_insert(h); +- *data = h->data; +- return len; +-} +- +-void Cache::swap_index(int i, int j) +-{ +- if(i==j) return; +- +- if(head[i].len) lru_delete(&head[i]); +- if(head[j].len) lru_delete(&head[j]); +- swap(head[i].data,head[j].data); +- swap(head[i].len,head[j].len); +- if(head[i].len) lru_insert(&head[i]); +- if(head[j].len) lru_insert(&head[j]); +- +- if(i>j) swap(i,j); +- for(head_t *h = lru_head.next; h!=&lru_head; h=h->next) +- { +- if(h->len > i) +- { +- if(h->len > j) +- swap(h->data[i],h->data[j]); +- else +- { +- // give up +- lru_delete(h); +- free(h->data); +- size += h->len; +- h->data = 0; +- h->len = 0; +- } +- } +- } +-} +- +-// +-// Kernel evaluation +-// +-// the static method k_function is for doing single kernel evaluation +-// the constructor of Kernel prepares to calculate the l*l kernel matrix +-// the member function get_Q is for getting one column from the Q Matrix +-// +-class QMatrix { +-public: +- virtual Qfloat *get_Q(int column, int len) const = 0; +- virtual Qfloat *get_QD() const = 0; +- virtual void swap_index(int i, int j) const = 0; +-}; +- +-class Kernel: public QMatrix { +-public: +- Kernel(int l, svm_node * const * x, const svm_parameter& param); +- virtual ~Kernel(); +- +- static double k_function(const svm_node *x, const svm_node *y, +- const svm_parameter& param); +- virtual Qfloat *get_Q(int column, int len) const = 0; +- virtual Qfloat *get_QD() const = 0; +- virtual void swap_index(int i, int j) const // no so const... +- { +- swap(x[i],x[j]); +- if(x_square) swap(x_square[i],x_square[j]); +- } +-protected: +- +- double (Kernel::*kernel_function)(int i, int j) const; +- +-private: +- const svm_node **x; +- double *x_square; +- +- // svm_parameter +- const int kernel_type; +- const double degree; +- const double gamma; +- const double coef0; +- +- static double dot(const svm_node *px, const svm_node *py); +- double kernel_linear(int i, int j) const +- { +- return dot(x[i],x[j]); +- } +- double kernel_poly(int i, int j) const +- { +- return pow(gamma*dot(x[i],x[j])+coef0,degree); +- } +- double kernel_rbf(int i, int j) const +- { +- return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j]))); +- } +- double kernel_sigmoid(int i, int j) const +- { +- return tanh(gamma*dot(x[i],x[j])+coef0); +- } +-}; +- +-Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param) +-:kernel_type(param.kernel_type), degree(param.degree), +- gamma(param.gamma), coef0(param.coef0) +-{ +- switch(kernel_type) +- { +- case LINEAR: +- kernel_function = &Kernel::kernel_linear; +- break; +- case POLY: +- kernel_function = &Kernel::kernel_poly; +- break; +- case RBF: +- kernel_function = &Kernel::kernel_rbf; +- break; +- case SIGMOID: +- kernel_function = &Kernel::kernel_sigmoid; +- break; +- } +- +- clone(x,x_,l); +- +- if(kernel_type == RBF) +- { +- x_square = new double[l]; +- for(int i=0;i<l;i++) +- x_square[i] = dot(x[i],x[i]); +- } +- else +- x_square = 0; +-} +- +-Kernel::~Kernel() +-{ +- delete[] x; +- delete[] x_square; +-} +- +-double Kernel::dot(const svm_node *px, const svm_node *py) +-{ +- double sum = 0; +- while(px->index != -1 && py->index != -1) +- { +- if(px->index == py->index) +- { +- sum += px->value * py->value; +- ++px; +- ++py; +- } +- else +- { +- if(px->index > py->index) +- ++py; +- else +- ++px; +- } +- } +- return sum; +-} +- +-double Kernel::k_function(const svm_node *x, const svm_node *y, +- const svm_parameter& param) +-{ +- switch(param.kernel_type) +- { +- case LINEAR: +- return dot(x,y); +- case POLY: +- return pow(param.gamma*dot(x,y)+param.coef0,param.degree); +- case RBF: +- { +- double sum = 0; +- while(x->index != -1 && y->index !=-1) +- { +- if(x->index == y->index) +- { +- double d = x->value - y->value; +- sum += d*d; +- ++x; +- ++y; +- } +- else +- { +- if(x->index > y->index) +- { +- sum += y->value * y->value; +- ++y; +- } +- else +- { +- sum += x->value * x->value; +- ++x; +- } +- } +- } +- +- while(x->index != -1) +- { +- sum += x->value * x->value; +- ++x; +- } +- +- while(y->index != -1) +- { +- sum += y->value * y->value; +- ++y; +- } +- +- return exp(-param.gamma*sum); +- } +- case SIGMOID: +- return tanh(param.gamma*dot(x,y)+param.coef0); +- default: +- return 0; /* Unreachable */ +- } +-} +- +-// Generalized SMO+SVMlight algorithm +-// Solves: +-// +-// min 0.5(\alpha^T Q \alpha) + b^T \alpha +-// +-// y^T \alpha = \delta +-// y_i = +1 or -1 +-// 0 <= alpha_i <= Cp for y_i = 1 +-// 0 <= alpha_i <= Cn for y_i = -1 +-// +-// Given: +-// +-// Q, b, y, Cp, Cn, and an initial feasible point \alpha +-// l is the size of vectors and matrices +-// eps is the stopping criterion +-// +-// solution will be put in \alpha, objective value will be put in obj +-// +-class Solver { +-public: +- Solver() {}; +- virtual ~Solver() {}; +- +- struct SolutionInfo { +- double obj; +- double rho; +- double upper_bound_p; +- double upper_bound_n; +- double r; // for Solver_NU +- }; +- +- void Solve(int l, const QMatrix& Q, const double *b_, const schar *y_, +- double *alpha_, double Cp, double Cn, double eps, +- SolutionInfo* si, int shrinking); +-protected: +- int active_size; +- schar *y; +- double *G; // gradient of objective function +- enum { LOWER_BOUND, UPPER_BOUND, FREE }; +- char *alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE +- double *alpha; +- const QMatrix *Q; +- const Qfloat *QD; +- double eps; +- double Cp,Cn; +- double *b; +- int *active_set; +- double *G_bar; // gradient, if we treat free variables as 0 +- int l; +- bool unshrinked; // XXX +- +- double get_C(int i) +- { +- return (y[i] > 0)? Cp : Cn; +- } +- void update_alpha_status(int i) +- { +- if(alpha[i] >= get_C(i)) +- alpha_status[i] = UPPER_BOUND; +- else if(alpha[i] <= 0) +- alpha_status[i] = LOWER_BOUND; +- else alpha_status[i] = FREE; +- } +- bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; } +- bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; } +- bool is_free(int i) { return alpha_status[i] == FREE; } +- void swap_index(int i, int j); +- void reconstruct_gradient(); +- virtual int select_working_set(int &i, int &j); +- virtual int max_violating_pair(int &i, int &j); +- virtual double calculate_rho(); +- virtual void do_shrinking(); +-}; +- +-void Solver::swap_index(int i, int j) +-{ +- Q->swap_index(i,j); +- swap(y[i],y[j]); +- swap(G[i],G[j]); +- swap(alpha_status[i],alpha_status[j]); +- swap(alpha[i],alpha[j]); +- swap(b[i],b[j]); +- swap(active_set[i],active_set[j]); +- swap(G_bar[i],G_bar[j]); +-} +- +-void Solver::reconstruct_gradient() +-{ +- // reconstruct inactive elements of G from G_bar and free variables +- +- if(active_size == l) return; +- +- int i; +- for(i=active_size;i<l;i++) +- G[i] = G_bar[i] + b[i]; +- +- for(i=0;i<active_size;i++) +- if(is_free(i)) +- { +- const Qfloat *Q_i = Q->get_Q(i,l); +- double alpha_i = alpha[i]; +- for(int j=active_size;j<l;j++) +- G[j] += alpha_i * Q_i[j]; +- } +-} +- +-void Solver::Solve(int l, const QMatrix& Q, const double *b_, const schar *y_, +- double *alpha_, double Cp, double Cn, double eps, +- SolutionInfo* si, int shrinking) +-{ +- this->l = l; +- this->Q = &Q; +- QD=Q.get_QD(); +- clone(b, b_,l); +- clone(y, y_,l); +- clone(alpha,alpha_,l); +- this->Cp = Cp; +- this->Cn = Cn; +- this->eps = eps; +- unshrinked = false; +- +- // initialize alpha_status +- { +- alpha_status = new char[l]; +- for(int i=0;i<l;i++) +- update_alpha_status(i); +- } +- +- // initialize active set (for shrinking) +- { +- active_set = new int[l]; +- for(int i=0;i<l;i++) +- active_set[i] = i; +- active_size = l; +- } +- +- // initialize gradient +- { +- G = new double[l]; +- G_bar = new double[l]; +- int i; +- for(i=0;i<l;i++) +- { +- G[i] = b[i]; +- G_bar[i] = 0; +- } +- for(i=0;i<l;i++) +- if(!is_lower_bound(i)) +- { +- const Qfloat *Q_i = Q.get_Q(i,l); +- double alpha_i = alpha[i]; +- int j; +- for(j=0;j<l;j++) +- G[j] += alpha_i*Q_i[j]; +- if(is_upper_bound(i)) +- for(j=0;j<l;j++) +- G_bar[j] += get_C(i) * Q_i[j]; +- } +- } +- +- // optimization step +- +- int iter = 0; +- int counter = min(l,1000)+1; +- +- while(1) +- { +- // show progress and do shrinking +- +- if(--counter == 0) +- { +- counter = min(l,1000); +- if(shrinking) do_shrinking(); +- info("."); info_flush(); +- } +- +- int i,j; +- if(select_working_set(i,j)!=0) +- { +- // reconstruct the whole gradient +- reconstruct_gradient(); +- // reset active set size and check +- active_size = l; +- info("*"); info_flush(); +- if(select_working_set(i,j)!=0) +- break; +- else +- counter = 1; // do shrinking next iteration +- } +- +- ++iter; +- +- // update alpha[i] and alpha[j], handle bounds carefully +- +- const Qfloat *Q_i = Q.get_Q(i,active_size); +- const Qfloat *Q_j = Q.get_Q(j,active_size); +- +- double C_i = get_C(i); +- double C_j = get_C(j); +- +- double old_alpha_i = alpha[i]; +- double old_alpha_j = alpha[j]; +- +- if(y[i]!=y[j]) +- { +- double quad_coef = Q_i[i]+Q_j[j]+2*Q_i[j]; +- if (quad_coef <= 0) +- quad_coef = TAU; +- double delta = (-G[i]-G[j])/quad_coef; +- double diff = alpha[i] - alpha[j]; +- alpha[i] += delta; +- alpha[j] += delta; +- +- if(diff > 0) +- { +- if(alpha[j] < 0) +- { +- alpha[j] = 0; +- alpha[i] = diff; +- } +- } +- else +- { +- if(alpha[i] < 0) +- { +- alpha[i] = 0; +- alpha[j] = -diff; +- } +- } +- if(diff > C_i - C_j) +- { +- if(alpha[i] > C_i) +- { +- alpha[i] = C_i; +- alpha[j] = C_i - diff; +- } +- } +- else +- { +- if(alpha[j] > C_j) +- { +- alpha[j] = C_j; +- alpha[i] = C_j + diff; +- } +- } +- } +- else +- { +- double quad_coef = Q_i[i]+Q_j[j]-2*Q_i[j]; +- if (quad_coef <= 0) +- quad_coef = TAU; +- double delta = (G[i]-G[j])/quad_coef; +- double sum = alpha[i] + alpha[j]; +- alpha[i] -= delta; +- alpha[j] += delta; +- +- if(sum > C_i) +- { +- if(alpha[i] > C_i) +- { +- alpha[i] = C_i; +- alpha[j] = sum - C_i; +- } +- } +- else +- { +- if(alpha[j] < 0) +- { +- alpha[j] = 0; +- alpha[i] = sum; +- } +- } +- if(sum > C_j) +- { +- if(alpha[j] > C_j) +- { +- alpha[j] = C_j; +- alpha[i] = sum - C_j; +- } +- } +- else +- { +- if(alpha[i] < 0) +- { +- alpha[i] = 0; +- alpha[j] = sum; +- } +- } +- } +- +- // update G +- +- double delta_alpha_i = alpha[i] - old_alpha_i; +- double delta_alpha_j = alpha[j] - old_alpha_j; +- +- for(int k=0;k<active_size;k++) +- { +- G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j; +- } +- +- // update alpha_status and G_bar +- +- { +- bool ui = is_upper_bound(i); +- bool uj = is_upper_bound(j); +- update_alpha_status(i); +- update_alpha_status(j); +- int k; +- if(ui != is_upper_bound(i)) +- { +- Q_i = Q.get_Q(i,l); +- if(ui) +- for(k=0;k<l;k++) +- G_bar[k] -= C_i * Q_i[k]; +- else +- for(k=0;k<l;k++) +- G_bar[k] += C_i * Q_i[k]; +- } +- +- if(uj != is_upper_bound(j)) +- { +- Q_j = Q.get_Q(j,l); +- if(uj) +- for(k=0;k<l;k++) +- G_bar[k] -= C_j * Q_j[k]; +- else +- for(k=0;k<l;k++) +- G_bar[k] += C_j * Q_j[k]; +- } +- } +- } +- +- // calculate rho +- +- si->rho = calculate_rho(); +- +- // calculate objective value +- { +- double v = 0; +- int i; +- for(i=0;i<l;i++) +- v += alpha[i] * (G[i] + b[i]); +- +- si->obj = v/2; +- } +- +- // put back the solution +- { +- for(int i=0;i<l;i++) +- alpha_[active_set[i]] = alpha[i]; +- } +- +- // juggle everything back +- /*{ +- for(int i=0;i<l;i++) +- while(active_set[i] != i) +- swap_index(i,active_set[i]); +- // or Q.swap_index(i,active_set[i]); +- }*/ +- +- si->upper_bound_p = Cp; +- si->upper_bound_n = Cn; +- +- info("\noptimization finished, #iter = %d\n",iter); +- +- delete[] b; +- delete[] y; +- delete[] alpha; +- delete[] alpha_status; +- delete[] active_set; +- delete[] G; +- delete[] G_bar; +-} +- +-// return 1 if already optimal, return 0 otherwise +-int Solver::select_working_set(int &out_i, int &out_j) +-{ +- // return i,j such that +- // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha) +- // j: minimizes the decrease of obj value +- // (if quadratic coefficeint <= 0, replace it with tau) +- // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha) +- +- double Gmax = -INF; +- int Gmax_idx = -1; +- int Gmin_idx = -1; +- double obj_diff_min = INF; +- +- for(int t=0;t<active_size;t++) +- if(y[t]==+1) +- { +- if(!is_upper_bound(t)) +- if(-G[t] >= Gmax) +- { +- Gmax = -G[t]; +- Gmax_idx = t; +- } +- } +- else +- { +- if(!is_lower_bound(t)) +- if(G[t] >= Gmax) +- { +- Gmax = G[t]; +- Gmax_idx = t; +- } +- } +- +- int i = Gmax_idx; +- const Qfloat *Q_i = NULL; +- if(i != -1) // NULL Q_i not accessed: Gmax=-INF if i=-1 +- Q_i = Q->get_Q(i,active_size); +- +- for(int j=0;j<active_size;j++) +- { +- if(y[j]==+1) +- { +- if (!is_lower_bound(j)) +- { +- double grad_diff=Gmax+G[j]; +- if (grad_diff >= eps) +- { +- double obj_diff; +- double quad_coef=Q_i[i]+QD[j]-2*y[i]*Q_i[j]; +- if (quad_coef > 0) +- obj_diff = -(grad_diff*grad_diff)/quad_coef; +- else +- obj_diff = -(grad_diff*grad_diff)/TAU; +- +- if (obj_diff <= obj_diff_min) +- { +- Gmin_idx=j; +- obj_diff_min = obj_diff; +- } +- } +- } +- } +- else +- { +- if (!is_upper_bound(j)) +- { +- double grad_diff= Gmax-G[j]; +- if (grad_diff >= eps) +- { +- double obj_diff; +- double quad_coef=Q_i[i]+QD[j]+2*y[i]*Q_i[j]; +- if (quad_coef > 0) +- obj_diff = -(grad_diff*grad_diff)/quad_coef; +- else +- obj_diff = -(grad_diff*grad_diff)/TAU; +- +- if (obj_diff <= obj_diff_min) +- { +- Gmin_idx=j; +- obj_diff_min = obj_diff; +- } +- } +- } +- } +- } +- +- if(Gmin_idx == -1) +- return 1; +- +- out_i = Gmax_idx; +- out_j = Gmin_idx; +- return 0; +-} +- +-// return 1 if already optimal, return 0 otherwise +-int Solver::max_violating_pair(int &out_i, int &out_j) +-{ +- // return i,j: maximal violating pair +- +- double Gmax1 = -INF; // max { -y_i * grad(f)_i | i in I_up(\alpha) } +- int Gmax1_idx = -1; +- +- double Gmax2 = -INF; // max { y_i * grad(f)_i | i in I_low(\alpha) } +- int Gmax2_idx = -1; +- +- for(int i=0;i<active_size;i++) +- { +- if(y[i]==+1) // y = +1 +- { +- if(!is_upper_bound(i)) // d = +1 +- { +- if(-G[i] >= Gmax1) +- { +- Gmax1 = -G[i]; +- Gmax1_idx = i; +- } +- } +- if(!is_lower_bound(i)) // d = -1 +- { +- if(G[i] >= Gmax2) +- { +- Gmax2 = G[i]; +- Gmax2_idx = i; +- } +- } +- } +- else // y = -1 +- { +- if(!is_upper_bound(i)) // d = +1 +- { +- if(-G[i] >= Gmax2) +- { +- Gmax2 = -G[i]; +- Gmax2_idx = i; +- } +- } +- if(!is_lower_bound(i)) // d = -1 +- { +- if(G[i] >= Gmax1) +- { +- Gmax1 = G[i]; +- Gmax1_idx = i; +- } +- } +- } +- } +- +- if(Gmax1+Gmax2 < eps) +- return 1; +- +- out_i = Gmax1_idx; +- out_j = Gmax2_idx; +- return 0; +-} +- +-void Solver::do_shrinking() +-{ +- int i,j,k; +- if(max_violating_pair(i,j)!=0) return; +- double Gm1 = -y[j]*G[j]; +- double Gm2 = y[i]*G[i]; +- +- // shrink +- +- for(k=0;k<active_size;k++) +- { +- if(is_lower_bound(k)) +- { +- if(y[k]==+1) +- { +- if(-G[k] >= Gm1) continue; +- } +- else if(-G[k] >= Gm2) continue; +- } +- else if(is_upper_bound(k)) +- { +- if(y[k]==+1) +- { +- if(G[k] >= Gm2) continue; +- } +- else if(G[k] >= Gm1) continue; +- } +- else continue; +- +- --active_size; +- swap_index(k,active_size); +- --k; // look at the newcomer +- } +- +- // unshrink, check all variables again before final iterations +- +- if(unshrinked || -(Gm1 + Gm2) > eps*10) return; +- +- unshrinked = true; +- reconstruct_gradient(); +- +- for(k=l-1;k>=active_size;k--) +- { +- if(is_lower_bound(k)) +- { +- if(y[k]==+1) +- { +- if(-G[k] < Gm1) continue; +- } +- else if(-G[k] < Gm2) continue; +- } +- else if(is_upper_bound(k)) +- { +- if(y[k]==+1) +- { +- if(G[k] < Gm2) continue; +- } +- else if(G[k] < Gm1) continue; +- } +- else continue; +- +- swap_index(k,active_size); +- active_size++; +- ++k; // look at the newcomer +- } +-} +- +-double Solver::calculate_rho() +-{ +- double r; +- int nr_free = 0; +- double ub = INF, lb = -INF, sum_free = 0; +- for(int i=0;i<active_size;i++) +- { +- double yG = y[i]*G[i]; +- +- if(is_lower_bound(i)) +- { +- if(y[i] > 0) +- ub = min(ub,yG); +- else +- lb = max(lb,yG); +- } +- else if(is_upper_bound(i)) +- { +- if(y[i] < 0) +- ub = min(ub,yG); +- else +- lb = max(lb,yG); +- } +- else +- { +- ++nr_free; +- sum_free += yG; +- } +- } +- +- if(nr_free>0) +- r = sum_free/nr_free; +- else +- r = (ub+lb)/2; +- +- return r; +-} +- +-// +-// Solver for nu-svm classification and regression +-// +-// additional constraint: e^T \alpha = constant +-// +-class Solver_NU : public Solver +-{ +-public: +- Solver_NU() {} +- void Solve(int l, const QMatrix& Q, const double *b, const schar *y, +- double *alpha, double Cp, double Cn, double eps, +- SolutionInfo* si, int shrinking) +- { +- this->si = si; +- Solver::Solve(l,Q,b,y,alpha,Cp,Cn,eps,si,shrinking); +- } +-private: +- SolutionInfo *si; +- int select_working_set(int &i, int &j); +- double calculate_rho(); +- void do_shrinking(); +-}; +- +-// return 1 if already optimal, return 0 otherwise +-int Solver_NU::select_working_set(int &out_i, int &out_j) +-{ +- // return i,j such that y_i = y_j and +- // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha) +- // j: minimizes the decrease of obj value +- // (if quadratic coefficeint <= 0, replace it with tau) +- // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha) +- +- double Gmaxp = -INF; +- int Gmaxp_idx = -1; +- +- double Gmaxn = -INF; +- int Gmaxn_idx = -1; +- +- int Gmin_idx = -1; +- double obj_diff_min = INF; +- +- for(int t=0;t<active_size;t++) +- if(y[t]==+1) +- { +- if(!is_upper_bound(t)) +- if(-G[t] >= Gmaxp) +- { +- Gmaxp = -G[t]; +- Gmaxp_idx = t; +- } +- } +- else +- { +- if(!is_lower_bound(t)) +- if(G[t] >= Gmaxn) +- { +- Gmaxn = G[t]; +- Gmaxn_idx = t; +- } +- } +- +- int ip = Gmaxp_idx; +- int in = Gmaxn_idx; +- const Qfloat *Q_ip = NULL; +- const Qfloat *Q_in = NULL; +- if(ip != -1) // NULL Q_ip not accessed: Gmaxp=-INF if ip=-1 +- Q_ip = Q->get_Q(ip,active_size); +- if(in != -1) +- Q_in = Q->get_Q(in,active_size); +- +- for(int j=0;j<active_size;j++) +- { +- if(y[j]==+1) +- { +- if (!is_lower_bound(j)) +- { +- double grad_diff=Gmaxp+G[j]; +- if (grad_diff >= eps) +- { +- double obj_diff; +- double quad_coef = Q_ip[ip]+QD[j]-2*Q_ip[j]; +- if (quad_coef > 0) +- obj_diff = -(grad_diff*grad_diff)/quad_coef; +- else +- obj_diff = -(grad_diff*grad_diff)/TAU; +- +- if (obj_diff <= obj_diff_min) +- { +- Gmin_idx=j; +- obj_diff_min = obj_diff; +- } +- } +- } +- } +- else +- { +- if (!is_upper_bound(j)) +- { +- double grad_diff=Gmaxn-G[j]; +- if (grad_diff >= eps) +- { +- double obj_diff; +- double quad_coef = Q_in[in]+QD[j]-2*Q_in[j]; +- if (quad_coef > 0) +- obj_diff = -(grad_diff*grad_diff)/quad_coef; +- else +- obj_diff = -(grad_diff*grad_diff)/TAU; +- +- if (obj_diff <= obj_diff_min) +- { +- Gmin_idx=j; +- obj_diff_min = obj_diff; +- } +- } +- } +- } +- } +- +- if(Gmin_idx == -1) +- return 1; +- +- if (y[Gmin_idx] == +1) +- out_i = Gmaxp_idx; +- else +- out_i = Gmaxn_idx; +- out_j = Gmin_idx; +- +- return 0; +-} +- +-void Solver_NU::do_shrinking() +-{ +- double Gmax1 = -INF; // max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) } +- double Gmax2 = -INF; // max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) } +- double Gmax3 = -INF; // max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) } +- double Gmax4 = -INF; // max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) } +- +- // find maximal violating pair first +- int k; +- for(k=0;k<active_size;k++) +- { +- if(!is_upper_bound(k)) +- { +- if(y[k]==+1) +- { +- if(-G[k] > Gmax1) Gmax1 = -G[k]; +- } +- else if(-G[k] > Gmax3) Gmax3 = -G[k]; +- } +- if(!is_lower_bound(k)) +- { +- if(y[k]==+1) +- { +- if(G[k] > Gmax2) Gmax2 = G[k]; +- } +- else if(G[k] > Gmax4) Gmax4 = G[k]; +- } +- } +- +- // shrinking +- +- double Gm1 = -Gmax2; +- double Gm2 = -Gmax1; +- double Gm3 = -Gmax4; +- double Gm4 = -Gmax3; +- +- for(k=0;k<active_size;k++) +- { +- if(is_lower_bound(k)) +- { +- if(y[k]==+1) +- { +- if(-G[k] >= Gm1) continue; +- } +- else if(-G[k] >= Gm3) continue; +- } +- else if(is_upper_bound(k)) +- { +- if(y[k]==+1) +- { +- if(G[k] >= Gm2) continue; +- } +- else if(G[k] >= Gm4) continue; +- } +- else continue; +- +- --active_size; +- swap_index(k,active_size); +- --k; // look at the newcomer +- } +- +- // unshrink, check all variables again before final iterations +- +- if(unshrinked || max(-(Gm1+Gm2),-(Gm3+Gm4)) > eps*10) return; +- +- unshrinked = true; +- reconstruct_gradient(); +- +- for(k=l-1;k>=active_size;k--) +- { +- if(is_lower_bound(k)) +- { +- if(y[k]==+1) +- { +- if(-G[k] < Gm1) continue; +- } +- else if(-G[k] < Gm3) continue; +- } +- else if(is_upper_bound(k)) +- { +- if(y[k]==+1) +- { +- if(G[k] < Gm2) continue; +- } +- else if(G[k] < Gm4) continue; +- } +- else continue; +- +- swap_index(k,active_size); +- active_size++; +- ++k; // look at the newcomer +- } +-} +- +-double Solver_NU::calculate_rho() +-{ +- int nr_free1 = 0,nr_free2 = 0; +- double ub1 = INF, ub2 = INF; +- double lb1 = -INF, lb2 = -INF; +- double sum_free1 = 0, sum_free2 = 0; +- +- for(int i=0;i<active_size;i++) +- { +- if(y[i]==+1) +- { +- if(is_lower_bound(i)) +- ub1 = min(ub1,G[i]); +- else if(is_upper_bound(i)) +- lb1 = max(lb1,G[i]); +- else +- { +- ++nr_free1; +- sum_free1 += G[i]; +- } +- } +- else +- { +- if(is_lower_bound(i)) +- ub2 = min(ub2,G[i]); +- else if(is_upper_bound(i)) +- lb2 = max(lb2,G[i]); +- else +- { +- ++nr_free2; +- sum_free2 += G[i]; +- } +- } +- } +- +- double r1,r2; +- if(nr_free1 > 0) +- r1 = sum_free1/nr_free1; +- else +- r1 = (ub1+lb1)/2; +- +- if(nr_free2 > 0) +- r2 = sum_free2/nr_free2; +- else +- r2 = (ub2+lb2)/2; +- +- si->r = (r1+r2)/2; +- return (r1-r2)/2; +-} +- +-// +-// Q matrices for various formulations +-// +-class SVC_Q: public Kernel +-{ +-public: +- SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_) +- :Kernel(prob.l, prob.x, param) +- { +- clone(y,y_,prob.l); +- cache = new Cache(prob.l,(int)(param.cache_size*(1<<20))); +- QD = new Qfloat[prob.l]; +- for(int i=0;i<prob.l;i++) +- QD[i]= (Qfloat)(this->*kernel_function)(i,i); +- } +- +- Qfloat *get_Q(int i, int len) const +- { +- Qfloat *data; +- int start; +- if((start = cache->get_data(i,&data,len)) < len) +- { +- for(int j=start;j<len;j++) +- data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j)); +- } +- return data; +- } +- +- Qfloat *get_QD() const +- { +- return QD; +- } +- +- void swap_index(int i, int j) const +- { +- cache->swap_index(i,j); +- Kernel::swap_index(i,j); +- swap(y[i],y[j]); +- swap(QD[i],QD[j]); +- } +- +- ~SVC_Q() +- { +- delete[] y; +- delete cache; +- delete[] QD; +- } +-private: +- schar *y; +- Cache *cache; +- Qfloat *QD; +-}; +- +-class ONE_CLASS_Q: public Kernel +-{ +-public: +- ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param) +- :Kernel(prob.l, prob.x, param) +- { +- cache = new Cache(prob.l,(int)(param.cache_size*(1<<20))); +- QD = new Qfloat[prob.l]; +- for(int i=0;i<prob.l;i++) +- QD[i]= (Qfloat)(this->*kernel_function)(i,i); +- } +- +- Qfloat *get_Q(int i, int len) const +- { +- Qfloat *data; +- int start; +- if((start = cache->get_data(i,&data,len)) < len) +- { +- for(int j=start;j<len;j++) +- data[j] = (Qfloat)(this->*kernel_function)(i,j); +- } +- return data; +- } +- +- Qfloat *get_QD() const +- { +- return QD; +- } +- +- void swap_index(int i, int j) const +- { +- cache->swap_index(i,j); +- Kernel::swap_index(i,j); +- swap(QD[i],QD[j]); +- } +- +- ~ONE_CLASS_Q() +- { +- delete cache; +- delete[] QD; +- } +-private: +- Cache *cache; +- Qfloat *QD; +-}; +- +-class SVR_Q: public Kernel +-{ +-public: +- SVR_Q(const svm_problem& prob, const svm_parameter& param) +- :Kernel(prob.l, prob.x, param) +- { +- l = prob.l; +- cache = new Cache(l,(int)(param.cache_size*(1<<20))); +- QD = new Qfloat[2*l]; +- sign = new schar[2*l]; +- index = new int[2*l]; +- for(int k=0;k<l;k++) +- { +- sign[k] = 1; +- sign[k+l] = -1; +- index[k] = k; +- index[k+l] = k; +- QD[k]= (Qfloat)(this->*kernel_function)(k,k); +- QD[k+l]=QD[k]; +- } +- buffer[0] = new Qfloat[2*l]; +- buffer[1] = new Qfloat[2*l]; +- next_buffer = 0; +- } +- +- void swap_index(int i, int j) const +- { +- swap(sign[i],sign[j]); +- swap(index[i],index[j]); +- swap(QD[i],QD[j]); +- } +- +- Qfloat *get_Q(int i, int len) const +- { +- Qfloat *data; +- int real_i = index[i]; +- if(cache->get_data(real_i,&data,l) < l) +- { +- for(int j=0;j<l;j++) +- data[j] = (Qfloat)(this->*kernel_function)(real_i,j); +- } +- +- // reorder and copy +- Qfloat *buf = buffer[next_buffer]; +- next_buffer = 1 - next_buffer; +- schar si = sign[i]; +- for(int j=0;j<len;j++) +- buf[j] = si * sign[j] * data[index[j]]; +- return buf; +- } +- +- Qfloat *get_QD() const +- { +- return QD; +- } +- +- ~SVR_Q() +- { +- delete cache; +- delete[] sign; +- delete[] index; +- delete[] buffer[0]; +- delete[] buffer[1]; +- delete[] QD; +- } +-private: +- int l; +- Cache *cache; +- schar *sign; +- int *index; +- mutable int next_buffer; +- Qfloat *buffer[2]; +- Qfloat *QD; +-}; +- +-// +-// construct and solve various formulations +-// +-static void solve_c_svc( +- const svm_problem *prob, const svm_parameter* param, +- double *alpha, Solver::SolutionInfo* si, double Cp, double Cn) +-{ +- int l = prob->l; +- double *minus_ones = new double[l]; +- schar *y = new schar[l]; +- +- int i; +- +- for(i=0;i<l;i++) +- { +- alpha[i] = 0; +- minus_ones[i] = -1; +- if(prob->y[i] > 0) y[i] = +1; else y[i]=-1; +- } +- +- Solver s; +- s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y, +- alpha, Cp, Cn, param->eps, si, param->shrinking); +- +- double sum_alpha=0; +- for(i=0;i<l;i++) +- sum_alpha += alpha[i]; +- +- if (Cp==Cn) +- info("nu = %f\n", sum_alpha/(Cp*prob->l)); +- +- for(i=0;i<l;i++) +- alpha[i] *= y[i]; +- +- delete[] minus_ones; +- delete[] y; +-} +- +-static void solve_nu_svc( +- const svm_problem *prob, const svm_parameter *param, +- double *alpha, Solver::SolutionInfo* si) +-{ +- int i; +- int l = prob->l; +- double nu = param->nu; +- +- schar *y = new schar[l]; +- +- for(i=0;i<l;i++) +- if(prob->y[i]>0) +- y[i] = +1; +- else +- y[i] = -1; +- +- double sum_pos = nu*l/2; +- double sum_neg = nu*l/2; +- +- for(i=0;i<l;i++) +- if(y[i] == +1) +- { +- alpha[i] = min(1.0,sum_pos); +- sum_pos -= alpha[i]; +- } +- else +- { +- alpha[i] = min(1.0,sum_neg); +- sum_neg -= alpha[i]; +- } +- +- double *zeros = new double[l]; +- +- for(i=0;i<l;i++) +- zeros[i] = 0; +- +- Solver_NU s; +- s.Solve(l, SVC_Q(*prob,*param,y), zeros, y, +- alpha, 1.0, 1.0, param->eps, si, param->shrinking); +- double r = si->r; +- +- info("C = %f\n",1/r); +- +- for(i=0;i<l;i++) +- alpha[i] *= y[i]/r; +- +- si->rho /= r; +- si->obj /= (r*r); +- si->upper_bound_p = 1/r; +- si->upper_bound_n = 1/r; +- +- delete[] y; +- delete[] zeros; +-} +- +-static void solve_one_class( +- const svm_problem *prob, const svm_parameter *param, +- double *alpha, Solver::SolutionInfo* si) +-{ +- int l = prob->l; +- double *zeros = new double[l]; +- schar *ones = new schar[l]; +- int i; +- +- int n = (int)(param->nu*prob->l); // # of alpha's at upper bound +- +- for(i=0;i<n;i++) +- alpha[i] = 1; +- if(n<prob->l) +- alpha[n] = param->nu * prob->l - n; +- for(i=n+1;i<l;i++) +- alpha[i] = 0; +- +- for(i=0;i<l;i++) +- { +- zeros[i] = 0; +- ones[i] = 1; +- } +- +- Solver s; +- s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones, +- alpha, 1.0, 1.0, param->eps, si, param->shrinking); +- +- delete[] zeros; +- delete[] ones; +-} +- +-static void solve_epsilon_svr( +- const svm_problem *prob, const svm_parameter *param, +- double *alpha, Solver::SolutionInfo* si) +-{ +- int l = prob->l; +- double *alpha2 = new double[2*l]; +- double *linear_term = new double[2*l]; +- schar *y = new schar[2*l]; +- int i; +- +- for(i=0;i<l;i++) +- { +- alpha2[i] = 0; +- linear_term[i] = param->p - prob->y[i]; +- y[i] = 1; +- +- alpha2[i+l] = 0; +- linear_term[i+l] = param->p + prob->y[i]; +- y[i+l] = -1; +- } +- +- Solver s; +- s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y, +- alpha2, param->C, param->C, param->eps, si, param->shrinking); +- +- double sum_alpha = 0; +- for(i=0;i<l;i++) +- { +- alpha[i] = alpha2[i] - alpha2[i+l]; +- sum_alpha += fabs(alpha[i]); +- } +- info("nu = %f\n",sum_alpha/(param->C*l)); +- +- delete[] alpha2; +- delete[] linear_term; +- delete[] y; +-} +- +-static void solve_nu_svr( +- const svm_problem *prob, const svm_parameter *param, +- double *alpha, Solver::SolutionInfo* si) +-{ +- int l = prob->l; +- double C = param->C; +- double *alpha2 = new double[2*l]; +- double *linear_term = new double[2*l]; +- schar *y = new schar[2*l]; +- int i; +- +- double sum = C * param->nu * l / 2; +- for(i=0;i<l;i++) +- { +- alpha2[i] = alpha2[i+l] = min(sum,C); +- sum -= alpha2[i]; +- +- linear_term[i] = - prob->y[i]; +- y[i] = 1; +- +- linear_term[i+l] = prob->y[i]; +- y[i+l] = -1; +- } +- +- Solver_NU s; +- s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y, +- alpha2, C, C, param->eps, si, param->shrinking); +- +- info("epsilon = %f\n",-si->r); +- +- for(i=0;i<l;i++) +- alpha[i] = alpha2[i] - alpha2[i+l]; +- +- delete[] alpha2; +- delete[] linear_term; +- delete[] y; +-} +- +-// +-// decision_function +-// +-struct decision_function +-{ +- double *alpha; +- double rho; +-}; +- +-decision_function svm_train_one( +- const svm_problem *prob, const svm_parameter *param, +- double Cp, double Cn) +-{ +- double *alpha = Malloc(double,prob->l); +- Solver::SolutionInfo si; +- switch(param->svm_type) +- { +- case C_SVC: +- solve_c_svc(prob,param,alpha,&si,Cp,Cn); +- break; +- case NU_SVC: +- solve_nu_svc(prob,param,alpha,&si); +- break; +- case ONE_CLASS: +- solve_one_class(prob,param,alpha,&si); +- break; +- case EPSILON_SVR: +- solve_epsilon_svr(prob,param,alpha,&si); +- break; +- case NU_SVR: +- solve_nu_svr(prob,param,alpha,&si); +- break; +- } +- +- info("obj = %f, rho = %f\n",si.obj,si.rho); +- +- // output SVs +- +- int nSV = 0; +- int nBSV = 0; +- for(int i=0;i<prob->l;i++) +- { +- if(fabs(alpha[i]) > 0) +- { +- ++nSV; +- if(prob->y[i] > 0) +- { +- if(fabs(alpha[i]) >= si.upper_bound_p) +- ++nBSV; +- } +- else +- { +- if(fabs(alpha[i]) >= si.upper_bound_n) +- ++nBSV; +- } +- } +- } +- +- info("nSV = %d, nBSV = %d\n",nSV,nBSV); +- +- decision_function f; +- f.alpha = alpha; +- f.rho = si.rho; +- return f; +-} +- +-// +-// svm_model +-// +-struct svm_model +-{ +- svm_parameter param; // parameter +- int nr_class; // number of classes, = 2 in regression/one class svm +- int l; // total #SV +- svm_node **SV; // SVs (SV[l]) +- double **sv_coef; // coefficients for SVs in decision functions (sv_coef[n-1][l]) +- double *rho; // constants in decision functions (rho[n*(n-1)/2]) +- double *probA; // pariwise probability information +- double *probB; +- +- // for classification only +- +- int *label; // label of each class (label[n]) +- int *nSV; // number of SVs for each class (nSV[n]) +- // nSV[0] + nSV[1] + ... + nSV[n-1] = l +- // XXX +- int free_sv; // 1 if svm_model is created by svm_load_model +- // 0 if svm_model is created by svm_train +-}; +- +-// Platt's binary SVM Probablistic Output: an improvement from Lin et al. +-void sigmoid_train( +- int l, const double *dec_values, const double *labels, +- double& A, double& B) +-{ +- double prior1=0, prior0 = 0; +- int i; +- +- for (i=0;i<l;i++) +- if (labels[i] > 0) prior1+=1; +- else prior0+=1; +- +- int max_iter=100; // Maximal number of iterations +- double min_step=1e-10; // Minimal step taken in line search +- double sigma=1e-3; // For numerically strict PD of Hessian +- double eps=1e-5; +- double hiTarget=(prior1+1.0)/(prior1+2.0); +- double loTarget=1/(prior0+2.0); +- double *t=Malloc(double,l); +- double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize; +- double newA,newB,newf,d1,d2; +- int iter; +- +- // Initial Point and Initial Fun Value +- A=0.0; B=log((prior0+1.0)/(prior1+1.0)); +- double fval = 0.0; +- +- for (i=0;i<l;i++) +- { +- if (labels[i]>0) t[i]=hiTarget; +- else t[i]=loTarget; +- fApB = dec_values[i]*A+B; +- if (fApB>=0) +- fval += t[i]*fApB + log(1+exp(-fApB)); +- else +- fval += (t[i] - 1)*fApB +log(1+exp(fApB)); +- } +- for (iter=0;iter<max_iter;iter++) +- { +- // Update Gradient and Hessian (use H' = H + sigma I) +- h11=sigma; // numerically ensures strict PD +- h22=sigma; +- h21=0.0;g1=0.0;g2=0.0; +- for (i=0;i<l;i++) +- { +- fApB = dec_values[i]*A+B; +- if (fApB >= 0) +- { +- p=exp(-fApB)/(1.0+exp(-fApB)); +- q=1.0/(1.0+exp(-fApB)); +- } +- else +- { +- p=1.0/(1.0+exp(fApB)); +- q=exp(fApB)/(1.0+exp(fApB)); +- } +- d2=p*q; +- h11+=dec_values[i]*dec_values[i]*d2; +- h22+=d2; +- h21+=dec_values[i]*d2; +- d1=t[i]-p; +- g1+=dec_values[i]*d1; +- g2+=d1; +- } +- +- // Stopping Criteria +- if (fabs(g1)<eps && fabs(g2)<eps) +- break; +- +- // Finding Newton direction: -inv(H') * g +- det=h11*h22-h21*h21; +- dA=-(h22*g1 - h21 * g2) / det; +- dB=-(-h21*g1+ h11 * g2) / det; +- gd=g1*dA+g2*dB; +- +- +- stepsize = 1; // Line Search +- while (stepsize >= min_step) +- { +- newA = A + stepsize * dA; +- newB = B + stepsize * dB; +- +- // New function value +- newf = 0.0; +- for (i=0;i<l;i++) +- { +- fApB = dec_values[i]*newA+newB; +- if (fApB >= 0) +- newf += t[i]*fApB + log(1+exp(-fApB)); +- else +- newf += (t[i] - 1)*fApB +log(1+exp(fApB)); +- } +- // Check sufficient decrease +- if (newf<fval+0.0001*stepsize*gd) +- { +- A=newA;B=newB;fval=newf; +- break; +- } +- else +- stepsize = stepsize / 2.0; +- } +- +- if (stepsize < min_step) +- { +- info("Line search fails in two-class probability estimates\n"); +- break; +- } +- } +- +- if (iter>=max_iter) +- info("Reaching maximal iterations in two-class probability estimates\n"); +- free(t); +-} +- +-double sigmoid_predict(double decision_value, double A, double B) +-{ +- double fApB = decision_value*A+B; +- if (fApB >= 0) +- return exp(-fApB)/(1.0+exp(-fApB)); +- else +- return 1.0/(1+exp(fApB)) ; +-} +- +-// Method 2 from the multiclass_prob paper by Wu, Lin, and Weng +-void multiclass_probability(int k, double **r, double *p) +-{ +- int t,j; +- int iter = 0, max_iter=100; +- double **Q=Malloc(double *,k); +- double *Qp=Malloc(double,k); +- double pQp, eps=0.005/k; +- +- for (t=0;t<k;t++) +- { +- p[t]=1.0/k; // Valid if k = 1 +- Q[t]=Malloc(double,k); +- Q[t][t]=0; +- for (j=0;j<t;j++) +- { +- Q[t][t]+=r[j][t]*r[j][t]; +- Q[t][j]=Q[j][t]; +- } +- for (j=t+1;j<k;j++) +- { +- Q[t][t]+=r[j][t]*r[j][t]; +- Q[t][j]=-r[j][t]*r[t][j]; +- } +- } +- for (iter=0;iter<max_iter;iter++) +- { +- // stopping condition, recalculate QP,pQP for numerical accuracy +- pQp=0; +- for (t=0;t<k;t++) +- { +- Qp[t]=0; +- for (j=0;j<k;j++) +- Qp[t]+=Q[t][j]*p[j]; +- pQp+=p[t]*Qp[t]; +- } +- double max_error=0; +- for (t=0;t<k;t++) +- { +- double error=fabs(Qp[t]-pQp); +- if (error>max_error) +- max_error=error; +- } +- if (max_error<eps) break; +- +- for (t=0;t<k;t++) +- { +- double diff=(-Qp[t]+pQp)/Q[t][t]; +- p[t]+=diff; +- pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff); +- for (j=0;j<k;j++) +- { +- Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff); +- p[j]/=(1+diff); +- } +- } +- } +- if (iter>=max_iter) +- info("Exceeds max_iter in multiclass_prob\n"); +- for(t=0;t<k;t++) free(Q[t]); +- free(Q); +- free(Qp); +-} +- +-// Cross-validation decision values for probability estimates +-void svm_binary_svc_probability( +- const svm_problem *prob, const svm_parameter *param, +- double Cp, double Cn, double& probA, double& probB) +-{ +- int i; +- int nr_fold = 5; +- int *perm = Malloc(int,prob->l); +- double *dec_values = Malloc(double,prob->l); +- +- // random shuffle +- for(i=0;i<prob->l;i++) perm[i]=i; +- for(i=0;i<prob->l;i++) +- { +- int j = i+rand()%(prob->l-i); +- swap(perm[i],perm[j]); +- } +- for(i=0;i<nr_fold;i++) +- { +- int begin = i*prob->l/nr_fold; +- int end = (i+1)*prob->l/nr_fold; +- int j,k; +- struct svm_problem subprob; +- +- subprob.l = prob->l-(end-begin); +- subprob.x = Malloc(struct svm_node*,subprob.l); +- subprob.y = Malloc(double,subprob.l); +- +- k=0; +- for(j=0;j<begin;j++) +- { +- subprob.x[k] = prob->x[perm[j]]; +- subprob.y[k] = prob->y[perm[j]]; +- ++k; +- } +- for(j=end;j<prob->l;j++) +- { +- subprob.x[k] = prob->x[perm[j]]; +- subprob.y[k] = prob->y[perm[j]]; +- ++k; +- } +- int p_count=0,n_count=0; +- for(j=0;j<k;j++) +- if(subprob.y[j]>0) +- p_count++; +- else +- n_count++; +- +- if(p_count==0 && n_count==0) +- for(j=begin;j<end;j++) +- dec_values[perm[j]] = 0; +- else if(p_count > 0 && n_count == 0) +- for(j=begin;j<end;j++) +- dec_values[perm[j]] = 1; +- else if(p_count == 0 && n_count > 0) +- for(j=begin;j<end;j++) +- dec_values[perm[j]] = -1; +- else +- { +- svm_parameter subparam = *param; +- subparam.probability=0; +- subparam.C=1.0; +- subparam.nr_weight=2; +- subparam.weight_label = Malloc(int,2); +- subparam.weight = Malloc(double,2); +- subparam.weight_label[0]=+1; +- subparam.weight_label[1]=-1; +- subparam.weight[0]=Cp; +- subparam.weight[1]=Cn; +- struct svm_model *submodel = svm_train(&subprob,&subparam); +- for(j=begin;j<end;j++) +- { +- svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]])); +- // ensure +1 -1 order; reason not using CV subroutine +- dec_values[perm[j]] *= submodel->label[0]; +- } +- svm_destroy_model(submodel); +- svm_destroy_param(&subparam); +- free(subprob.x); +- free(subprob.y); +- } +- } +- sigmoid_train(prob->l,dec_values,prob->y,probA,probB); +- free(dec_values); +- free(perm); +-} +- +-// Return parameter of a Laplace distribution +-double svm_svr_probability( +- const svm_problem *prob, const svm_parameter *param) +-{ +- int i; +- int nr_fold = 5; +- double *ymv = Malloc(double,prob->l); +- double mae = 0; +- +- svm_parameter newparam = *param; +- newparam.probability = 0; +- svm_cross_validation(prob,&newparam,nr_fold,ymv); +- for(i=0;i<prob->l;i++) +- { +- ymv[i]=prob->y[i]-ymv[i]; +- mae += fabs(ymv[i]); +- } +- mae /= prob->l; +- double std=sqrt(2*mae*mae); +- int count=0; +- mae=0; +- for(i=0;i<prob->l;i++) +- if (fabs(ymv[i]) > 5*std) +- count=count+1; +- else +- mae+=fabs(ymv[i]); +- mae /= (prob->l-count); +- info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae); +- free(ymv); +- return mae; +-} +- +- +-// label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data +-// perm, length l, must be allocated before calling this subroutine +-void svm_group_classes(const svm_problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm) +-{ +- int l = prob->l; +- int max_nr_class = 16; +- int nr_class = 0; +- int *label = Malloc(int,max_nr_class); +- int *count = Malloc(int,max_nr_class); +- int *data_label = Malloc(int,l); +- int i; +- +- for(i=0;i<l;i++) +- { +- int this_label = (int)prob->y[i]; +- int j; +- for(j=0;j<nr_class;j++) +- { +- if(this_label == label[j]) +- { +- ++count[j]; +- break; +- } +- } +- data_label[i] = j; +- if(j == nr_class) +- { +- if(nr_class == max_nr_class) +- { +- max_nr_class *= 2; +- label = (int *)realloc(label,max_nr_class*sizeof(int)); +- count = (int *)realloc(count,max_nr_class*sizeof(int)); +- } +- label[nr_class] = this_label; +- count[nr_class] = 1; +- ++nr_class; +- } +- } +- +- int *start = Malloc(int,nr_class); +- start[0] = 0; +- for(i=1;i<nr_class;i++) +- start[i] = start[i-1]+count[i-1]; +- for(i=0;i<l;i++) +- { +- perm[start[data_label[i]]] = i; +- ++start[data_label[i]]; +- } +- start[0] = 0; +- for(i=1;i<nr_class;i++) +- start[i] = start[i-1]+count[i-1]; +- +- *nr_class_ret = nr_class; +- *label_ret = label; +- *start_ret = start; +- *count_ret = count; +- free(data_label); +-} +- +-// +-// Interface functions +-// +-svm_model *svm_train(const svm_problem *prob, const svm_parameter *param) +-{ +- svm_model *model = Malloc(svm_model,1); +- model->param = *param; +- model->free_sv = 0; // XXX +- +- if(param->svm_type == ONE_CLASS || +- param->svm_type == EPSILON_SVR || +- param->svm_type == NU_SVR) +- { +- // regression or one-class-svm +- model->nr_class = 2; +- model->label = NULL; +- model->nSV = NULL; +- model->probA = NULL; model->probB = NULL; +- model->sv_coef = Malloc(double *,1); +- +- if(param->probability && +- (param->svm_type == EPSILON_SVR || +- param->svm_type == NU_SVR)) +- { +- model->probA = Malloc(double,1); +- model->probA[0] = svm_svr_probability(prob,param); +- } +- +- decision_function f = svm_train_one(prob,param,0,0); +- model->rho = Malloc(double,1); +- model->rho[0] = f.rho; +- +- int nSV = 0; +- int i; +- for(i=0;i<prob->l;i++) +- if(fabs(f.alpha[i]) > 0) ++nSV; +- model->l = nSV; +- model->SV = Malloc(svm_node *,nSV); +- model->sv_coef[0] = Malloc(double,nSV); +- int j = 0; +- for(i=0;i<prob->l;i++) +- if(fabs(f.alpha[i]) > 0) +- { +- model->SV[j] = prob->x[i]; +- model->sv_coef[0][j] = f.alpha[i]; +- ++j; +- } +- +- free(f.alpha); +- } +- else +- { +- // classification +- int l = prob->l; +- int nr_class; +- int *label = NULL; +- int *start = NULL; +- int *count = NULL; +- int *perm = Malloc(int,l); +- +- // group training data of the same class +- svm_group_classes(prob,&nr_class,&label,&start,&count,perm); +- svm_node **x = Malloc(svm_node *,l); +- int i; +- for(i=0;i<l;i++) +- x[i] = prob->x[perm[i]]; +- +- // calculate weighted C +- +- double *weighted_C = Malloc(double, nr_class); +- for(i=0;i<nr_class;i++) +- weighted_C[i] = param->C; +- for(i=0;i<param->nr_weight;i++) +- { +- int j; +- for(j=0;j<nr_class;j++) +- if(param->weight_label[i] == label[j]) +- break; +- if(j == nr_class) +- fprintf(stderr,"warning: class label %d specified in weight is not found\n", param->weight_label[i]); +- else +- weighted_C[j] *= param->weight[i]; +- } +- +- // train k*(k-1)/2 models +- +- bool *nonzero = Malloc(bool,l); +- for(i=0;i<l;i++) +- nonzero[i] = false; +- decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2); +- +- double *probA=NULL,*probB=NULL; +- if (param->probability) +- { +- probA=Malloc(double,nr_class*(nr_class-1)/2); +- probB=Malloc(double,nr_class*(nr_class-1)/2); +- } +- +- int p = 0; +- for(i=0;i<nr_class;i++) +- for(int j=i+1;j<nr_class;j++) +- { +- svm_problem sub_prob; +- int si = start[i], sj = start[j]; +- int ci = count[i], cj = count[j]; +- sub_prob.l = ci+cj; +- sub_prob.x = Malloc(svm_node *,sub_prob.l); +- sub_prob.y = Malloc(double,sub_prob.l); +- int k; +- for(k=0;k<ci;k++) +- { +- sub_prob.x[k] = x[si+k]; +- sub_prob.y[k] = +1; +- } +- for(k=0;k<cj;k++) +- { +- sub_prob.x[ci+k] = x[sj+k]; +- sub_prob.y[ci+k] = -1; +- } +- +- if(param->probability) +- svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p]); +- +- f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]); +- for(k=0;k<ci;k++) +- if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0) +- nonzero[si+k] = true; +- for(k=0;k<cj;k++) +- if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0) +- nonzero[sj+k] = true; +- free(sub_prob.x); +- free(sub_prob.y); +- ++p; +- } +- +- // build output +- +- model->nr_class = nr_class; +- +- model->label = Malloc(int,nr_class); +- for(i=0;i<nr_class;i++) +- model->label[i] = label[i]; +- +- model->rho = Malloc(double,nr_class*(nr_class-1)/2); +- for(i=0;i<nr_class*(nr_class-1)/2;i++) +- model->rho[i] = f[i].rho; +- +- if(param->probability) +- { +- model->probA = Malloc(double,nr_class*(nr_class-1)/2); +- model->probB = Malloc(double,nr_class*(nr_class-1)/2); +- for(i=0;i<nr_class*(nr_class-1)/2;i++) +- { +- model->probA[i] = probA[i]; +- model->probB[i] = probB[i]; +- } +- } +- else +- { +- model->probA=NULL; +- model->probB=NULL; +- } +- +- int total_sv = 0; +- int *nz_count = Malloc(int,nr_class); +- model->nSV = Malloc(int,nr_class); +- for(i=0;i<nr_class;i++) +- { +- int nSV = 0; +- for(int j=0;j<count[i];j++) +- if(nonzero[start[i]+j]) +- { +- ++nSV; +- ++total_sv; +- } +- model->nSV[i] = nSV; +- nz_count[i] = nSV; +- } +- +- info("Total nSV = %d\n",total_sv); +- +- model->l = total_sv; +- model->SV = Malloc(svm_node *,total_sv); +- p = 0; +- for(i=0;i<l;i++) +- if(nonzero[i]) model->SV[p++] = x[i]; +- +- int *nz_start = Malloc(int,nr_class); +- nz_start[0] = 0; +- for(i=1;i<nr_class;i++) +- nz_start[i] = nz_start[i-1]+nz_count[i-1]; +- +- model->sv_coef = Malloc(double *,nr_class-1); +- for(i=0;i<nr_class-1;i++) +- model->sv_coef[i] = Malloc(double,total_sv); +- +- p = 0; +- for(i=0;i<nr_class;i++) +- for(int j=i+1;j<nr_class;j++) +- { +- // classifier (i,j): coefficients with +- // i are in sv_coef[j-1][nz_start[i]...], +- // j are in sv_coef[i][nz_start[j]...] +- +- int si = start[i]; +- int sj = start[j]; +- int ci = count[i]; +- int cj = count[j]; +- +- int q = nz_start[i]; +- int k; +- for(k=0;k<ci;k++) +- if(nonzero[si+k]) +- model->sv_coef[j-1][q++] = f[p].alpha[k]; +- q = nz_start[j]; +- for(k=0;k<cj;k++) +- if(nonzero[sj+k]) +- model->sv_coef[i][q++] = f[p].alpha[ci+k]; +- ++p; +- } +- +- free(label); +- free(probA); +- free(probB); +- free(count); +- free(perm); +- free(start); +- free(x); +- free(weighted_C); +- free(nonzero); +- for(i=0;i<nr_class*(nr_class-1)/2;i++) +- free(f[i].alpha); +- free(f); +- free(nz_count); +- free(nz_start); +- } +- return model; +-} +- +-// Stratified cross validation +-void svm_cross_validation(const svm_problem *prob, const svm_parameter *param, int nr_fold, double *target) +-{ +- int i; +- int *fold_start = Malloc(int,nr_fold+1); +- int l = prob->l; +- int *perm = Malloc(int,l); +- int nr_class; +- +- // stratified cv may not give leave-one-out rate +- // Each class to l folds -> some folds may have zero elements +- if((param->svm_type == C_SVC || +- param->svm_type == NU_SVC) && nr_fold < l) +- { +- int *start = NULL; +- int *label = NULL; +- int *count = NULL; +- svm_group_classes(prob,&nr_class,&label,&start,&count,perm); +- +- // random shuffle and then data grouped by fold using the array perm +- int *fold_count = Malloc(int,nr_fold); +- int c; +- int *index = Malloc(int,l); +- for(i=0;i<l;i++) +- index[i]=perm[i]; +- for (c=0; c<nr_class; c++) +- for(i=0;i<count[c];i++) +- { +- int j = i+rand()%(count[c]-i); +- swap(index[start[c]+j],index[start[c]+i]); +- } +- for(i=0;i<nr_fold;i++) +- { +- fold_count[i] = 0; +- for (c=0; c<nr_class;c++) +- fold_count[i]+=(i+1)*count[c]/nr_fold-i*count[c]/nr_fold; +- } +- fold_start[0]=0; +- for (i=1;i<=nr_fold;i++) +- fold_start[i] = fold_start[i-1]+fold_count[i-1]; +- for (c=0; c<nr_class;c++) +- for(i=0;i<nr_fold;i++) +- { +- int begin = start[c]+i*count[c]/nr_fold; +- int end = start[c]+(i+1)*count[c]/nr_fold; +- for(int j=begin;j<end;j++) +- { +- perm[fold_start[i]] = index[j]; +- fold_start[i]++; +- } +- } +- fold_start[0]=0; +- for (i=1;i<=nr_fold;i++) +- fold_start[i] = fold_start[i-1]+fold_count[i-1]; +- free(start); +- free(label); +- free(count); +- free(index); +- free(fold_count); +- } +- else +- { +- for(i=0;i<l;i++) perm[i]=i; +- for(i=0;i<l;i++) +- { +- int j = i+rand()%(l-i); +- swap(perm[i],perm[j]); +- } +- for(i=0;i<=nr_fold;i++) +- fold_start[i]=i*l/nr_fold; +- } +- +- for(i=0;i<nr_fold;i++) +- { +- int begin = fold_start[i]; +- int end = fold_start[i+1]; +- int j,k; +- struct svm_problem subprob; +- +- subprob.l = l-(end-begin); +- subprob.x = Malloc(struct svm_node*,subprob.l); +- subprob.y = Malloc(double,subprob.l); +- +- k=0; +- for(j=0;j<begin;j++) +- { +- subprob.x[k] = prob->x[perm[j]]; +- subprob.y[k] = prob->y[perm[j]]; +- ++k; +- } +- for(j=end;j<l;j++) +- { +- subprob.x[k] = prob->x[perm[j]]; +- subprob.y[k] = prob->y[perm[j]]; +- ++k; +- } +- struct svm_model *submodel = svm_train(&subprob,param); +- if(param->probability && +- (param->svm_type == C_SVC || param->svm_type == NU_SVC)) +- { +- double *prob_estimates=Malloc(double,svm_get_nr_class(submodel)); +- for(j=begin;j<end;j++) +- target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates); +- free(prob_estimates); +- } +- else +- for(j=begin;j<end;j++) +- target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]); +- svm_destroy_model(submodel); +- free(subprob.x); +- free(subprob.y); +- } +- free(fold_start); +- free(perm); +-} +- +- +-int svm_get_svm_type(const svm_model *model) +-{ +- return model->param.svm_type; +-} +- +-int svm_get_nr_class(const svm_model *model) +-{ +- return model->nr_class; +-} +- +-void svm_get_labels(const svm_model *model, int* label) +-{ +- if (model->label != NULL) +- for(int i=0;i<model->nr_class;i++) +- label[i] = model->label[i]; +-} +- +-double svm_get_svr_probability(const svm_model *model) +-{ +- if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) && +- model->probA!=NULL) +- return model->probA[0]; +- else +- { +- info("Model doesn't contain information for SVR probability inference\n"); +- return 0; +- } +-} +- +-void svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values) +-{ +- if(model->param.svm_type == ONE_CLASS || +- model->param.svm_type == EPSILON_SVR || +- model->param.svm_type == NU_SVR) +- { +- double *sv_coef = model->sv_coef[0]; +- double sum = 0; +- for(int i=0;i<model->l;i++) +- sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param); +- sum -= model->rho[0]; +- *dec_values = sum; +- } +- else +- { +- int i; +- int nr_class = model->nr_class; +- int l = model->l; +- +- double *kvalue = Malloc(double,l); +- for(i=0;i<l;i++) +- kvalue[i] = Kernel::k_function(x,model->SV[i],model->param); +- +- int *start = Malloc(int,nr_class); +- start[0] = 0; +- for(i=1;i<nr_class;i++) +- start[i] = start[i-1]+model->nSV[i-1]; +- +- int p=0; +- int pos=0; +- for(i=0;i<nr_class;i++) +- for(int j=i+1;j<nr_class;j++) +- { +- double sum = 0; +- int si = start[i]; +- int sj = start[j]; +- int ci = model->nSV[i]; +- int cj = model->nSV[j]; +- +- int k; +- double *coef1 = model->sv_coef[j-1]; +- double *coef2 = model->sv_coef[i]; +- for(k=0;k<ci;k++) +- sum += coef1[si+k] * kvalue[si+k]; +- for(k=0;k<cj;k++) +- sum += coef2[sj+k] * kvalue[sj+k]; +- sum -= model->rho[p++]; +- dec_values[pos++] = sum; +- } +- +- free(kvalue); +- free(start); +- } +-} +- +-double svm_predict(const svm_model *model, const svm_node *x) +-{ +- if(model->param.svm_type == ONE_CLASS || +- model->param.svm_type == EPSILON_SVR || +- model->param.svm_type == NU_SVR) +- { +- double res; +- svm_predict_values(model, x, &res); +- +- if(model->param.svm_type == ONE_CLASS) +- return (res>0)?1:-1; +- else +- return res; +- } +- else +- { +- int i; +- int nr_class = model->nr_class; +- double *dec_values = Malloc(double, nr_class*(nr_class-1)/2); +- svm_predict_values(model, x, dec_values); +- +- int *vote = Malloc(int,nr_class); +- for(i=0;i<nr_class;i++) +- vote[i] = 0; +- int pos=0; +- for(i=0;i<nr_class;i++) +- for(int j=i+1;j<nr_class;j++) +- { +- if(dec_values[pos++] > 0) +- ++vote[i]; +- else +- ++vote[j]; +- } +- +- int vote_max_idx = 0; +- for(i=1;i<nr_class;i++) +- if(vote[i] > vote[vote_max_idx]) +- vote_max_idx = i; +- free(vote); +- free(dec_values); +- return model->label[vote_max_idx]; +- } +-} +- +-double svm_predict_probability( +- const svm_model *model, const svm_node *x, double *prob_estimates) +-{ +- if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) && +- model->probA!=NULL && model->probB!=NULL) +- { +- int i; +- int nr_class = model->nr_class; +- double *dec_values = Malloc(double, nr_class*(nr_class-1)/2); +- svm_predict_values(model, x, dec_values); +- +- double min_prob=1e-7; +- double **pairwise_prob=Malloc(double *,nr_class); +- for(i=0;i<nr_class;i++) +- pairwise_prob[i]=Malloc(double,nr_class); +- int k=0; +- for(i=0;i<nr_class;i++) +- for(int j=i+1;j<nr_class;j++) +- { +- pairwise_prob[i][j]=min(max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob); +- pairwise_prob[j][i]=1-pairwise_prob[i][j]; +- k++; +- } +- multiclass_probability(nr_class,pairwise_prob,prob_estimates); +- +- int prob_max_idx = 0; +- for(i=1;i<nr_class;i++) +- if(prob_estimates[i] > prob_estimates[prob_max_idx]) +- prob_max_idx = i; +- for(i=0;i<nr_class;i++) +- free(pairwise_prob[i]); +- free(dec_values); +- free(pairwise_prob); +- return model->label[prob_max_idx]; +- } +- else +- return svm_predict(model, x); +-} +- +-const char *svm_type_table[] = +-{ +- "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL +-}; +- +-const char *kernel_type_table[]= +-{ +- "linear","polynomial","rbf","sigmoid",NULL +-}; +- +-int svm_save_model(const char *model_file_name, const svm_model *model) +-{ +- FILE *fp = fopen(model_file_name,"w"); +- if(fp==NULL) return -1; +- +- const svm_parameter& param = model->param; +- +- fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]); +- fprintf(fp,"kernel_type %s\n", kernel_type_table[param.kernel_type]); +- +- if(param.kernel_type == POLY) +- fprintf(fp,"degree %g\n", param.degree); +- +- if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID) +- fprintf(fp,"gamma %g\n", param.gamma); +- +- if(param.kernel_type == POLY || param.kernel_type == SIGMOID) +- fprintf(fp,"coef0 %g\n", param.coef0); +- +- int nr_class = model->nr_class; +- int l = model->l; +- fprintf(fp, "nr_class %d\n", nr_class); +- fprintf(fp, "total_sv %d\n",l); +- +- { +- fprintf(fp, "rho"); +- for(int i=0;i<nr_class*(nr_class-1)/2;i++) +- fprintf(fp," %g",model->rho[i]); +- fprintf(fp, "\n"); +- } +- +- if(model->label) +- { +- fprintf(fp, "label"); +- for(int i=0;i<nr_class;i++) +- fprintf(fp," %d",model->label[i]); +- fprintf(fp, "\n"); +- } +- +- if(model->probA) // regression has probA only +- { +- fprintf(fp, "probA"); +- for(int i=0;i<nr_class*(nr_class-1)/2;i++) +- fprintf(fp," %g",model->probA[i]); +- fprintf(fp, "\n"); +- } +- if(model->probB) +- { +- fprintf(fp, "probB"); +- for(int i=0;i<nr_class*(nr_class-1)/2;i++) +- fprintf(fp," %g",model->probB[i]); +- fprintf(fp, "\n"); +- } +- +- if(model->nSV) +- { +- fprintf(fp, "nr_sv"); +- for(int i=0;i<nr_class;i++) +- fprintf(fp," %d",model->nSV[i]); +- fprintf(fp, "\n"); +- } +- +- fprintf(fp, "SV\n"); +- const double * const *sv_coef = model->sv_coef; +- const svm_node * const *SV = model->SV; +- +- for(int i=0;i<l;i++) +- { +- for(int j=0;j<nr_class-1;j++) +- fprintf(fp, "%.16g ",sv_coef[j][i]); +- +- const svm_node *p = SV[i]; +- while(p->index != -1) +- { +- fprintf(fp,"%d:%.8g ",p->index,p->value); +- p++; +- } +- fprintf(fp, "\n"); +- } +- +- fclose(fp); +- return 0; +-} +- +-svm_model *svm_load_model(const char *model_file_name) +-{ +- FILE *fp = fopen(model_file_name,"rb"); +- if(fp==NULL) return NULL; +- +- // read parameters +- +- svm_model *model = Malloc(svm_model,1); +- svm_parameter& param = model->param; +- model->rho = NULL; +- model->probA = NULL; +- model->probB = NULL; +- model->label = NULL; +- model->nSV = NULL; +- +- char cmd[81]; +- while(1) +- { +- fscanf(fp,"%80s",cmd); +- +- if(strcmp(cmd,"svm_type")==0) +- { +- fscanf(fp,"%80s",cmd); +- int i; +- for(i=0;svm_type_table[i];i++) +- { +- if(strcmp(svm_type_table[i],cmd)==0) +- { +- param.svm_type=i; +- break; +- } +- } +- if(svm_type_table[i] == NULL) +- { +- fprintf(stderr,"unknown svm type.\n"); +- free(model->rho); +- free(model->label); +- free(model->nSV); +- free(model); +- return NULL; +- } +- } +- else if(strcmp(cmd,"kernel_type")==0) +- { +- fscanf(fp,"%80s",cmd); +- int i; +- for(i=0;kernel_type_table[i];i++) +- { +- if(strcmp(kernel_type_table[i],cmd)==0) +- { +- param.kernel_type=i; +- break; +- } +- } +- if(kernel_type_table[i] == NULL) +- { +- fprintf(stderr,"unknown kernel function.\n"); +- free(model->rho); +- free(model->label); +- free(model->nSV); +- free(model); +- return NULL; +- } +- } +- else if(strcmp(cmd,"degree")==0) +- fscanf(fp,"%lf",¶m.degree); +- else if(strcmp(cmd,"gamma")==0) +- fscanf(fp,"%lf",¶m.gamma); +- else if(strcmp(cmd,"coef0")==0) +- fscanf(fp,"%lf",¶m.coef0); +- else if(strcmp(cmd,"nr_class")==0) +- fscanf(fp,"%d",&model->nr_class); +- else if(strcmp(cmd,"total_sv")==0) +- fscanf(fp,"%d",&model->l); +- else if(strcmp(cmd,"rho")==0) +- { +- int n = model->nr_class * (model->nr_class-1)/2; +- model->rho = Malloc(double,n); +- for(int i=0;i<n;i++) +- fscanf(fp,"%lf",&model->rho[i]); +- } +- else if(strcmp(cmd,"label")==0) +- { +- int n = model->nr_class; +- model->label = Malloc(int,n); +- for(int i=0;i<n;i++) +- fscanf(fp,"%d",&model->label[i]); +- } +- else if(strcmp(cmd,"probA")==0) +- { +- int n = model->nr_class * (model->nr_class-1)/2; +- model->probA = Malloc(double,n); +- for(int i=0;i<n;i++) +- fscanf(fp,"%lf",&model->probA[i]); +- } +- else if(strcmp(cmd,"probB")==0) +- { +- int n = model->nr_class * (model->nr_class-1)/2; +- model->probB = Malloc(double,n); +- for(int i=0;i<n;i++) +- fscanf(fp,"%lf",&model->probB[i]); +- } +- else if(strcmp(cmd,"nr_sv")==0) +- { +- int n = model->nr_class; +- model->nSV = Malloc(int,n); +- for(int i=0;i<n;i++) +- fscanf(fp,"%d",&model->nSV[i]); +- } +- else if(strcmp(cmd,"SV")==0) +- { +- while(1) +- { +- int c = getc(fp); +- if(c==EOF || c=='\n') break; +- } +- break; +- } +- else +- { +- fprintf(stderr,"unknown text in model file\n"); +- free(model->rho); +- free(model->label); +- free(model->nSV); +- free(model); +- return NULL; +- } +- } +- +- // read sv_coef and SV +- +- int elements = 0; +- long pos = ftell(fp); +- +- while(1) +- { +- int c = fgetc(fp); +- switch(c) +- { +- case '\n': +- // count the '-1' element +- case ':': +- ++elements; +- break; +- case EOF: +- goto out; +- default: +- ; +- } +- } +-out: +- fseek(fp,pos,SEEK_SET); +- +- int m = model->nr_class - 1; +- int l = model->l; +- model->sv_coef = Malloc(double *,m); +- int i; +- for(i=0;i<m;i++) +- model->sv_coef[i] = Malloc(double,l); +- model->SV = Malloc(svm_node*,l); +- svm_node *x_space=NULL; +- if(l>0) x_space = Malloc(svm_node,elements); +- +- int j=0; +- for(i=0;i<l;i++) +- { +- model->SV[i] = &x_space[j]; +- for(int k=0;k<m;k++) +- fscanf(fp,"%lf",&model->sv_coef[k][i]); +- while(1) +- { +- int c; +- do { +- c = getc(fp); +- if(c=='\n') goto out2; +- } while(isspace(c)); +- ungetc(c,fp); +- fscanf(fp,"%d:%lf",&(x_space[j].index),&(x_space[j].value)); +- ++j; +- } +-out2: +- x_space[j++].index = -1; +- } +- +- fclose(fp); +- +- model->free_sv = 1; // XXX +- return model; +-} +- +-void svm_destroy_model(svm_model* model) +-{ +- if(model->free_sv && model->l > 0) +- free((void *)(model->SV[0])); +- for(int i=0;i<model->nr_class-1;i++) +- free(model->sv_coef[i]); +- free(model->SV); +- free(model->sv_coef); +- free(model->rho); +- free(model->label); +- free(model->probA); +- free(model->probB); +- free(model->nSV); +- free(model); +-} +- +-void svm_destroy_param(svm_parameter* param) +-{ +- free(param->weight_label); +- free(param->weight); +-} +- +-const char *svm_check_parameter(const svm_problem *prob, const svm_parameter *param) +-{ +- // svm_type +- +- int svm_type = param->svm_type; +- if(svm_type != C_SVC && +- svm_type != NU_SVC && +- svm_type != ONE_CLASS && +- svm_type != EPSILON_SVR && +- svm_type != NU_SVR) +- return "unknown svm type"; +- +- // kernel_type +- +- int kernel_type = param->kernel_type; +- if(kernel_type != LINEAR && +- kernel_type != POLY && +- kernel_type != RBF && +- kernel_type != SIGMOID) +- return "unknown kernel type"; +- +- // cache_size,eps,C,nu,p,shrinking +- +- if(param->cache_size <= 0) +- return "cache_size <= 0"; +- +- if(param->eps <= 0) +- return "eps <= 0"; +- +- if(svm_type == C_SVC || +- svm_type == EPSILON_SVR || +- svm_type == NU_SVR) +- if(param->C <= 0) +- return "C <= 0"; +- +- if(svm_type == NU_SVC || +- svm_type == ONE_CLASS || +- svm_type == NU_SVR) +- if(param->nu < 0 || param->nu > 1) +- return "nu < 0 or nu > 1"; +- +- if(svm_type == EPSILON_SVR) +- if(param->p < 0) +- return "p < 0"; +- +- if(param->shrinking != 0 && +- param->shrinking != 1) +- return "shrinking != 0 and shrinking != 1"; +- +- if(param->probability != 0 && +- param->probability != 1) +- return "probability != 0 and probability != 1"; +- +- if(param->probability == 1 && +- svm_type == ONE_CLASS) +- return "one-class SVM probability output not supported yet"; +- +- +- // check whether nu-svc is feasible +- +- if(svm_type == NU_SVC) +- { +- int l = prob->l; +- int max_nr_class = 16; +- int nr_class = 0; +- int *label = Malloc(int,max_nr_class); +- int *count = Malloc(int,max_nr_class); +- +- int i; +- for(i=0;i<l;i++) +- { +- int this_label = (int)prob->y[i]; +- int j; +- for(j=0;j<nr_class;j++) +- if(this_label == label[j]) +- { +- ++count[j]; +- break; +- } +- if(j == nr_class) +- { +- if(nr_class == max_nr_class) +- { +- max_nr_class *= 2; +- label = (int *)realloc(label,max_nr_class*sizeof(int)); +- count = (int *)realloc(count,max_nr_class*sizeof(int)); +- } +- label[nr_class] = this_label; +- count[nr_class] = 1; +- ++nr_class; +- } +- } +- +- for(i=0;i<nr_class;i++) +- { +- int n1 = count[i]; +- for(int j=i+1;j<nr_class;j++) +- { +- int n2 = count[j]; +- if(param->nu*(n1+n2)/2 > min(n1,n2)) +- { +- free(label); +- free(count); +- return "specified nu is infeasible"; +- } +- } +- } +- free(label); +- free(count); +- } +- +- return NULL; +-} +- +-int svm_check_probability_model(const svm_model *model) +-{ +- return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) && +- model->probA!=NULL && model->probB!=NULL) || +- ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) && +- model->probA!=NULL); +-} +--- a/bio-tools-psort-svmloc/libsvm.h ++++ /dev/null +@@ -1,70 +0,0 @@ +-#ifndef _LIBSVM_H +-#define _LIBSVM_H +- +-#ifdef __cplusplus +-extern "C" { +-#endif +- +-struct svm_node +-{ +- int index; +- double value; +-}; +- +-struct svm_problem +-{ +- int l; +- double *y; +- struct svm_node **x; +-}; +- +-enum { C_SVC, NU_SVC, ONE_CLASS, EPSILON_SVR, NU_SVR }; /* svm_type */ +-enum { LINEAR, POLY, RBF, SIGMOID }; /* kernel_type */ +- +-struct svm_parameter +-{ +- int svm_type; +- int kernel_type; +- double degree; /* for poly */ +- double gamma; /* for poly/rbf/sigmoid */ +- double coef0; /* for poly/sigmoid */ +- +- /* these are for training only */ +- double cache_size; /* in MB */ +- double eps; /* stopping criteria */ +- double C; /* for C_SVC, EPSILON_SVR and NU_SVR */ +- int nr_weight; /* for C_SVC */ +- int *weight_label; /* for C_SVC */ +- double* weight; /* for C_SVC */ +- double nu; /* for NU_SVC, ONE_CLASS, and NU_SVR */ +- double p; /* for EPSILON_SVR */ +- int shrinking; /* use the shrinking heuristics */ +- int probability; /* do probability estimates */ +-}; +- +-struct svm_model *svm_train(const struct svm_problem *prob, const struct svm_parameter *param); +-void svm_cross_validation(const struct svm_problem *prob, const struct svm_parameter *param, int nr_fold, double *target); +- +-int svm_save_model(const char *model_file_name, const struct svm_model *model); +-struct svm_model *svm_load_model(const char *model_file_name); +- +-int svm_get_svm_type(const struct svm_model *model); +-int svm_get_nr_class(const struct svm_model *model); +-void svm_get_labels(const struct svm_model *model, int *label); +-double svm_get_svr_probability(const struct svm_model *model); +- +-void svm_predict_values(const struct svm_model *model, const struct svm_node *x, double* dec_values); +-double svm_predict(const struct svm_model *model, const struct svm_node *x); +-double svm_predict_probability(const struct svm_model *model, const struct svm_node *x, double* prob_estimates); +- +-void svm_destroy_model(struct svm_model *model); +-void svm_destroy_param(struct svm_parameter *param); +- +-const char *svm_check_parameter(const struct svm_problem *prob, const struct svm_parameter *param); +-int svm_check_probability_model(const struct svm_model *model); +- +-#ifdef __cplusplus +-} +-#endif +- +-#endif /* _LIBSVM_H */ diff --git a/debian/patches/series b/debian/patches/series index bcc14e2..64eaee1 100644 --- a/debian/patches/series +++ b/debian/patches/series @@ -1,3 +1,4 @@ biosquid_includes.patch fix_hmmer_interface.patch fix_build_Algorithm-SVM.patch +remove-bundled-party-libsvm.patch -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/psortb.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
