This is an automated email from the git hooks/post-receive script. tille pushed a commit to branch master in repository psortb.
commit 70430347b223592cb39e22f6f215e9159e0a8f39 Author: Andreas Tille <[email protected]> Date: Mon Apr 24 16:28:35 2017 +0200 libsvm is removed from source now --- debian/patches/remove-bundled-party-libsvm.patch | 3116 +--------------------- 1 file changed, 2 insertions(+), 3114 deletions(-) diff --git a/debian/patches/remove-bundled-party-libsvm.patch b/debian/patches/remove-bundled-party-libsvm.patch index 960ba43..461b43b 100644 --- a/debian/patches/remove-bundled-party-libsvm.patch +++ b/debian/patches/remove-bundled-party-libsvm.patch @@ -3,7 +3,7 @@ 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 @@ +@@ -8,5 +8,4 @@ bindings.h lib/Bio/Tools/PSort/SVMLoc.pm sample.model fre_patterns.txt @@ -11,7 +11,7 @@ Last-Update: 2017-04-20 typemap --- a/bio-tools-psort-svmloc/Makefile.PL +++ b/bio-tools-psort-svmloc/Makefile.PL -@@ -9,7 +9,7 @@ +@@ -9,7 +9,7 @@ $CC = 'g++'; @libs = qw/svmloc/; %paths = (); @@ -33,3115 +33,3 @@ Last-Update: 2017-04-20 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 */ -- 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
